- Email: [email protected]
A traction-recovery method for evaluating boundary stresses on thermal elasticity problems of FGMs
Engineering Analysis with Boundary Elements 61 (2015) 226–231
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A traction-recovery method for evaluating boundary stresses on thermal elasticity problems of FGMs Jian Liu, Hai-Feng Peng, Xiao-Wei Gao, Miao Cui n State Key Laboratory of Structural Analysis for Industrial Equipment, School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China
art ic l e i nf o
a b s t r a c t
Article history: Received 13 November 2014 Received in revised form 9 June 2015 Accepted 22 July 2015 Available online 27 August 2015
A new approach is presented to calculate boundary stresses in thermal stress analysis of structures consisting of functionally grades materials (FGMs) based on the traction-recovery method. In this approach, the in-plane strains are calculated first using the computed nodal displacements by simply differentiating shape functions at the point of interest, and then the boundary stresses are recovered by Hooke's law together with the known tractions on the boundary. This approach has the advantage of without need to evaluate strongly singular boundary integrals. With the comparison to the FEM software ANSYS, two numerical examples for plane stress and 3D problems are presented to verify the correctness of the proposed method in evaluating boundary thermal-stresses of FGMs. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Traction-recovery method Boundary stresses Thermo-elasticity Functionally graded material Boundary element method
1. Introduction In recent decades, the boundary element method (BEM) has attracted growing interest and appears to be an effective alternative to the finite element method (FEM) for certain engineering problems [1–4], thanks to the reduction of dimensionality, the superior in solving infinite or semi-infinite field problems, and the highly accurate solution on the boundary, etc. However, one of the issues that have still received a great deal of attention is how to accurately determine the boundary stresses for some complicated thermal-stress problems such as non-homogenous and nonisotropic problems [5–8]. It has been found that the determinations of the stresses inside the domain and on the boundary are quite different [2–5]. The stresses inside the domain can be directly evaluated using the stress boundary integral equations. However, the stresses on the boundary are not easy to be evaluated in the same way, due to the hyper-singularities existing in the stress boundary integrals. Therefore, for many years the computation of boundary stresses is an important research issue in the BEM community. Two main methods are used to evaluate the boundary stresses or related high order singular boundary integrals. The first one is the direct method [6–13], which follows the same procedure as computing the internal stresses. Although this method is
n
Corresponding author. E-mail address: [email protected] (M. Cui).
http://dx.doi.org/10.1016/j.enganabound.2015.07.016 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
considered to be more accurate, it does require a substantially much more computational time due to the cost of numerically evaluating all boundary integrals [14]. What's worse, evaluating the high order singular boundary integrals is really an enormous challenge for the BEM researchers. The second method, the most popular one, to evaluate boundary stresses in BEM is the traction-recovery method (TRM) [3,5,15,16], which computes the boundary stresses from the computed displacements and tractions by differentiating shape functions of the boundary element under consideration. Due to the differentiation of shape functions with respect to intrinsic coordinates, the accuracy using TRM is lower than the direct method. On the other hand, due to the use of the precisely computed traction on the surface to recover normal direction related stresses, the accuracy of TRM is higher than the finite element method (FEM). For this reason, this method has been widely used. For the traction-recovery method, Gao and Daivies [3] gives a complete set of formulations and Fortran codes for isotropic materials without consideration of the thermal effect. However, up to date, no detailed formulations are published for the thermal-stress analysis of FGMs. This paper tries to close this gap for the first time.
2. Strain–stress relationship in thermo-elasticity for FGMs Considering the volumetric expansion caused by temperature change over a functionally graded material, the relationship between the total strain εij and the stress σ ij in thermo-elasticity
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
can be expressed as [17,18] 1 ν σ ij δij σ kk þδij kðxÞθ εij ¼ 2μðxÞ 1þν
ð1Þ
where the shear modulus μ and thermal expansion coefficient k are the functions of the coordinates x, Poisson's ratio ν is constant, and θ is the temperature change. Eq. (1) can be converted to νδij εkk ~ ð2Þ δij kðxÞθ σ ij ¼ 2μðxÞ εij þ 1 2ν where 2μðxÞð1 þ νÞkðxÞ ~ kðxÞ ¼ 1 2ν
ð3Þ
By eliminating the local normal strain ε33 in Eq. (2), we can obtain the relationship between the stresses and strains in the tangential plane as follows: 2μðxÞ ν ~ ~ ðσ 33 þ kðxÞθÞ ðε11 þ νε22 Þ þ kðxÞθ 1ν 1ν 2μðxÞ ν ~ ~ ðσ 33 þ kðxÞθÞ σ 22 ¼ ðε22 þ νε11 Þ þ kðxÞθ 1ν 1ν σ 12 ¼ 2μðxÞε12 σ 11 ¼
the shape functions of the boundary element [3]. For 2-noded linear boundary element, N l can be expressed as: 1 N1 ¼ ð1 ξÞ; 2
1 N2 ¼ ð1 þ ξÞ 2
where ξ is the local coordinate ½ 1; þ 1. Using the above Eqs. (5)–(10), the local tangential strains ε0IJ can be obtained. Substituting the results into Eq. (4) yields ∂uj 2μðxÞ ∂ξK ∂ξK ν ~ ~ L t þ kðxÞθ L þ ν L þ kðxÞθ σ 011 ¼ 1j 2j 1 ν ∂x01 ∂x01 ∂ξK 1 ν 3j j ∂uj 2μðxÞ ∂ξK ∂ξ ν ~ ~ L3j t j þ kðxÞθ σ 022 ¼ L2j þ ν K0 L1j þ kðxÞθ 0 1 ν ∂x2 ∂x1 ∂ξK 1 ν ∂uj ∂ξK ∂ξ σ 012 ¼ μðxÞ L2j þ K0 L1j ð12Þ 0 ∂x1 ∂x2 ∂ξK According to the equilibrium equation, the other local stress components related to the normal direction can be recovered from the tractions. σ 023 ¼ t 02 ¼ L2j t j
ð4Þ
σ 013 ¼ t 01 ¼ L1j t j
where the upper-case subscripts I, J range from 1 to 2 and indicate the values refer to the local axis system. The local derivatives of displacement can be obtained by chain rule ∂u0I ∂u0I ∂ξK ¼ ∂x0J ∂ξK ∂x0J
ð6Þ
σ mn ¼ Lrm Lsn σ 0rs
σ mn ¼ Amnjl ulj þBmnj t j C mn θ where the coefficients Amnjl , Bmnj and C mn are 1 ∂ξK ∂ξ L1m L1n Amnjl ¼ 2μðxÞ L1j þ ν K0 L2j 0 1ν ∂x1 ∂x2
∂ξK ∂ξK þ L2m L2n L2j þ ν 0 L1j : ∂x02 ∂x1 1 ∂ξK ∂ξK ∂N l L þ L þ ðL1m L2n þ L2m L1n Þ 2j 1j 2 ∂x01 ∂x02 ∂ξK Bmnj ¼ ðL3m L1n þ L1m L3n ÞL1j þ ðL2m L3n þ L3m L2n ÞL2j ν 1 2ν þ δmn þ L3m L3n L3j 1ν 1ν
∂ξ1 1 ∂ξ1 cos θ ; ; ¼ ¼ ∂x01 jm1 j ∂x02 jm1 j sin θ
C mn ¼
∂ξ2 1 ¼ ∂x02 jm2 j sin θ
where θ is the angle defined in Fig. 1, and sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 ∂x1 ∂x2 ∂x3 jmK j ¼ þ þ ∂ξk ∂ξk ∂ξk
ð7Þ
ð8Þ ð9Þ
The local displacement components u0I in Eq. (6) can be expressed in terms of the nodal displacements referred to the global axis system, i.e., u0I ¼ LIj uj ¼ LIj
M X
N l ulj
ð10Þ
ð14Þ
After some rather complicate reorganization and simplification, the stresses in terms of the global Cartesian equivalents can be expressed as follows:
where the subscript K ranges from 1 to 2. ∂ξK =∂x0J are the derivatives of the intrinsic coordinates ξ with respect to the local coordinates x0J , the detailed expressions of which can be seen in Refs. [3,13], i.e.,
1 ∂xi ∂xi jm1 jjm2 j ∂ξ1 ∂ξ2
ð13Þ
Now, we have obtained the expressions for the stresses in the local coordinate system. Finally, employ the transformation relationships
In order to get the expressions for stresses in terms of known displacements and tractions instead of the strains in Eq. (4), first, we introduce a local Cartesian coordinate system xI' as shown in Fig. 1. The coordinate axis x3' is directed in the normal direction while x1' and x2' are tangential to the surface. The local tangential strains can be expressed in terms of the differentials of the local displacements: ! 0 1 ∂u0I ∂uJ ε0IJ ¼ þ ð5Þ 2 ∂x0J ∂x0I
cos θ ¼
ð11Þ
σ 033 ¼ t 03 ¼ L3j t j
3. Traction-recovery method for 3D
∂ξ2 ¼ 0; ∂x01
227
1 2ν ~ kðxÞðδmn L3m L3n Þ 1ν
ð15Þ
ð16Þ
ð17Þ ð18Þ
From Eq. (15), it can be seen that the evaluation of the thermal stresses at boundary points depends on the nodal displacements, tractions and temperature change of the boundary element. Usually, displacements and tractions are evaluated using the displacement integral equation, while the temperature changes are given by thermal analysis with BEM or other means. Based on the equilibrium equation σ ij;j ¼ 0 without considering body forces and introducing the weight function U ij ðxp ; xÞ, the following weak form of the equilibrium equation can be written out as Z U ij ðxp ; xÞσ jk;k ðxÞdΩðxÞ ¼ 0 ð19Þ Ω
l¼1
where LIj are the direction cosines of the local coordinate system with regard to the global coordinate system [3], M is the total number of boundary nodes over a boundary element, and N l are
Substituting Eq. (2) into Eq. (19) and using Gauss's divergence theorem, the displacement boundary-domain integral equations for non-homogeneous thermoelasticity can be derived as follows
228
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
[16,18]:
Z
coefficients Amnjl , Bmnj and C mn can be written as
Z
T ij ðxp ; xÞu~ j ðxÞdΓðxÞ Z Z ~ þ V ij ðxp ; xÞu~ j ðxÞdΩðxÞ þ U ij;j ðxp ; xÞθðxÞdΩðxÞ
u~ i ðxp Þ ¼
Γ
U ij ðxp ; xÞt j ðxÞdΓðxÞ
Amnjl ¼
Γ
Ω
Ω
ð20Þ
where, xp is the source point. The normalized displacement u~ i , temperature θ~ and fundamental solutions U ij ; T ij ; V ij and U ij;j can be found in Refs. [3,17]. The domain integral appearing in Eq. (20) is transformed into the equivalent boundary integral by employing the radial integration method [19] such that the Eq. (20) is converted to the pure boundary integral equation Z Z u~ i ðxp Þ ¼ U ij ðxp ; xÞt j ðxÞdΓðxÞ T ij ðxp ; xÞu~ j ðxÞdΓðxÞ Γ
Γ
Z
1 ∂r u p F ðx ; yÞdΓðyÞ þ þ α p Γ r ðx ; yÞ ∂n
Z
1 ∂r θ p F ðx ; yÞdΓðyÞ α p Γ r ðx ; yÞ ∂n ð21Þ
u
p
θ
p
1 2ν L1m L1n L2j Bmnj ¼ ðL1m L2n þ L2m L1n ÞL1j þ δmn 1 ν α βν ~ kðxÞL1m L1n 1ν
C mn ¼
F ðx ; yÞ ¼ F ðx ; yÞ ¼
R rðxp ;yÞ 0
V ij ðx ; xÞu~ j ðxÞr dr
0
α ~ U ij;j ðxp ; xÞθðxÞr dr
R rðxp ;yÞ
p
∂ξ 1 ¼ ∂x01 J ðξÞ
In Eqs. (21) and (22), rðxp ; yÞis the distance between the source point xp and the boundary point y, and ∂r=∂n ¼ r ;i ni . After the boundary Γ is discretized into boundary elements, the displacements and tractions can be obtained by solving the algebraic equations of the system from the integral equation Eq. (21). Then, substituting the boundary displacement and traction into Eq. (15), the thermal stress over the boundary Γ can be also computed.
ð24Þ ð25Þ
ð26Þ
where J ðξÞ is the Jacobian of the transformation from the global system ðx1 ; x2 Þ to the local system ðξÞ. In plane strain, the above coefficients need to be supplemented by the additional terms 2μðxÞν ∂ξ ∂Nl L 1 ν 1j ∂x'1 ∂ξ
ð27Þ
B33j ¼
ν L 1 ν 2j
ð28Þ
C 33 ¼
1 2ν ~ kðxÞ 1ν
ð29Þ
α
ð22Þ
ð23Þ
where
A33jl ¼
where
2μðxÞ ∂ξ ∂N l L1m L1n L1j 0 1ν ∂xJ ∂ξ
5. Numerical examples To verify the correctness and the accuracy of the present method for solving boundary thermal-stresses in FGMs, two numerical examples for plane stress and 3D problems of FGMs are presented in the following.
4. Traction-recovery method for 2D
5.1. 2D rectangular FGMs plate
For 2D problems, the analogous problem arises as well and can be solved in a similar manner. In this case, the local coordinate ξ is tangential to the boundary element, while x'1 and x'2 are the local Cartesian axes through the point of interest, as shown in Fig. 2. In the two-dimensional case, the thermal stresses at boundary points can also be expressed by Eq. (15), just only the subscripts ðm; n; j; lÞ in Eqs. (15)–(18) range from 1 to 2. Therefore, the
The first example is concerned with a rectangular functionally graded plate subjected to an inhomogeneous temperature distribution. As is shown in Fig. 3, the left and right edges of the plate are given the temperate 100 K and 30 K, respectively, while the top and bottom edges are q¼ 0. The bottom edge is fixed and the other three edges are traction-free. This problem can be reduced to a plane stress problem and the geometry of the plate is shown in Fig. 3. The shear modulus and thermal expansion coefficient vary in the y direction, which are described by μ ¼ ð50 þ 300yÞGPa and k ¼ ð1:0 þ 15yÞn10 5 K 1 , respectively. The thermal conductivity and Poisson ratio are assumed to be constant with λ ¼ 0:015 W=ðmKÞ and ν ¼ 0:3. For the BEM model, both top and bottom edges are discretized into 40 equally spaced linear boundary elements while the left and right edges into 20 elements as shown in Fig. 4a. Totally, there are 120 boundary elements, 120 boundary points and 741 interior points. For the purpose of verification, the problem is also computed using the FEM software ANSYS. The FEM model is shown in Fig. 4b.
Fig. 1. Local orthogonal set of axes over a boundary element.
Fig. 2. Local orthogonal set of axes over a line boundary element.
Fig. 3. Rectangular functionally graded plate.
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
229
Fig. 6. Normal stress σ yy over bottom edge.
Fig. 4. Computational model of rectangular plate: (a) BEM, and (b) FEM.
Fig. 7. Shear stress σ xy over bottom edge.
Fig. 5. Normal stress σ xx over bottom edge.
The computed normal stresses σ xx and σ yy at boundary points along the line of y¼ 0 are given in Figs. 5 and 6 while the shear stress σ xy is shown in Fig. 7. From Figs. 5–7, we can find that the results computed by traction-recovery method are in good agreement with that obtained by the software ASNSY. This indicates the present method for evaluating boundary thermal stress is correct and can efficiently solve the plane stress problem. 5.2. 3D thick-walled cylinder with FGMs The second example to be considered is a thick-walled cylinder subjected to an internal pressure P ¼ 10 Pa. The inner and outer radiuses are 10 cm and 20 cm as depicted in Fig. 8. Both the shear modulus and thermal expansion coefficient vary along the radial direction, which are described by μ ¼ ð20 þ 300RÞGPa and k ¼ ð 0:5 þ 15RÞn10 5 K 1 . The Poisson ratio is selected as υ ¼ 0:3. The uniform temperature change of the whole cylinder is θ ¼ 100 K.
Fig. 8. Section of thick-walled cylinder.
The height of the cylinder is assumed z¼ 20 cm. Two end sections (z ¼0.0 cm and z¼ 20.0 cm) are constrained in the zdirection. Due to the symmetry, only a quarter of the cylinder is selected for analyzing as shown in Fig. 9a. The internal and external surfaces are discretized into 12 equally-spaced quadratic elements in hoop direction and 6 equally-spaced elements in the z-direction, and the other surfaces along the radial direction are discretized into 4 equally-spaced quadratic boundary elements. Totally, there are 866 boundary nodes and 288 boundary elements in BEM model.
230
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
Fig. 11. Normal stress σ yy over internal hoop line.
Fig. 12. Shear stress σ xy over internal hoop line.
Fig. 9. Computational model of 3D thick-walled cylinder: (a) BEM, and (b) FEM.
Along the internal hoop line of the z ¼10 cm plane, boundary stresses σ xx , σ yy and σ xy are given in Figs. 10–12. From Figs. 10–12 we can see that both the normal and shear stresses computed by the approach proposed in this paper coincide well with those of ANSYS. This means that the method presented in this paper can precisely evaluate the boundary thermal-stress in 3D problems with FGMs.
6. Concluding discussions A robust approach to calculate the boundary stress in thermal stress analysis of structures consisting of FGMs is proposed based on the computed displacements and tractions. Two numerical examples for plane stress and 3D problems are given. The computed results show that the method presented in the paper can accurately evaluate the boundary thermal-stress for functionally graded material structures.
Acknowledgments The authors gratefully acknowledge the National Natural Science Foundation of China (11172055, 51206014). Fig. 10. Normal stress σ xx over internal hoop line.
Due to the condition described above, the deformation of any intersection plane paralleling to the x–y plane is the same. Therefore, a 2D plane strain ANSYS results can be used to verify the current results in the x–y plane, whose model is shown in Fig. 9b.
References [1] Banerjee PK. The boundary element methods in engineering. London: McGraw-Hill Book Co.; 1994. [2] Beer G. Programming the boundary element method. An introduction for engineers. Chichester: John Wiley & Sons, Ltd.; 2001.
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
[3] Gao XW, Davies TG. Boundary element programming in mechanics. Cambridge: Cambridge University Press; 2002. [4] Beer G, Smith I, Duenser C. The boundary element method with programming. Wien: Springer-Verlag; 2008. [5] Rizzo FJ. An integral equation approach to boundary value problems of classical elastostatics. Q Appl Math 1967;25:83–95. [6] Guiggiani M. Hyper singular formulation for boundary stresses evaluation. Eng Anal Bound Elem 1994;14:169–79. [7] Sladek V, Sladek J. Singular integrals and boundary elements. Comput Methods Appl Mech Eng 1998;157:251–66. [8] Matsumoto T, Tanaka M, Hirata H. Boundary stress calculation using regularized boundary integral equation for displacement gradients. Int J Numer Methods Eng 1993;36:783–97. [9] Mantic V. On computing boundary limiting values of boundary integrals with strongly singular and hypersingular kernels in 3d BEM for elastostatics. Eng Anal Bound Elem 1992;13:115–34. [10] Tanaka M, Sladek V, Sladek J. Regularization techniques applied to boundary element methods. Appl Mech Rev 1994;47:457–99. [11] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary elements. Comput Methods Appl Mech Eng 2010;199:2856–64.
231
[12] Wang J, Tsay TK. Analytical evaluation and application of the singularities in boundary element method. Eng Anal Bound Elem 2005;29:241–56. [13] Gao XW, Feng WZ, Yang K, Cui M. Projection plane method for evaluation of arbitrary high order singular boundary integrals. Eng Anal Bound Elem 2015;50:265–74. [14] Becker AA. The boundary element method in engineering. London: McGrawHill Book Co.; 1992. [15] Banerjee PK, Raveendra ST. Advanced boundary element analysis of two-and three-dimensional problems of elasto-plasticity. Int J Numer Methods Eng 1986;23:985–1002. [16] Brebbia CA, Dominguez J. Boundary elements: an introductory course. Southampton and Boston: Computational Mechanics Publicatons and Mc Graw-Hill Book Co.; 1992. [17] Gao XW, Yang K. Thermal stress analysis of functionally graded material structures using boundary element method. Chin J Theor Appl Mech 2011;43:136–43. [18] Lachat JC. A further development of the boundary integral technique for elastostatics. UK: University of Southampton; 1975 Ph.D. thesis. [19] Gao XW. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Bound Elem 2002;26:905–16.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A traction-recovery method for evaluating boundary stresses on thermal elasticity problems of FGMs Jian Liu, Hai-Feng Peng, Xiao-Wei Gao, Miao Cui n State Key Laboratory of Structural Analysis for Industrial Equipment, School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China
art ic l e i nf o
a b s t r a c t
Article history: Received 13 November 2014 Received in revised form 9 June 2015 Accepted 22 July 2015 Available online 27 August 2015
A new approach is presented to calculate boundary stresses in thermal stress analysis of structures consisting of functionally grades materials (FGMs) based on the traction-recovery method. In this approach, the in-plane strains are calculated first using the computed nodal displacements by simply differentiating shape functions at the point of interest, and then the boundary stresses are recovered by Hooke's law together with the known tractions on the boundary. This approach has the advantage of without need to evaluate strongly singular boundary integrals. With the comparison to the FEM software ANSYS, two numerical examples for plane stress and 3D problems are presented to verify the correctness of the proposed method in evaluating boundary thermal-stresses of FGMs. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Traction-recovery method Boundary stresses Thermo-elasticity Functionally graded material Boundary element method
1. Introduction In recent decades, the boundary element method (BEM) has attracted growing interest and appears to be an effective alternative to the finite element method (FEM) for certain engineering problems [1–4], thanks to the reduction of dimensionality, the superior in solving infinite or semi-infinite field problems, and the highly accurate solution on the boundary, etc. However, one of the issues that have still received a great deal of attention is how to accurately determine the boundary stresses for some complicated thermal-stress problems such as non-homogenous and nonisotropic problems [5–8]. It has been found that the determinations of the stresses inside the domain and on the boundary are quite different [2–5]. The stresses inside the domain can be directly evaluated using the stress boundary integral equations. However, the stresses on the boundary are not easy to be evaluated in the same way, due to the hyper-singularities existing in the stress boundary integrals. Therefore, for many years the computation of boundary stresses is an important research issue in the BEM community. Two main methods are used to evaluate the boundary stresses or related high order singular boundary integrals. The first one is the direct method [6–13], which follows the same procedure as computing the internal stresses. Although this method is
n
Corresponding author. E-mail address: [email protected] (M. Cui).
http://dx.doi.org/10.1016/j.enganabound.2015.07.016 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
considered to be more accurate, it does require a substantially much more computational time due to the cost of numerically evaluating all boundary integrals [14]. What's worse, evaluating the high order singular boundary integrals is really an enormous challenge for the BEM researchers. The second method, the most popular one, to evaluate boundary stresses in BEM is the traction-recovery method (TRM) [3,5,15,16], which computes the boundary stresses from the computed displacements and tractions by differentiating shape functions of the boundary element under consideration. Due to the differentiation of shape functions with respect to intrinsic coordinates, the accuracy using TRM is lower than the direct method. On the other hand, due to the use of the precisely computed traction on the surface to recover normal direction related stresses, the accuracy of TRM is higher than the finite element method (FEM). For this reason, this method has been widely used. For the traction-recovery method, Gao and Daivies [3] gives a complete set of formulations and Fortran codes for isotropic materials without consideration of the thermal effect. However, up to date, no detailed formulations are published for the thermal-stress analysis of FGMs. This paper tries to close this gap for the first time.
2. Strain–stress relationship in thermo-elasticity for FGMs Considering the volumetric expansion caused by temperature change over a functionally graded material, the relationship between the total strain εij and the stress σ ij in thermo-elasticity
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
can be expressed as [17,18] 1 ν σ ij δij σ kk þδij kðxÞθ εij ¼ 2μðxÞ 1þν
ð1Þ
where the shear modulus μ and thermal expansion coefficient k are the functions of the coordinates x, Poisson's ratio ν is constant, and θ is the temperature change. Eq. (1) can be converted to νδij εkk ~ ð2Þ δij kðxÞθ σ ij ¼ 2μðxÞ εij þ 1 2ν where 2μðxÞð1 þ νÞkðxÞ ~ kðxÞ ¼ 1 2ν
ð3Þ
By eliminating the local normal strain ε33 in Eq. (2), we can obtain the relationship between the stresses and strains in the tangential plane as follows: 2μðxÞ ν ~ ~ ðσ 33 þ kðxÞθÞ ðε11 þ νε22 Þ þ kðxÞθ 1ν 1ν 2μðxÞ ν ~ ~ ðσ 33 þ kðxÞθÞ σ 22 ¼ ðε22 þ νε11 Þ þ kðxÞθ 1ν 1ν σ 12 ¼ 2μðxÞε12 σ 11 ¼
the shape functions of the boundary element [3]. For 2-noded linear boundary element, N l can be expressed as: 1 N1 ¼ ð1 ξÞ; 2
1 N2 ¼ ð1 þ ξÞ 2
where ξ is the local coordinate ½ 1; þ 1. Using the above Eqs. (5)–(10), the local tangential strains ε0IJ can be obtained. Substituting the results into Eq. (4) yields ∂uj 2μðxÞ ∂ξK ∂ξK ν ~ ~ L t þ kðxÞθ L þ ν L þ kðxÞθ σ 011 ¼ 1j 2j 1 ν ∂x01 ∂x01 ∂ξK 1 ν 3j j ∂uj 2μðxÞ ∂ξK ∂ξ ν ~ ~ L3j t j þ kðxÞθ σ 022 ¼ L2j þ ν K0 L1j þ kðxÞθ 0 1 ν ∂x2 ∂x1 ∂ξK 1 ν ∂uj ∂ξK ∂ξ σ 012 ¼ μðxÞ L2j þ K0 L1j ð12Þ 0 ∂x1 ∂x2 ∂ξK According to the equilibrium equation, the other local stress components related to the normal direction can be recovered from the tractions. σ 023 ¼ t 02 ¼ L2j t j
ð4Þ
σ 013 ¼ t 01 ¼ L1j t j
where the upper-case subscripts I, J range from 1 to 2 and indicate the values refer to the local axis system. The local derivatives of displacement can be obtained by chain rule ∂u0I ∂u0I ∂ξK ¼ ∂x0J ∂ξK ∂x0J
ð6Þ
σ mn ¼ Lrm Lsn σ 0rs
σ mn ¼ Amnjl ulj þBmnj t j C mn θ where the coefficients Amnjl , Bmnj and C mn are 1 ∂ξK ∂ξ L1m L1n Amnjl ¼ 2μðxÞ L1j þ ν K0 L2j 0 1ν ∂x1 ∂x2
∂ξK ∂ξK þ L2m L2n L2j þ ν 0 L1j : ∂x02 ∂x1 1 ∂ξK ∂ξK ∂N l L þ L þ ðL1m L2n þ L2m L1n Þ 2j 1j 2 ∂x01 ∂x02 ∂ξK Bmnj ¼ ðL3m L1n þ L1m L3n ÞL1j þ ðL2m L3n þ L3m L2n ÞL2j ν 1 2ν þ δmn þ L3m L3n L3j 1ν 1ν
∂ξ1 1 ∂ξ1 cos θ ; ; ¼ ¼ ∂x01 jm1 j ∂x02 jm1 j sin θ
C mn ¼
∂ξ2 1 ¼ ∂x02 jm2 j sin θ
where θ is the angle defined in Fig. 1, and sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 ∂x1 ∂x2 ∂x3 jmK j ¼ þ þ ∂ξk ∂ξk ∂ξk
ð7Þ
ð8Þ ð9Þ
The local displacement components u0I in Eq. (6) can be expressed in terms of the nodal displacements referred to the global axis system, i.e., u0I ¼ LIj uj ¼ LIj
M X
N l ulj
ð10Þ
ð14Þ
After some rather complicate reorganization and simplification, the stresses in terms of the global Cartesian equivalents can be expressed as follows:
where the subscript K ranges from 1 to 2. ∂ξK =∂x0J are the derivatives of the intrinsic coordinates ξ with respect to the local coordinates x0J , the detailed expressions of which can be seen in Refs. [3,13], i.e.,
1 ∂xi ∂xi jm1 jjm2 j ∂ξ1 ∂ξ2
ð13Þ
Now, we have obtained the expressions for the stresses in the local coordinate system. Finally, employ the transformation relationships
In order to get the expressions for stresses in terms of known displacements and tractions instead of the strains in Eq. (4), first, we introduce a local Cartesian coordinate system xI' as shown in Fig. 1. The coordinate axis x3' is directed in the normal direction while x1' and x2' are tangential to the surface. The local tangential strains can be expressed in terms of the differentials of the local displacements: ! 0 1 ∂u0I ∂uJ ε0IJ ¼ þ ð5Þ 2 ∂x0J ∂x0I
cos θ ¼
ð11Þ
σ 033 ¼ t 03 ¼ L3j t j
3. Traction-recovery method for 3D
∂ξ2 ¼ 0; ∂x01
227
1 2ν ~ kðxÞðδmn L3m L3n Þ 1ν
ð15Þ
ð16Þ
ð17Þ ð18Þ
From Eq. (15), it can be seen that the evaluation of the thermal stresses at boundary points depends on the nodal displacements, tractions and temperature change of the boundary element. Usually, displacements and tractions are evaluated using the displacement integral equation, while the temperature changes are given by thermal analysis with BEM or other means. Based on the equilibrium equation σ ij;j ¼ 0 without considering body forces and introducing the weight function U ij ðxp ; xÞ, the following weak form of the equilibrium equation can be written out as Z U ij ðxp ; xÞσ jk;k ðxÞdΩðxÞ ¼ 0 ð19Þ Ω
l¼1
where LIj are the direction cosines of the local coordinate system with regard to the global coordinate system [3], M is the total number of boundary nodes over a boundary element, and N l are
Substituting Eq. (2) into Eq. (19) and using Gauss's divergence theorem, the displacement boundary-domain integral equations for non-homogeneous thermoelasticity can be derived as follows
228
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
[16,18]:
Z
coefficients Amnjl , Bmnj and C mn can be written as
Z
T ij ðxp ; xÞu~ j ðxÞdΓðxÞ Z Z ~ þ V ij ðxp ; xÞu~ j ðxÞdΩðxÞ þ U ij;j ðxp ; xÞθðxÞdΩðxÞ
u~ i ðxp Þ ¼
Γ
U ij ðxp ; xÞt j ðxÞdΓðxÞ
Amnjl ¼
Γ
Ω
Ω
ð20Þ
where, xp is the source point. The normalized displacement u~ i , temperature θ~ and fundamental solutions U ij ; T ij ; V ij and U ij;j can be found in Refs. [3,17]. The domain integral appearing in Eq. (20) is transformed into the equivalent boundary integral by employing the radial integration method [19] such that the Eq. (20) is converted to the pure boundary integral equation Z Z u~ i ðxp Þ ¼ U ij ðxp ; xÞt j ðxÞdΓðxÞ T ij ðxp ; xÞu~ j ðxÞdΓðxÞ Γ
Γ
Z
1 ∂r u p F ðx ; yÞdΓðyÞ þ þ α p Γ r ðx ; yÞ ∂n
Z
1 ∂r θ p F ðx ; yÞdΓðyÞ α p Γ r ðx ; yÞ ∂n ð21Þ
u
p
θ
p
1 2ν L1m L1n L2j Bmnj ¼ ðL1m L2n þ L2m L1n ÞL1j þ δmn 1 ν α βν ~ kðxÞL1m L1n 1ν
C mn ¼
F ðx ; yÞ ¼ F ðx ; yÞ ¼
R rðxp ;yÞ 0
V ij ðx ; xÞu~ j ðxÞr dr
0
α ~ U ij;j ðxp ; xÞθðxÞr dr
R rðxp ;yÞ
p
∂ξ 1 ¼ ∂x01 J ðξÞ
In Eqs. (21) and (22), rðxp ; yÞis the distance between the source point xp and the boundary point y, and ∂r=∂n ¼ r ;i ni . After the boundary Γ is discretized into boundary elements, the displacements and tractions can be obtained by solving the algebraic equations of the system from the integral equation Eq. (21). Then, substituting the boundary displacement and traction into Eq. (15), the thermal stress over the boundary Γ can be also computed.
ð24Þ ð25Þ
ð26Þ
where J ðξÞ is the Jacobian of the transformation from the global system ðx1 ; x2 Þ to the local system ðξÞ. In plane strain, the above coefficients need to be supplemented by the additional terms 2μðxÞν ∂ξ ∂Nl L 1 ν 1j ∂x'1 ∂ξ
ð27Þ
B33j ¼
ν L 1 ν 2j
ð28Þ
C 33 ¼
1 2ν ~ kðxÞ 1ν
ð29Þ
α
ð22Þ
ð23Þ
where
A33jl ¼
where
2μðxÞ ∂ξ ∂N l L1m L1n L1j 0 1ν ∂xJ ∂ξ
5. Numerical examples To verify the correctness and the accuracy of the present method for solving boundary thermal-stresses in FGMs, two numerical examples for plane stress and 3D problems of FGMs are presented in the following.
4. Traction-recovery method for 2D
5.1. 2D rectangular FGMs plate
For 2D problems, the analogous problem arises as well and can be solved in a similar manner. In this case, the local coordinate ξ is tangential to the boundary element, while x'1 and x'2 are the local Cartesian axes through the point of interest, as shown in Fig. 2. In the two-dimensional case, the thermal stresses at boundary points can also be expressed by Eq. (15), just only the subscripts ðm; n; j; lÞ in Eqs. (15)–(18) range from 1 to 2. Therefore, the
The first example is concerned with a rectangular functionally graded plate subjected to an inhomogeneous temperature distribution. As is shown in Fig. 3, the left and right edges of the plate are given the temperate 100 K and 30 K, respectively, while the top and bottom edges are q¼ 0. The bottom edge is fixed and the other three edges are traction-free. This problem can be reduced to a plane stress problem and the geometry of the plate is shown in Fig. 3. The shear modulus and thermal expansion coefficient vary in the y direction, which are described by μ ¼ ð50 þ 300yÞGPa and k ¼ ð1:0 þ 15yÞn10 5 K 1 , respectively. The thermal conductivity and Poisson ratio are assumed to be constant with λ ¼ 0:015 W=ðmKÞ and ν ¼ 0:3. For the BEM model, both top and bottom edges are discretized into 40 equally spaced linear boundary elements while the left and right edges into 20 elements as shown in Fig. 4a. Totally, there are 120 boundary elements, 120 boundary points and 741 interior points. For the purpose of verification, the problem is also computed using the FEM software ANSYS. The FEM model is shown in Fig. 4b.
Fig. 1. Local orthogonal set of axes over a boundary element.
Fig. 2. Local orthogonal set of axes over a line boundary element.
Fig. 3. Rectangular functionally graded plate.
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
229
Fig. 6. Normal stress σ yy over bottom edge.
Fig. 4. Computational model of rectangular plate: (a) BEM, and (b) FEM.
Fig. 7. Shear stress σ xy over bottom edge.
Fig. 5. Normal stress σ xx over bottom edge.
The computed normal stresses σ xx and σ yy at boundary points along the line of y¼ 0 are given in Figs. 5 and 6 while the shear stress σ xy is shown in Fig. 7. From Figs. 5–7, we can find that the results computed by traction-recovery method are in good agreement with that obtained by the software ASNSY. This indicates the present method for evaluating boundary thermal stress is correct and can efficiently solve the plane stress problem. 5.2. 3D thick-walled cylinder with FGMs The second example to be considered is a thick-walled cylinder subjected to an internal pressure P ¼ 10 Pa. The inner and outer radiuses are 10 cm and 20 cm as depicted in Fig. 8. Both the shear modulus and thermal expansion coefficient vary along the radial direction, which are described by μ ¼ ð20 þ 300RÞGPa and k ¼ ð 0:5 þ 15RÞn10 5 K 1 . The Poisson ratio is selected as υ ¼ 0:3. The uniform temperature change of the whole cylinder is θ ¼ 100 K.
Fig. 8. Section of thick-walled cylinder.
The height of the cylinder is assumed z¼ 20 cm. Two end sections (z ¼0.0 cm and z¼ 20.0 cm) are constrained in the zdirection. Due to the symmetry, only a quarter of the cylinder is selected for analyzing as shown in Fig. 9a. The internal and external surfaces are discretized into 12 equally-spaced quadratic elements in hoop direction and 6 equally-spaced elements in the z-direction, and the other surfaces along the radial direction are discretized into 4 equally-spaced quadratic boundary elements. Totally, there are 866 boundary nodes and 288 boundary elements in BEM model.
230
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
Fig. 11. Normal stress σ yy over internal hoop line.
Fig. 12. Shear stress σ xy over internal hoop line.
Fig. 9. Computational model of 3D thick-walled cylinder: (a) BEM, and (b) FEM.
Along the internal hoop line of the z ¼10 cm plane, boundary stresses σ xx , σ yy and σ xy are given in Figs. 10–12. From Figs. 10–12 we can see that both the normal and shear stresses computed by the approach proposed in this paper coincide well with those of ANSYS. This means that the method presented in this paper can precisely evaluate the boundary thermal-stress in 3D problems with FGMs.
6. Concluding discussions A robust approach to calculate the boundary stress in thermal stress analysis of structures consisting of FGMs is proposed based on the computed displacements and tractions. Two numerical examples for plane stress and 3D problems are given. The computed results show that the method presented in the paper can accurately evaluate the boundary thermal-stress for functionally graded material structures.
Acknowledgments The authors gratefully acknowledge the National Natural Science Foundation of China (11172055, 51206014). Fig. 10. Normal stress σ xx over internal hoop line.
Due to the condition described above, the deformation of any intersection plane paralleling to the x–y plane is the same. Therefore, a 2D plane strain ANSYS results can be used to verify the current results in the x–y plane, whose model is shown in Fig. 9b.
References [1] Banerjee PK. The boundary element methods in engineering. London: McGraw-Hill Book Co.; 1994. [2] Beer G. Programming the boundary element method. An introduction for engineers. Chichester: John Wiley & Sons, Ltd.; 2001.
J. Liu et al. / Engineering Analysis with Boundary Elements 61 (2015) 226–231
[3] Gao XW, Davies TG. Boundary element programming in mechanics. Cambridge: Cambridge University Press; 2002. [4] Beer G, Smith I, Duenser C. The boundary element method with programming. Wien: Springer-Verlag; 2008. [5] Rizzo FJ. An integral equation approach to boundary value problems of classical elastostatics. Q Appl Math 1967;25:83–95. [6] Guiggiani M. Hyper singular formulation for boundary stresses evaluation. Eng Anal Bound Elem 1994;14:169–79. [7] Sladek V, Sladek J. Singular integrals and boundary elements. Comput Methods Appl Mech Eng 1998;157:251–66. [8] Matsumoto T, Tanaka M, Hirata H. Boundary stress calculation using regularized boundary integral equation for displacement gradients. Int J Numer Methods Eng 1993;36:783–97. [9] Mantic V. On computing boundary limiting values of boundary integrals with strongly singular and hypersingular kernels in 3d BEM for elastostatics. Eng Anal Bound Elem 1992;13:115–34. [10] Tanaka M, Sladek V, Sladek J. Regularization techniques applied to boundary element methods. Appl Mech Rev 1994;47:457–99. [11] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary elements. Comput Methods Appl Mech Eng 2010;199:2856–64.
231
[12] Wang J, Tsay TK. Analytical evaluation and application of the singularities in boundary element method. Eng Anal Bound Elem 2005;29:241–56. [13] Gao XW, Feng WZ, Yang K, Cui M. Projection plane method for evaluation of arbitrary high order singular boundary integrals. Eng Anal Bound Elem 2015;50:265–74. [14] Becker AA. The boundary element method in engineering. London: McGrawHill Book Co.; 1992. [15] Banerjee PK, Raveendra ST. Advanced boundary element analysis of two-and three-dimensional problems of elasto-plasticity. Int J Numer Methods Eng 1986;23:985–1002. [16] Brebbia CA, Dominguez J. Boundary elements: an introductory course. Southampton and Boston: Computational Mechanics Publicatons and Mc Graw-Hill Book Co.; 1992. [17] Gao XW, Yang K. Thermal stress analysis of functionally graded material structures using boundary element method. Chin J Theor Appl Mech 2011;43:136–43. [18] Lachat JC. A further development of the boundary integral technique for elastostatics. UK: University of Southampton; 1975 Ph.D. thesis. [19] Gao XW. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng Anal Bound Elem 2002;26:905–16.