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.


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 Sections 2 and 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.

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 h 1 and h 2 , respectively. Also, let the modulus of elasticity of the material in the tensile and compressive zones be E + (z) and E − (z), respectively, while the Poisson's ratios remain the same. 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 1 h and 2 h , respectively. Also, let the modulus of elasticity of the material in the tensile and compressive zones be  () Ez and  () Ez , respectively, while the Poisson's ratios remain the same.  If an exponential function is used to express the function grade of the material, E + (z) and E − (z) may be expressed as where α 1 and α 2 are two grade indexes. E + (z) = E − (z) = E 0 when z = 0, 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 ε x = z/ρ if the plane section assumption holds. Thus, according to Ambartsumyan's bimodular model the bending stress in the tensile and compressive zones, σ + x and σ − x , are also the tensile and compressive principal stress and they are, respectively, and Let the normal resultant at any section be N, thus N = 0 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 If D * 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.

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

Determination of the Neutral Layer
It should be noted here that the two important parameters h 1 and h 2 have still not been determined. From Equations (6) and (7), we may have where α 1 and α 2 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 α 2 1 α 2 2 , Equation (22) may be transformed into a dimensionless form, such that in which H 1 and H 2 are the basic variables and satisfy H 1 + H 2 = 1. 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 e α 1 H 1 and e −α 2 H 2 may be spread with respect to H 1 and H 2 , 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 which is an algebra equation of H 1 and is easily solved either by an analytical method or by a numerical technique once the numerical values of α 1 and α 2 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 H 1 and the position of the unknown neutral layer is determined analytically.

Stress
Let the stress components in the two-dimensional beam problem shown in Figure 1 be σ x , σ z , and τ xz , let the strain components be ε x , ε z , and γ xz , 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 α i (i = 1, 2) 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 ϕ +/− (x, z) depends only on z, that is where f +/− (z) 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 C +/− 1 and C +/− 2 are four undetermined constants. Continuously integrating with respect to z, we obtain where C +/− 3 and C +/− 4 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 C +/− 1 and C +/− 2 .
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.

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 g 1 (z) and g 2 (x) are two undermined functions. Substituting u and w into the third expression in Equation (50), we have where a is a rigid displacement item. Integrating the above expression with respect to z and x, we have, respectively, and where and d 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 w(x, z)| z=0 , 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).

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 f +/− (z) is an unknown function, and it may be determined by satisfying the consistency relation. The strain components expressed in term of f +/− (z) are where       / / / 1 2 3 , and C C C are six undetermined constants; the constant item has been neglected.
Thus, the stress function now has the following form Satisfying the consistency relation for any x gives Continuously integrating with respect to z, we have where C +/− 1 , C +/− 2 and C +/− 3 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 (67) 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 introduced beforehand, we have a simple expression which gives C + 1 = C − 1 due to A + 1 + A − 1 = 0. Second, integrating Equation (73), substituting Equations (70) and (71) into it, and also considering A + 2 and A − 2 introduced beforehand, Equation (73) may be simplified as Combining Equations (74) and (75) will solve C + 1 and C − 1 , and substituting them into Equations (70) and (71), we finally obtain Substituting the determined C +/− 1 , C +/− 2 and C +/− 3 into Equation (67), the stress components are obtained as follows It is easily found that the item Px in σ +/− x 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 u = w = ∂w/∂x = 0 at x = l, z = 0, we have Thus, the final displacements are determined.

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 α 1 = α 2 , the quantities in a bimodular FGM beam regress to the corresponding quantities in a standard FGM beam; when α 1 = α 2 = 0, the regression continues up to the classical problem.

Quantities A Classical Beam A FGM Beam A Bimodular FGM Beam
Modulus of elasticity Moment of inertia I y bh 3 12 Bending stress Shearing stress

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.

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 h/l in the expressions of u * and w * . 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., h/l 1; this makes the magnitude of the u * value much less than the value of w * . Thus, in one-dimensional beam theory the horizontal displacement u * is generally neglected without much error. On the other hand, if h/l 1, also the term (h/l) 2 1 and 0 < µ < 0.5 for common materials; thus, the second term µ(h/l) 2 ζ 2 in w * 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 E + (ζ) = E 0 e α 1 ζ where 0 ≤ ζ ≤ H 1 and E − (ζ) = E 0 e α 2 ζ where −H 2 ≤ ζ ≤ 0, it is easily found from Figure 3 that if the grade indexes α 1 > 0 and α 2 > 0, E + (ζ) > E − (ζ) holds; if α 1 < 0 and α 2 < 0, E + (ζ) < E − (ζ) holds; obviously, α 1 = α 2 = 0 corresponds to the classical problem. Therefore, 13 representative examples concerning the taken values of α 1 and α 2 are selected, including ±0.5, ±1.0, and ±2.0. Some relative parameters, including H 1 and H 2 (from Equation (24)) and a (from Equation (83)), are computed and listed in Table 2.  (24)) and a (from Equation (83)), are computed and listed in Table 2.
(a) (b) From Table 2, it is easily found that for , as the values of  1 and  2 increase, the tensile height decreases while the compressive height increases, which means that the neutral axis is tending downward (see Figures 1 and 2, z axis is down); for , as the absolute values of  1 and  2 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 where h/l is taken as 1/10. For the main three types of cases listed in Table 2  From Table 2, it is easily found that for E + (ζ) > E − (ζ), as the values of α 1 and α 2 increase, the tensile height decreases while the compressive height increases, which means that the neutral axis is tending downward (see Figures 1 and 2, z axis is down); for E + (ζ) < E − (ζ), as the absolute values of α 1 and α 2 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 H 1 and H 2 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., x = l/2 or η = 0.5) of a beam under pure bending is considered, u * 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   Similarly, we may use the midspan stress formulas (η = 0.5) 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 l/h = 10. For the main three cases listed in Table 2, the variation of stresses with ζ(= z/h) are plotted in Figures 6 and 7, in which the shearing stress curve t/p is directly from Equation (88).
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 ζ = 0, which may be easily seen from Figures 4, 6 and 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 E + (ζ) > E − (ζ) and at the edge of the tensile area for E + (ζ) < E − (ζ), while the maximum displacement is equal for E + (ζ) = E − (ζ). Figure 5 shows that, for any point on the neutral layer, the deflection value when E + (ζ) > E − (ζ) is always less than the corresponding value when E + (ζ) < E − (ζ).      , the maximum compressive stress still takes place at the compressive edge 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 E + (ζ) > E − (ζ), 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 E + (ζ) < E − (ζ), 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 E + (ζ) = E − (ζ), 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 M(x) = M for pure bending and M(x) = Px for lateral bending. Via extreme conditions ∂σ +/− x /∂z = 0, we have e α i z/h (1 + z e α i z/h > 0 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 E + (ζ) > E − (ζ), the maximum compressive stress takes place at z = −h/α 2 ; for E + (ζ) < E − (ζ), the maximum tensile stress takes place at z = −h/α 1 . 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 (ζ = 0) and takes zero at the top and bottom of the beam. For E + (ζ) = E − (ζ), the shearing stress in tension and compression is symmetrical with respect to ζ = 0, while for the other two cases the rule does not hold. Moreover, the maximum shearing stress in the case of E + (ζ) > E − (ζ) is less than the maximum stress in the case of E + (ζ) < E − (ζ).

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.
Author Contributions: X.H. and J.S. proposed the studied problem and the corresponding solving method; X.L. and J.D. conducted the theoretical derivation and the computation; X.H. and X.L. wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.