Heat Mass Transfer (2012) 48:283–290 DOI 10.1007/s00231-011-0881-x
ORIGINAL
An iterative method to solve the heat transfer problem under the non-linear boundary conditions Zhenggang Zhu • Michael Kaliske
Received: 19 August 2009 / Accepted: 31 July 2011 / Published online: 18 August 2011 Ó Springer-Verlag 2011
Abstract The aim of the paper is to determine the approximation of the tangential matrix for solving the nonlinear heat transfer problem. Numerical model of the strongly non-linear heat transfer problem based on the theory of the finite element method is presented. The tangential matrix of the Newton method is formulated. A method to solve the heat transfer with the non-linear boundary conditions, based on the secant slope of a reference function, is developed. The contraction mapping principle is introduced to verify the convergence of this method. The application of the method is shown by two examples. Numerical results of these examples are comparable to the ones solved with the Newton method and the commercial software COMSOL for the heat transfer problem under the radiative boundary conditions.
1 Introduction Iterative methods are always used to solve strongly nonlinear problems. Various methods, based on the approximation of the tangential matrix, have been developed. One of the most common used methods is the Newton method [29], which has the most rapidly convergent rate only if the initial values used in iterations are within the zone of attraction and, thus, divergence does not occur. However, the tangential matrix has to be updated for each iteration, which makes it necessary to use a large memory space Z. Zhu (&) M. Kaliske Institute for Structural Analysis, Technische Universita¨t Dresden, 01062 Dresden, Germany e-mail:
[email protected] M. Kaliske e-mail:
[email protected]
during the iterative process and computer may even be overburdened during the exchange of large amounts of data between the memory and the hard disks. Furthermore, the tangential matrix remains tough for a strongly non-linear temperature field [16]. In order to find the tangent, the iterative methods, such as a return mapping algorithm [22] for an elastoplasticity in the field of structural mechanics, are always needed. In this case, a constant approximation of the tangential matrix, known as the modified NewtonRaphson method, may be more effective in an economical usage of the computer memory [25, 29]. A relaxation factor, which avoids the direct evaluation of tangential matrix, is used to solve the heat conduction under the nonlinear boundary conditions [8, 12], although the number of iteration cycles and the convergence are affected greatly by the used relaxation factor. The Picard method (also known as the successive substitution method) is introduced by [23] to solve the non-linear heat conduction. These methods, which are more computationally economical than the Newton method, play important roles in solving non-linear problems. However, they may diverge or even fail to find the unique solutions. Analytical or semi-analytical methods are developed for both linear and non-linear heat transfer problems. The linear heat conduction has been widely investigated in [5, 19, 20, 24], where numerical methods, such as the explicit method, the implicit method, the finite element method and the finite difference method, are introduced. A onedimensional heat conduction [5, 18] and a two-dimensional one [17] under radiation boundary conditions is solved with these methods. Analytical methods are introduced for the simply non-linear steady heat conduction with linear boundary conditions [1]. A three-layer scheme with temperature evaluated at central timepoint is presented in [4], although it involves the repeated computation of corrected
123
284
value of the radiation coefficient. The reduction of the nonlinear problem to a linear one is extremely beneficial [6]. A simplified analytical solution method is developed for onedimensional steady thermal field [3], although it is only applicable to those non-linear problems with thermal conductivity linearly varying with temperature. The separation of variables method is presented in [9] for transient heat conduction with the thermal conductivity and heat capacity as the power functions of location. Two set of equations in linear form and exponential one are used for thermal conductivity in order to change the non-linear problem to be linear [6]. However, a general application of these methods to non-linear problems seems impossible. The non-linear heat transfer problem may be solved by a linearization process. Conventionally, the heat transfer coefficient is assumed to be constant within a time step, which, nonetheless, might result in ill accuracy or numerical instability [12]. An analytical method is outlined in [11] by dividing the whole region into a large number of small sub-regions, where the thermal properties can be assumed constant. In this way, the exact solutions within each sub-region are obtained. Nevertheless, the method is applicable only for some special cases, especially when the temperature variation is not so much [11]. The governing equations may also be solved in an explicit manner [2, 7, 10, 14], where the physical coefficients are evaluated on the basis of the solutions from the last time step. In this case, no iterations are needed even in the presence of strong non-linearity, although it is necessary to control the time increment. An alternative way to avoid iterations is presented in [26, 27], where temperature is expanded in a similar form as Taylor series with the unknown coefficients. Then, the resulting equations are solved with the finite element method. A similar way is used in [28]. By expanding an energy density function as a Taylor polynomial in a spatial domain, the partial differential equation is reduced to a set of first-order ordinary differential equations in time [28]. These methods may be theoretically meaningful. Nonetheless, the numbers of unknowns remains to be determined in each time step since equating powers need to be iteratively computed [27], and even worse, these methods may involve a large number of degrees of freedom even for simple structures. The unconditionally stable, implicit backward difference method is more favourable for time discretization [23] with physical coefficients evaluated in the current time step, which necessitates solving iteratively the resulting nonlinear equations. This paper aims to introduce an alternative way for the approximation of the tangential matrix in solving the nonlinear heat conduction under the radiative conditions. The tangential matrix and a modified tangent method, based on a reference function, are presented in the paper.
123
Heat Mass Transfer (2012) 48:283–290
2 Theory Energy balance of the traditional heat conduction [13] can be described by qc
oT ¼ r ðkrTÞ þ Q ot
ð1Þ
where, q is the density, c the heat capacity, k the heat conductivity, Q the heat source, and T the temperature. Using the method of Galerkin weighted residual, the heat balance equation (Eq. 1) without considering the heat source ðQ ¼ 0Þ can be transformed into Z oT T RðT Þ ¼ ½N qc r ðkrT Þ dX ot X Z oT þ ½rN T krT dX ¼ ½N T qc ð2Þ ot X Z oT ½N T k dC on C
where, RðT Þ is the weighted residual dependent on the temperature T evaluated by the nodal values ðT ¼ ½N fT gÞ; ½N the matrix of shape functions, X the structure, C the boundary, oT on the temperature gradient, n the normal direction of surface, r the divergence operator. The first part of the right side is the energy gained ðM Þ in the structure, the second is the heat conduction ðSÞ and the third is the exchange of heat energy ðH Þ of the structure with the surroundings, Z Z oT M ¼ ½N T qc dX; S ¼ ½rN T krTdX; ot X X Z ð3Þ T oT H ¼ ½N k dC on C
As a strongly non-linear problem, the radiation on the boundary [15] is calculated by oT 4 k ¼ er T 4 T1 on ð4Þ 2 ðT þ T1 ÞðT T1 Þ ¼ er T 2 þ T1 ¼ hqr ðT T1 Þ where e is thermal emissivity, r the Stefan–Boltzmann constant, T1 the temperature of the environment. The first part in the right-hand side of Eq. 4 defines the target function ðerT 4 Þ and the reference function hqr in Fig. 1 as 2 hqr ¼ er T 2 þ T1 ð T þ T1 Þ ð5Þ Energy balance on the radiative boundary defined by Eq. 4 yields
Heat Mass Transfer (2012) 48:283–290
285
q0 ¼ erT 4
ð6Þ
where the heat flux q0 in Fig. 1 is calculated by q0 ¼ k
oT 4 þ erT1 on
ð7Þ
Energy balance for the whole structure is thought to be achieved if the weighted residual (Eq. 2) approaches zero. The tangential matrix K used for the Newton method is given by K¼
(a) The Newton method
oR oM oS oH oM oS ¼ þ þ ¼ þ þ HT oT oT oT oT oT oT
ð8Þ
where HT is defined here as the tangent of heat exchange between the structure and its surroundings. A sequence of the approximations of the unknown temperature T obtained during the iterative process are denoted by T ð0Þ ; T ð1Þ ; T ð2Þ ; and T ð3Þ in Fig. 1. With the time derivative substituted by the backward difference method, we get Z oM oðqcÞ T T ð0Þ qc ¼ ½N T þ ½N dX ð9Þ oT oT Mt Mt X
(b) Modified Newton-Raphson method
oS ¼ oT
Z X
oH ¼ oT
Z
0
2
31 N oT ox B T T ok 6 oT 7C @½rN k½rN þ ½rN 4 N oy 5AdX oT N oT oz
ð10Þ
½N T 4erT 3 ½N dC
ð11Þ
C
Given the trial value T ð0Þ evaluated from the results of the last time point, the temperature can be solved by an iterative process MT ðnÞ ¼ K 1 RðnÞ ;
(c) Modified tangent (hqr) method
T ðnþ1Þ ¼ T ðnÞ þ MT ðnÞ
ð12Þ
where the superscript n is the index of iteration number, RðnÞ the weighted residual of the nth iteration, and MT ðnÞ the increment. Various methods can be chosen according to the selection of K used in the iterative process (Eq. 12), for example, K is the tangential matrix (Eq. 8) for the Newton method (Fig. 1a). A constant value of the tangential matrix K can be used during iteration for the modified Newton– Raphson method, though numerical process may fail to converge for such a problem as shown in Fig. 1b, where the tangential matrix K is assumed to be kept as the one obtained at the initial trial value T ð0Þ : 2.1 The contraction mapping principle
(d) Picard method Fig. 1 Iterative process of solving non-linear boundary conditions
The contraction-mapping principle [21], based on the norm of the Jacobian matrix of the recursive form, is applied to the methods introduced in Fig. 1 for a detailed explanation of the convergence of various methods. The recursive form
123
286
Heat Mass Transfer (2012) 48:283–290
for solving the unknown temperature can be formulated from Eq. 12 as T ðnþ1Þ ¼ G T ðnÞ ¼ T ðnÞ K 1 RðnÞ ð13Þ The starting iteration values T ð0Þ can be taken from the solutions obtained in the previous time step. The superscript (n) is omitted for simplification in the following derivations. The Jacobian matrix of Eq. 13 is evaluated by oG oR ¼ I K 1 oT oT
If the eigenvalues of the matrix K and the tangential matrix oR calculated by oT in Eq. 14 are denoted by l and k, respectively, the convergence condition given by inequation 15, substituted by Eq. 14, can be reformulated as 1 l1 k\1 ð16Þ The contraction mapping principle introduced above is used to determine the convergence of the modified tangent (hqr) method (Fig. 1c) under the radiative boundary conditions, where oH oT in Eq. 8 is approximated by Z ð17Þ HT ½N T hqr ½N dC C
The orders of magnitude of l and k may be approximated by the maximum value of those of the right side in Eq. 8 oM oS oH ½k ¼ max ; ; ð18Þ oT oT oT
oM oS ½l ¼ max ð19Þ ; ; ½HT oT oT where ½ denotes the order of magnitude and the derivatives with respect to the temperature are referred to Eqs. 9, 10, and 11. oS oH oH 1. If ½k ¼ max oM ; ; oT ¼ oT and ½l ¼ oM oS oT oT max oT ; oT ; ½HT ¼ ½HT ; then the order of magnitude of the last term in the left-hand side of inequation 16 may be evaluated by
4erT 3 T3 þ T3 þ T3 þ T3 l k ¼ 3 2 T þ T3 hqr T þ T 1 T 2 þ T1 1 1
ð20Þ
As a result, inequation 16 is closely related to two variables in Eq. 20, namely the environmental
123
l1 k
4 T1 2 T1 3 1þ þ T þ T 4 ¼ 1 þ c þ c2 þ c3
ð14Þ
where I is the unity matrix. A sufficient condition to ensure the recursive form of Eq. 13 converging to the unique solution is that the norm of the Jacobian matrix (Eq. 14) must satisfy the contraction-mapping principle [21], namely
oG
\1 ð15Þ
oT
temperature T1 and the surface temperature T. Equation 20, where the ratio of the environmental temperature T1 to the surface temperature T is defined by a variable c, may be reformulated into a simple form as T1 T
ð21Þ
The variable c is positive since it is really hard to achieve absolute zero in temperature in nature (0 K on the Kelvin scale), ð22Þ
c[0
Inequation 16, substituted by Eq. 21, may be recast into 4 1 \1 ð23Þ 2 3 1þcþc þc which yields c3 þ c2 þ c 1 [ 0
ð24Þ
The left side in inequation 24 may be redefined by f ð cÞ ¼ c3 þ c2 þ c 1
ð25Þ
The first derivative of the function f is given by df ¼ 3c2 þ 2c þ 1 [ 0 dc
ð26Þ
which leads to a conclusion that the function f is monotonically increasing. Since f ð0Þ ¼ 1;
f ðþ1Þ [ 0
ð27Þ
then the function f has the unique zero point, which may be iteratively solved with the Newton method. The iterative process is shown as follows given the initial trial value cð0Þ ¼ 0; 1 df f ðiÞ ; (b) loop cðiþ1Þ ¼ cðiÞ dc ðiþ1Þ (c) end loop if c cðiÞ \108 :
(a)
The zero point of the function f is roughly estimated to be 0.5437. Therefore, inequation 24 may yield c [ 0:5437
–
ð28Þ
If the environmental temperature T1 and the surface temperature T are at the same order of magnitude, c1 Equation 20 may follow
ð29Þ
Heat Mass Transfer (2012) 48:283–290
T3
T3 þ T3 þ T3 þ T3 1 2 T þ T3 þ T1 T 2 þ T 1 1
287
ð30Þ
Both, Eqs. 20 and 30, prove that the last term in the right-hand side of the tangent given by Eq. 8 may be roughly estimated by Eq. 5, namely Z oM oS þ þ ½N T hqr dC K ð31Þ oT oT C
–
The convergence condition given by inequation 16, substituted by Eqs. 20 and 30, is fully satisfied. If the environmental temperature is far greater than the surface temperature, c1
–
ð32Þ
inequation 28 is fully satisfied. If the environmental temperature is far less than the surface temperature, c1
ð33Þ
the condition of inequation 28 may be violated and convergence may not be achieved. As a result, the norm of the Jacobian matrix evaluated by the modified tangent (hqr) method does not violate the contraction mapping principle and the iteration based on the modified tangent (hqr) method converges to the unique solutions of the unknown temperature (Fig. 1c) if inequation 16 holds. Otherwise, this method may fail to predict the solutions. Similarly, the last term in the left-hand side of inequation 16 for the Picard method (also known as the successive substitution method) [23] in Fig. 1d may be given by
4erT 3 l1 k ¼ ¼4 erT 3
which yields
oG
[1
oT
ð34Þ
ð35Þ
Therefore, the tangent approximated by this method falls far away from the real tangent of the Newton method. The iteration based on the Picard method may diverge, which is also clearly verified in Fig. 1d. oS oH oH 2. If ½k ¼ max oM ; ; oT ¼ oT and ½l ¼ oM oS oM oS oT oT max oT ; oT ; ½HT ¼ max oT ; oT ; then oH 1 oT l k ¼ [1 ð36Þ oS max oM oT ; oT In this case, environmental temperature is far less than the surface temperature, which makes the smallest order of magnitude of ½HT : The approximation of
tangential matrix may not satisfy the convergence condition (inequation 16). Therefore, the modified tangent (hqr) method may have high risk divergence of iterative process. oS 3. If ½k ¼ max oM oT ; oT ; then oM oS oH 1 max oT ; oT ; oT oS l k ¼ max oM oT ; oT ; ½HT oS ð37Þ max oM ; oT oT oS ¼ max oM oT ; oT ; ½HT A large temperature difference between surface and surroundings may cause oS (a) ½l ¼ max oM oT ; oT ; ½HT ¼ ½HT Then oS 1 max oM oT ; oT \1 ð38Þ l k ¼ ½HT and 1 l1 k\1: oM oS oS (b) ½l ¼ max oM oT ; oT ; ½HT ¼ max oT ; oT Then oS 1 max oM ; oT oT ¼ 1 l k ð39Þ oS max oM oT ; oT and 1 l1 k\1: Therefore, the selection of hqr in Eq. 5 as the approximation of oH oT does not affect the convergence of iteration process. 3 Applications The conception of the tangential matrix K defined by Eq. (31) is used to solve a cube with the size of 1 m 9 1 m 9 1 m (Fig. 2a). Only one 8-node element is applied in order to determine whether the modified tangent (hqr) method is applicable in solving the non-linear radiative boundary conditions. The upper side of the model is set by a constant temperature ð373 KÞ; and the lower side is subjected to a radiative boundary condition according to Eq. 4. The left, the right, the front and the back surfaces are assumed to be heat-insulated. Physical coefficients used in the model are qc ¼ 10 J=ðm3 KÞ;
k ¼ 0:35 W=ðm KÞ;
e ¼ 0:9
ð40Þ
and Stefan–Boltzmann constant is r ¼ 5:67 108 W=ðm2 K4 Þ
ð41Þ
The initial temperature T0 and the environmental one T1 are given by T0 ¼ 293 K;
T1 ¼ 293 K
ð42Þ
respectively. In this case, a constant time step size (2 s) is used and the orders of magnitude (Eq. 18) may be calculated by
123
288
Heat Mass Transfer (2012) 48:283–290
(a) Model
Fig. 3 Convergence rate
T0 ¼ 373 K;
T1 ¼ 100 K
ð45Þ
In this case,
(b) Results of temperature of node at (1,0,0)
101 ¼ 0 10
Fig. 2 Numerical simulation of the three dimensional sample
h i oM oS oH qc : : : ½k : 4erT 3 oT oT oT Mt 0
¼ 10 : 10
1
: 10
ð43Þ
0
and
100 l1 k 0 ¼ 1 10
ð44Þ
Therefore, the modified tangent (hqr) method predicts the unique solutions to the unknown temperature. Temperature fields are solved when the weighted residual falls below the predefined tolerance, such as 10-3. The results (Fig. 2b), compared with the ones calculated by the Newton method, prove that the iterative method using the modified tangent (Eq. 31) is completely acceptable for the heat transfer subjected to the radiative boundary conditions. The convergence rate for the modified tangent method (Eq. 31) shown in Fig. 3 for 10 substeps, which is almost quadratic. However, the convergence of the method is only of first order (Fig. 3) since the matrix K in Eq. 31 is not really the tangent. If the same example is simulated except for a large difference between surface temperature and environmental temperature. The initial temperature T0 and the environmental one T1 are set to be
123
½4erT 3 l1 k hqr
ð46Þ
and 1 l1 k [ 1: Iterative process with the modified tangent (hqr) method fails to converge to the unique solutions. Another example is taken from [16], where a cooling process of a two dimensional sample is analyzed. All four boundaries are set to be radiative, and, therefore, only one quarter of the model (Fig. 4) is simulated. The initial temperature T0, the environmental temperature T1 ; and the physical coefficients, such as the density q, the heat capacity c, the heat conductivity k, the thermal emissivity e; and Stefan–Boltzmann constant r, are given by
Fig. 4 The two-dimensional sample
Heat Mass Transfer (2012) 48:283–290
289
converges to the unique solutions except when the environmental temperature is far less than surface temperature. This physical limitation is explained in the first example, where a large temperature difference between surface and surroundings is defined. The numerical results of both examples solved with the method developed in the paper are compared with the ones of the Newton method and the finite element software COMSOL, which proves that the method is applicable for most cases in solving the nonlinear heat transfer problem under the radiative boundary conditions. Acknowledgments The authors acknowledge gratefully the financial support of the German Academic Exchange Service (DAAD). Fig. 5 Comparison of temperature of node at (0.075, 0.04)
References T0 ¼ 1553 K;
T1 ¼ 298 K;
c ¼ 670 J=ðkg KÞ; e ¼ 0:8;
3
q ¼ 7; 800 kg=m ;
k ¼ 30 W=ðm KÞ;
r ¼ 5:67 108 W=ðm2 K4 Þ:
ð47Þ ð48Þ ð49Þ
In this case, the time step size is set to be 0.1 s. The orders of magnitude of this example (Eq. 18) may be calculated by h i oM oS oH qc : : : ½k : 4erT 3 oT oT oT Mt ð50Þ 7 1 0 ¼ 10 : 10 : 10 and ½HT ¼ 102 ;
107 l1 k 7 ¼ 1 10
ð51Þ
The tangential matrix is approximated by Eq. 31 in these examples. A large temperature difference between surface and surroundings is defined for this example (Eq. 47), although the convergence of the iteration is achieved since 1 l1 k\1: The numerical results of the model solved with the modified tangent (hqr) method are comparable with those obtained from the commercial software COMSOL (Fig. 5).
4 Conclusions How to solve the heat transfer problem under radiative boundary conditions is presented in this paper. Various methods, such as Picard method, Modified Newton-Raphson method, the Newton method and the modified tangent method, are introduced. The tangential matrix is derived. The orders of magnitude are analyzed and the convergence of the modified tangent method is outlined. The iteration with the tangent approximated by a reference function
1. Arslanturk C (2009) Correlation equations for optimum design of annular fins with temperature dependent thermal conductivity. Heat Mass Transf 45:519–525 2. Blobner J, Bialecki RA, Kuhn G (1999) Transient non-linear heat conduction-radiation problems: a boundary element formulation. Int J Numer Methods Eng 46:1865–1882 3. Bouaziz MN, Hanini S (2007) Efficiency and optimisation of fin with temperature-dependent thermal conductivity: a simplified solution. Heat Mass Transf 44:1–9 4. Budacova J (1977) Method of finite elements for solving some heat-conduction problems. J Eng Phys Thermophys 33:728–733 5. Campo A (1977) Unsteady heat transfer from a circular fin with nonlinear dissipation. Heat Mass Transf 10:203–210 6. Campo A (1982) Estimate of the transient conduction of heat in materials with linear thermal properties based on the solution for constant properties. Heat Mass Transf 17:1–9 7. Comini G, Guidice SD, Lewis RW, Zienkiewicz OC (1974) Finite element solution of non-linear heat conduction problems with special reference to phase change. Int J Numer Methods Eng 8:613–624 8. Gururaja Rao C, Nagabhushana Rao V, Krishna Das C (2008) Simulation studies on multi-mode heat transfer from an open cavity with a flush-mounted discrete heat source. Heat Mass Transf 44:727–737 9. Hosseini SM, Akhlaghi M, Shakeri M (2007) Transient heat conduction in functionally graded thick hollow cylinders by analytical method. Heat Mass Transf 43:669–675 10. Khattabi A, Steinhagen P (1993) Analysis of transient nonlinear heat conduction in wood using finite-difference solutions. Holz als Roh-und Werkstoff 51:272–278 11. Kim S, Huang CH (2006) The approximate solutions to the nonlinear heat conduction problems in a semi-infinite medium. Heat Mass Transf 42:727–738 12. Koizumi M, Utamura M, Kotani K (1985) Three-dimensional transient heat conduction analysis with non-linear boundary conditions by boundary element method. J Nucl Sci Technol 22:972–982 13. Kwon YW, Bang H (1997) The finite element method using matlab. CRC Press, Boca Raton 14. Lees M (1966) A linear three-level difference scheme for quasilinear parabolic equations. Math Comput 20:516–522 15. Lewis RW, Nithiarasu P, Seetharamu K (2004) Fundamentals of the finite element method for heat and fluid flow. Wiley, Chichester
123
290 16. Liu FL, Du RY (2005) Nonlinear fem for calculating temperature field. J Hebei Norm Univ (Nat Sci Edition) 29:21–24 17. Mehta RC, Jayachandran T (1988) Finite element analysis of conductive and radiative heating of a thin skin calorimeter. Heat Mass Transf 22:227–230 18. Onur N, Sivrioglu M (1993) Transient heat conduction with uniform heat generation in a slab subjected to convection and radiation cooling. Heat Mass Transf 28:345–349 19. Peaceman DW, Rachford HH Jr (1955) The numerical solution of parabolic and elliptic differential equations. J Soc Ind Appl Math 3:28–41 20. Qin QH (1994) Unconditionally stable fem for transient linear heat conduction analysis. Commun Numer Methods Eng 10:427–435 21. Rheinboldt WC (1998) Methods for solving systems of nonlinear equations. 2nd edn. SIAM, Philadelphia 22. Simo JC, Taylor RL (1986) A return mapping algorithm for plane stress elastioplasticity. Int J Numer Methods Eng 22:649–670
123
Heat Mass Transfer (2012) 48:283–290 23. Singh IV, Tanaka M, Endo M (2007) Meshless method for nonlinear heat conduction analysis of nano-composites. Heat Mass Transf 43:1097–1106 24. Vujicˇic´ MR (2006) Finite element solution of transient heat conduction using iterative solvers. Eng Comput 23:408–431 25. Wang XC (2003) Finite element method (in Chinese). Tsinghua University Press, Beijing 26. Yang H (1999) A new approach of time stepping for solving transfer problems. Commun Numer Methods Eng 15:325–334 27. Yang H, Gao Q (2003) A precise time stepping scheme to solve hyperbolic and parabolic heat transfer problems with radiative boundary condition. Heat Mass Transf 39:571–577 28. Yu J, Yang Y, Campo A (2010) Approximate solution of the nonlinear heat conduction equation in a semi-infinite domain. Math Probl Eng 2010:1–24 29. Zienkiewicz O, Taylor R (2000) The finite element method— volume 2—solid mechanics. 5th edn. Butterworth-Heinemann, Oxford