One-Dimensional and Two-Dimensional Analytical Solutions for Functionally Graded Beams with Different Moduli in Tension and Compression.

The material considered in this study not only has a functionally graded characteristic but also exhibits different tensile and compressive moduli of elasticity. One-dimensional and two-dimensional mechanical models for a functionally graded beam with a bimodular effect were established first. By taking the grade function as an exponential expression, the analytical solutions of a bimodular functionally graded beam under pure bending and lateral-force bending were obtained. The regression from a two-dimensional solution to a one-dimensional solution is verified. The physical quantities in a bimodular functionally graded beam are compared with their counterparts in a classical problem and a functionally graded beam without a bimodular effect. The validity of the plane section assumption under pure bending and lateral-force bending is analyzed. Three typical cases that the tensile modulus is greater than, equal to, or less than the compressive modulus are discussed. The result indicates that due to the introduction of the bimodular functionally graded effect of the materials, the maximum tensile and compressive bending stresses may not take place at the bottom and top of the beam. The real location at which the maximum bending stress takes place is determined via the extreme condition for the analytical solution.

1. Introduction

Most materials may exhibit different elastic responses in a state of tension and compression, but these characteristics are often neglected due to the complexity of their analysis. Materials that have apparently different moduli in tension and compression are known as bimodular materials [1], for example, ceramics, graphite, concrete, and some biological materials (nacre, for example [2]). During recent decades, many studies have described useful material models for studying bimodular materials. One is Bert’s model [3] based on the criterion of positive-negative signs of the strains in longitudinal fibers. This model is widely used in laminated composites [4,5,6,7,8]. Another is Ambartsumyan’s bimodular model [9] for isotropic materials, which has attracted the most attention in the engineering community. This model assesses different moduli in terms of tension and compression based on the positive-negative signs of principal stresses, which is especially important for the analysis and design of structures. It is well-known that the cracking direction of a concrete beam is always normal to the direction of principal tensile stresses in the beam. The difficulty in applying Ambartsumyan’s bimodular model is that the stress state of a point must be known in advance. However, with the exception of some fundamental problems, we must resort to finite element analysis to acquire the states of the stresses in a structure [10,11,12,13,14].

In addition to the bimodular effect in materials, it is also interesting to consider the functionally graded characteristic of materials. Functionally graded materials (FGMs) possess properties that vary gradually with the location within the material. The use of FGMs has many advantages in aerospace, automotive, and biomedical applications. There are many approximations that may be used to model the variation of properties in FGMs. One is the exponential variation, where the elastic constants vary according to the form of the exponential function. Many researchers have found this functional form to be convenient in solving elasticity problems. Sankar [15] obtained an elasticity solution for a functionally graded beam subjected to transverse loads in which the Young’s modulus is assumed to vary exponentially through the thickness and the Poisson ratio is held constant. Sankar and co-workers studied the relative issues of functionally graded beams, including thermal stresses [16], a sandwich beam with a functionally graded core [17], and a combined Fourier series–Galerkin method [18]. Without specifying the gradient variations of a material property, Zhong and co-workers presented a general solution of a functionally graded beam by the Airy stress function method [19] and a displacement function approach [20]. Daouadji et al. [21] employed the stress function approach to study the problem of a functionally graded cantilever beam subjected to a linearly distributed load, in which the Young’s modulus along the thickness direction varies with power-law functions or with exponential functions. Considering that there are many research works in this field, we do not review them in detail.

Recently, analytical studies of bimodular beams and plates have been performed. Among these works, the determination of the unknown neutral layer is a key issue because it opens up the possibility for the establishment of a mechanical model based on a subarea in tension and compression. Under the assumption that shearing stresses have no contribution to the neutral axis, Yao and Ye [22] obtained a one-dimensional analytical solution of a bimodular shallow beam. He et al. adopted the stress function method to find the elasticity solution of a bimodular straight beam [23] and curved beams [24]. Later, the classical Kirchhoff hypothesis was used to assess the existence of the elastic neutral layers of a thin plate during bending with a small deflection [25]. Consequently, a series of analytical solutions of plates is derived in rectangular and polar coordinate systems. More recently, He et al. [26] presented an elasticity solution of a bimodular FGM beam under uniformly distributed loads and discussed several concrete numerical examples. However, some basic problems are still unclear, which include the consistency between a one-dimensional solution and a two-dimensional solution, the validity of the plane section assumption, the corresponding relation among a classical beam, a standard FGM beam, and a bimodular FGM beam, as well as the bimodular effect on stresses and deformations in a general sense.

In this study, we will adopt a bimodular FGM beam theory to derive the one-dimensional and two-dimensional solutions. Theoretically speaking, any FGM beams may be suitable for this theory provided that the bimodular effect in tension and compression needs to be emphasized for a refined analysis; or, in other words, a certain constituent that forms functionally graded materials presents a relatively obvious bimodular effect which can not be ignored otherwise it will introduce much error into the analysis. The article is organized as follows. The corresponding analytical solutions under pure bending and lateral-force bending will be obtained in Section 2 and Section 3, respectively. Specifically, a perturbation method is adopted to solve the transcendental equation for the determination of the unknown neutral layer. The validity of the plane section assumption is discussed and some important physical quantities among a classical beam, a standard FGM beam, and a bimodular FGM beam are compared in Section 4. Besides this, without specifying the real magnitude of the external load and the geometrical dimension of the beam, the bimodular effect on the stress and deformation in a general sense will be investigated in Section 4. Some important conclusions and subsequent studies are given in the concluding remarks.

2. Functionally Graded Beams under Pure Bending

2.1. One-Dimensional Solution

2.1.1. Bending Stress

A bimodular functionally graded beam with a rectangular section dimension of h × b is subjected to a bending moment M alone as shown in Figure 1. This causes a bending of the beam in the plane coordinate system xoz. Note that due to the introduction of the bimodular effect in tension and compression as well as the functionally graded characteristic of the material, the neutral layer of the beam generally does not locate on the half height of the section. The x axis is established on the unknown neutral layer as shown in Figure 1. It is obvious that the zone below the neutral layer is in tension while the zone up the layer is in compression. Let the tensile and compressive section heights of the beam be and , respectively. Also, let the modulus of elasticity of the material in the tensile and compressive zones be and , respectively, while the Poisson’s ratios remain the same.

Figure 1

Scheme of a bimodular functionally graded material (FGM) beam under pure bending.

If an exponential function is used to express the function grade of the material, and may be expressed as where and are two grade indexes. when , that is, at the neutral layer the tensile modulus is equal to the compressive one. Let the curvature radius of the neutral layer be ; then, the bending strain along the x axis in the whole beam will be the uniform expression if the plane section assumption holds. Thus, according to Ambartsumyan’s bimodular model the bending stress in the tensile and compressive zones, and , are also the tensile and compressive principal stress and they are, respectively, and

Let the normal resultant at any section be N, thus yields

Substituting Equations (2) and (3) into Equation (4), we have

If we let Equation (4) will lead to the following relation which is used for the determination of the unknown neutral layer later.

Similarly, the bending moment at any section will give

Substituting Equations (2) and (3) into Equation (8), we have

If we let

Equation (8) will yield

If is introduced to denote the flexural stiffness of the bimodular functionally graded beam, that is, the deformation of the beam will follow the familiar form

Substituting the relation (11) into Equations (2) and (3), we obtain the one-dimensional solution of the bending stress in the tensile and compressive zones, respectively, and

It should be noted here that due to this being the pure bending case, only the bending stress may be obtained and the shearing stress can be derived in the lateral-force bending case, which will be discussed in Section 3.

2.1.2. Deflection Curve

Let the vertical displacement of any point on the neutral layer be w; then, Equation (11) may be expressed in terms of the second-order derivative of w to x as follows

Integrating twice with respect to x, we have where c and d are two undetermined constants. If a simply-supported beam is considered, the boundary conditions give where l is the span length of the beam. Thus, the deflection curve of the neutral axis is

If a cantilever beam with the right end fixed is considered, as shown in Figure 1, the displacement restriction is and the deflection curve of the neutral axis will be

2.1.3. Determination of the Neutral Layer

It should be noted here that the two important parameters and have still not been determined. From Equations (6) and (7), we may have where and are two indexes concerning the grade function as indicated above. If we introduce the following dimensionless variables and also multiply the two ends of the equation by , Equation (22) may be transformed into a dimensionless form, such that in which and are the basic variables and satisfy . It is a transcendental equation and is hard to solve analytically to some extent due to the existence of an exponential function. Next, we will adopt the perturbation idea to solve the transcendental equation.

The exponential items and may be spread with respect to and , respectively,

If the linear approximation is adopted, such that

substituting it into Equation (24) will yield which is exactly the solution of a classical problem without considering the functionally graded property and bimodular effect of the material. We call it the first-order approximation solution of the problem. If the second-order approximation is adopted, such that substituting it into Equation (24) and considering yields which is an algebra equation of and is easily solved either by an analytical method or by a numerical technique once the numerical values of and are known. The solution of Equation (29) may be called the second-order approximation solution. Similarly, if more items in Equation (25) are taken, we will obtain a high-order approximation solution according to the procedure indicated above. Thus, based on the perturbation idea, the transcendental equation may be gradually transformed into a nonlinear algebra equation of and the position of the unknown neutral layer is determined analytically.

2.2. Two-Dimensional Solution

2.2.1. Stress

Let the stress components in the two-dimensional beam problem shown in Figure 1 be , , and , let the strain components be , , and , and also let the displacement components in the same problem be u and w. Then, in the differential equation of equilibrium in which the body forces are neglected, the geometrical relation as well as the consistency equation are the same as those in the classical problem, and they are, respectively, and

The physical equation gives

After considering the different moduli in tension and compression as well as the functional grade of the material, the physical equation may take the following form where a superscript “” denotes a tensile (compressive) quantity and correspond to the cases of tension and compression, respectively. Equation (33) is in essence two sets of equations concerning tension and compression.

Next, the stress function method will be adopted to obtain the solution of this two-dimensional problem. Due to pure bending, here we still consider that the stress function depends only on z, that is where is an unknown function and “” still denotes a tensile (compressive) quantity. According to the relation between the stress function and the stress components,

Equation (33) may be changed as

Letting Equation (36) satisfy the consistency relation, we obtain

Integrating twice with respect to z, we have where and are four undetermined constants. Continuously integrating with respect to z, we obtain where and are four undetermined constants and may be neglected. The stress function is simplified as

Correspondingly, the stress expressions are

Next, we will use the boundary conditions as well as the continuity condition of stress to determine the four unknown constants and .

First, the continuity conditions of the stresses on the neutral layer give

According to Equation (41), it is easily found that the last two conditions are surely satisfied and the first condition yields

The stress boundary conditions on the two main sides of the beam are, respectively, which are surely satisfied due to pure bending. At the left end of the beam, the application of Saint-Venant’s Principle gives

It is easily found that the last condition is satisfied and the first two conditions will yield, respectively, and

Considering the Equations (6), (7), and (10), we solve

Thus, the final stress components are which is the same as the one-dimensional solution obtained in Section 2.1.1.

2.2.2. Displacement

After the determination of the stress components, the combination of the physical equations and the geometrical equations will give

Integrating the first two expressions with respect to x and z, we have, respectively, and where and are two undermined functions. Substituting and into the third expression in Equation (50), we have where is a rigid displacement item. Integrating the above expression with respect to z and x, we have, respectively, and where and are still rigid displacement items. Now, the displacement may be expressed as and

If we consider here a simply-supported beam, the corresponding boundary conditions give where l is the span length of the beam. Thus, the last displacement components are

The deflection curve of the neutral layer may be obtained by , which is the same as the one-dimensional solution, i.e., Equation (19). If a cantilever beam with the right end fixed is considered, the restriction conditions yield the last displacement components will be

Similarly, the deflection curve of the neutral layer is consistent with the result presented in Equation (21).

3. Bimodular Functionally Graded Beams under Latera-Force Bending

Let us consider the lateral-force bending problem of a bimodular functionally graded beam, as shown in Figure 2, in which the left end of the beam is subjected to the action of a concentrated force P and the right end is fixed. Due to the combined action of the bending moment and the shearing force, any point in the beam is in diagonal tension or diagonal compression; so, it is very difficult to determine the position and shape of the unknown neutral layer if the constitutive law defined in the principal stress direction is strictly followed. For this purpose, an important assumption that shearing stresses have no contribution to the neutral axis [22] is used to establish the simplified mechanical model. In the light of the assumption, the beam will deflect and develop a so-called tensile zone and compressive zone under the external load. The tension and compression of any point in the beam depend only on the direction of the bending stress and are independent of the shearing stress. Thus, similar to the case of pure bending shown in Figure 1, the mechanical model based on a subarea in tension and compression is still established in the case of lateral-force bending as shown in Figure 2. The basic equations of the problem are the same as those in Section 2.2.1, that is, Equations (30)–(33). According to the loading conditions, the stress function may be assumed to be where is an unknown function, and it may be determined by satisfying the consistency relation. The strain components expressed in term of are

Figure 2

Scheme of a bimodular FGM beam under lateral-force bending.

Satisfying the consistency relation for any x gives

Continuously integrating with respect to z, we have where are six undetermined constants; the constant item has been neglected. Thus, the stress function now has the following form

The stress components expressed in terms of the undetermined constants are

The continuity conditions of the stresses on the neutral layer under lateral-force bending are the same as those under pure bending; thus, applying Equation (42) yields and

Similarly, the stress boundary conditions on the two main sides of the beam are the same as those in Equation (44). Satisfying the conditions in the tensile and compressive zones yields, respectively, and

At the left end of the beam, the application of Saint-Venant’s Principle gives

It is easily found that the first two conditions are satisfied and the last condition gives

Equations (69), (70), (71), and (73) may be used for the solution of . First, substituting Equations (70) and (71) into Equation (69) and also considering and introduced beforehand, we have a simple expression which gives due to . Second, integrating Equation (73), substituting Equations (70) and (71) into it, and also considering and introduced beforehand, Equation (73) may be simplified as

Combining Equations (74) and (75) will solve and , and substituting them into Equations (70) and (71), we finally obtain

Substituting the determined into Equation (67), the stress components are obtained as follows

It is easily found that the item in is exactly the magnitude of the bending moment, which is consistent with Equations (14) and (15).

By use of the physical equation and the geometrical equation, the displacement components may be determined as where a, d, and c are the items concerning rigid displacement. Using the boundary condition at we have

Thus, the final displacements are determined.

4. Results and Discussions

4.1. Comparision among Three Types of Beam

As indicated before, the material considered in this study not only has a functionally graded characteristic but also presents different mechanical properties in tension and compression. It is valuable to compare physical quantities in a bimodular FGM beam and a standard FGM beam (without bimodular effect) with their counterparts in a classical problem. We should note that in a classical problem, there is no variation of material properties along the thickness direction; thus, the relevant integrals are usually done over the whole section height. The comparisons among the three types of beams are listed in Table 1. It is easily found that when the grade indexes , the quantities in a bimodular FGM beam regress to the corresponding quantities in a standard FGM beam; when , the regression continues up to the classical problem.

Table 1

Comparisons among a classical beam, an FGM beam, and a bimodular FGM beam.

Quantities	A Classical Beam	A FGM Beam	A Bimodular FGM Beam
Modulus of elasticity
E	E=Const.	E(z)=E_0e^αz/h	E^+(z)=E_0e^α_1z/h,E^−(z)=E_0e^α_2z/h
Moment of inertia
I_y	bh^312	∫_Ae^αz/hz^2dA	∫0h_1e^α_1z/hz^2bdz+∫−h_20e^α_2z/hz^2bdz=b(A2++A2−)
Bending stiffness
D	EI_y	E_0∫_Ae^αz/hz^2dA	E_0b(A2++A2−)
Curvature
1ρ	MEI_y	ME_0∫_Ae^αz/hz^2dA	ME_0b(A2++A2−)
Bending stress
σ_x	MI_yz	Mz∫_Ae^αz/hz^2dA	σx+=Mb(A2++A2−)ze^α_1z/hσx−=Mb(A2++A2−)ze^α_2z/h
Static moment when computing shearing stress
S_y	b2(h^24−z^2)	∫_Ae^αz/hzdA	{S^+=∫zh_1e^a_1z/hzbdz      for 0≤z≤h_1=b[(h^2a12−hza_1)e^a_1z/h−(h^2a12−hh_1a_1)e^a_1h_1/h]S^−=∫z−h_1e^a_2z/hzbdz      for −h_2≤z≤0=b[(h^2a22−hza_2)e^a_2z/h−(h^2a22+hh_2a_2)e^a_2h_2/h]
Shearing stress
τ_xz	PS_yI_yb	P∫_Ae^αz/hzdAb∫_Ae^αz/hz^2dA	{τzx+=      for 0≤z≤h_1Pb(A2++A2−)[(h^2a12−hza_1)e^a_1z/h−(h^2a12−hh_1a_1)e^a_1h_1/h],τzx+=      for −h_2≤z≤0Pb(A2++A2−)[(h^2a22−hza_2)e^a_2z/h−(h^2a22−hh_2a_2)e^a_2h_2/h],

4.2. Plane Section Assumption

For the pure bending problem, the rotation of a vertical element of the cross section, , may be obtained from Equation (61),

It is obvious that the rotation is not dependent on z, which shows that for the pure bending problem, the plane section assumption is surely satisfied. However, for the lateral-force bending problem, the rotation may be obtained from Equation (78), respectively, for the tensile area and for the compressive area

It is readily seen that the rotation is now the function of z. This means that on any cross section, a vertical element under bending will deviate from the original vertical direction and the deviated value varies with the distance from the neutral layer, i.e., z. Consequently, for the lateral-force bending problem the plane section assumption no longer holds. Moreover, unlike the pure bending problem, the rotation will not continuously develop at the neutral layer due to the difference in tension and compression.

4.3. Bimodular Effect on Stress and Displacement

The bimodular effect on stress and displacement may be analyzed by the use of the analytical results obtained. To avoid the inconvenience introduced by the dimension of physical quantities, besides Equation (23), we adopt the following dimensionless manner:

The two-dimensional solution for stress and displacement under pure bending, i.e., Equations (49) and (61), may be changed as, respectively, and

The above dimensionless displacement is helpful for analyzing the approximation degree from a two-dimensional solution to a one-dimensional one. We note that there exists a common factor in the expressions of and . If a typical shallow beam is considered here, the ratio of the beam height to span length will be much less than 1, i.e., ; this makes the magnitude of the value much less than the value of . Thus, in one-dimensional beam theory the horizontal displacement is generally neglected without much error. On the other hand, if , also the term and for common materials; thus, the second term in may be neglected comparing to other items. This yields which is exactly the dimensionless one-dimensional solution for displacement.

Similarly, the two-dimensional solution for stress under lateral-force bending, i.e., Equation (77), may be changed as

Considering the characteristics of the grade function where and where , it is easily found from Figure 3 that if the grade indexes and , holds; if and , holds; obviously, corresponds to the classical problem. Therefore, 13 representative examples concerning the taken values of and are selected, including , , and . Some relative parameters, including and (from Equation (24)) and a (from Equation (83)), are computed and listed in Table 2.

Figure 3

Variation of with the thickness direction. (a) and ; (b) and .

Table 2

Numerical values of , , and a in different cases.

Cases	Groups	α_1	α_2	H_1	H_2	a
E^+(ζ)>E^−(ζ)	(a)	1.0	2.0	0.3725	0.6275	0.0560
	(b)	2.0	1.0	0.3859	0.6141	0.0836
	(c)	1.0	1.0	0.4180	0.5820	0.0762
	(d)	1.0	0.5	0.4399	0.5601	0.0872
	(e)	0.5	0.5	0.4585	0.5415	0.0815
E^+(ζ)=E^−(ζ)	(f)	0	0	1/2	1/2	1/12
E^+(ζ)<E^−(ζ)	(g)	−2.0	−1.0	0.6275	0.3725	0.0560
	(h)	−1.0	−2.0	0.6141	0.3859	0.0836
	(i)	−1.0	−1.0	0.5820	0.4180	0.0762
	(j)	−1.0	−0.5	0.5638	0.4362	0.0720
	(k)	−0.5	−0.5	0.5415	0.4585	0.0815

From Table 2, it is easily found that for , as the values of and increase, the tensile height decreases while the compressive height increases, which means that the neutral axis is tending downward (see Figure 1 and Figure 2, z axis is down); for , as the absolute values of and increase, the tensile height increases while the compressive height decreases, which means the neutral axis is tending upward. Besides this, we also note another interesting phenomenon, which is that due to the characteristic of an exponential function, the heights in tension and compression and are exactly reversed in some cases, including groups (a) and (g), (b) and (h), (c) and (i), and (e) and (k). For the values of a, they are the same as in the combinations above.

If the midspan displacement (i.e., or ) of a beam under pure bending is considered, in Equation (85) may be changed as where h/l is taken as 1/10. For the main three types of cases listed in Table 2, i.e., the representative groups (d), (f), and (j), the varying curves of with as well as the deflection curve of the neutral layer (, see in Equation (86)) with are plotted in Figure 4 and Figure 5, respectively, in which the solid lines correspond to the case of , the dashed lines correspond to the case of , and the dotted lines correspond to the case of .

Figure 4

Variation of displacement at midspan () with in three cases.

Figure 5

Variation of deflection of the neutral layer () with in three cases.

Similarly, we may use the midspan stress formulas () of a beam under lateral-force bending to analyze the bimodular effect on the bending stress and shearing stress. Thus, Equation (87) is changed as where . For the main three cases listed in Table 2, the variation of stresses with are plotted in Figure 6 and Figure 7, in which the shearing stress curve is directly from Equation (88).

Figure 6

Variation of bending stress at midspan () with in three cases.

Figure 7

Variation of shearing stress with in three cases.

We should note such a fact that since the neutral layer is established on the x axis beforehand, the dividing line between tension and compression is always on , which may be easily seen from Figure 4, Figure 6 and Figure 7. Figure 4 shows that the horizontal displacement varies in a linear relation along the direction of the beam thickness as indicated in Equation (89). The maximum horizontal displacement takes place at the edge of the compressive area for and at the edge of the tensile area for , while the maximum displacement is equal for . Figure 5 shows that, for any point on the neutral layer, the deflection value when is always less than the corresponding value when .

Figure 6 presents a typical exponent relation of bending stress varying along the direction of the beam thickness. Due to the variation of elastic modulus with the thickness direction, the location at which the maximum stress takes place may be changed. For , the maximum tensile stress still takes place at the tensile edge of the beam while the maximum compressive stress will take place on a certain level between the compressive edge and the neutral layer; for , the maximum compressive stress still takes place at the compressive edge of the beam while the maximum tensile stress will take place on a certain level between the tensile edge and the neutral layer; for , the maximum tensile and compressive stress are equal and take place at the tensile and compressive edges of the beam, respectively, as we expected. This conclusion may be proved by the use of the extreme condition for an analytical solution of bending stress. We take the first-order derivative of bending stress with respect to the thickness direction, z, such that, where for pure bending and for lateral bending. Via extreme conditions , we have permanently holds true, we have which determines the location at which the maximum tensile or compression stress takes place. By referring to Figure 3, it is obvious that for , the maximum compressive stress takes place at ; for , the maximum tensile stress takes place at . This phenomenon is quite different from the classical problem.

For the three cases of different moduli in tension and compression, Figure 7 uniformly indicates that the maximum shearing stress takes place at the neutral layer () and takes zero at the top and bottom of the beam. For , the shearing stress in tension and compression is symmetrical with respect to , while for the other two cases the rule does not hold. Moreover, the maximum shearing stress in the case of is less than the maximum stress in the case of .

5. Concluding Remarks

In this study, one-dimensional and two-dimensional mechanical models for a functionally graded beam with different moduli in tension and compression were established. The corresponding analytical solutions under pure bending and lateral-force bending were obtained. The following three conclusions can be drawn.

(1) The mechanical models established on the one-dimensional and two-dimensional theory are consistent; the two-dimensional solution may regress to the corresponding one-dimensional solution.

(2) For pure bending problems, the plane section assumption still holds for a bimodular functionally graded beam; for lateral-force bending problems, the plane section assumption holds only in the case of a shallow beam.

(3) The introduction of the bimodular effect and functionally graded characteristic of materials will change the stress and deformation of the structure to some extent. Specifically, the maximum bending stress may take place at a certain level between the neutral layer and edge fibers of the beam, which should be given more attention in the analysis and design of similar structures.

The material considered in this study not only has a functionally graded characteristic but also exhibits different tensile and compressive moduli of elasticity, which further complicates the analysis of similar structures made from these materials. It will be worthwhile considering the plate model adopting classical plate theory for laminate (or higher order theory) to discretize the material properties along the direction of the plate thickness (or here along the beam height).

Moreover, since beams, plates, and shells can all be attributed to, from the point of view of loading and deformation, bending elements under external loads, this work may be extended to the static and dynamic responses of functionally graded beams [27], of functionally graded plates [28], as well as of functionally graded shells [29], in which the bimodular effect of the materials will be incorporated. At the same time, this work may also be extended to an investigation on the existing capabilities and limitations in numerical modeling of fracture problems in functionally graded materials by means of the well-known finite element code ABAQUS [30]. We will study these interesting issues in the future.