Standard staggered and staggered Newton schemes in thermo-hydro-mechanical problems

  • Published on

  • View

  • Download

Embed Size (px)


<ul><li><p>ELSEVIER </p><p>Computer methods in applied </p><p>mechanics and englneerlng </p><p>Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>Standard staggered and staggered Newton schemes in thermo-hydro-mechanical problems </p><p>B.A. Schreflera,*, L. Simoni, E. Turskab Istituto di Scienza e Tecnica delle Costruzioni, University of Padova, Via Marzolo 9, 35131 Padova, Italy </p><p>bPolish Academy of Sciences, IFTR, ul. Swietokrzyska 21, W-049 Warszawa, Poland </p><p>Received 11 May 1995; revised 12 March 1996 </p><p>Abstract </p><p>The aim of the paper is to compare the rates of iteration convergence of two staggered numerical schemes: one using successive substitutions and the other the Newton method. In the case of very weakly coupled problems the Newton method is one order faster, but there are specific cases when it is slower than the standard staggered algorithm. In both cases a convergence condition must be verified. The behaviour of the investigated staggered procedures in case of thermo-hydro-mechanical problems are shown on examples. </p><p>1. Introduction </p><p>Current engineering problems often involve very complex, coupled sets of equations for which numerical integration may lead to very large systems. In such cases numerical efficiency becomes a dominant concern. That is why, despite their conditional stability, the staggered schemes play an important role in the design of numerical models (see [l-4]). Also, to circumvent the problem of conditional stability a number of stabilization techniques has been developed [5,6]. </p><p>The standard stalggered procedure is described and analyzed in [2,3,7-91. It is an iteration method of the Gauss-Seidel type used to solve large, coupled, nonsymmetric sets of algebraic equations, which are a discretization of field equations (e.g. two- or three-phase consolidation or thermo-mechanical problems). Usually, the spatial integration has been performed by the finite element method and the time integration b:y the implicit generalized midpoint rule. </p><p>The advantageous feature of the staggered procedure is that it allows to solve the equations sequentially, so we are able to use available numerical codes for simpler problems. The main concept of the staggered strategy is to solve a first block of equations for the first set of field variables, with the other variables frozen. Then solve the remaining block of equations for the second set of field variables, with the updated first variables fixed. This is obtained by performing an appropriate partitioning of the matrices (operators) on the 1.h.s. and transferring one component to the r.h.s. of the equation. In the linear case the first set of equations can be substituted into the second, thus giving only one set of equations to be iterated (for examples see [l] and Refs. in [2]). The conditions for convergence and stability can be found in [2,9]. </p><p>* Corresponding author. </p><p>0045-7825/97/$17.00 C3 1997 Elsevier Science S.A. All rights reserved PZZ SOO45-7825(96)01170-X </p></li><li><p>94 B.A. Schrefler et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>If the equations obtained from an operator split are solved simultaneously we have on the other hand an inherent or physical parullelization of the problem. This means that the discretized equations of every single field can be solved on a separate processor [lo]. However, the splitted variables are not a priori labelled, hence the results of the analysis presented here (in particular the error formulas) are also valid for domain decomposition techniques. </p><p>The idea of partitioning has been associated with fractional step methods and in this context analyzed [ll]. The fractional step method has been introduced in Russian literature (e.g. see [12,13]). It is a scheme which reduces a complicated problem into a sequence of simpler ones by splitting-up (partitioning) the operator of the solved equations. In this way we obtain either differential equations of smaller order or of simpler form. To preserve the convergence and stability properties of the full operator the component operators must maintain appropriate group characteristics (see [13,14,16]). </p><p>As a modification of the standard staggered procedure, the Newton procedure was suggested. It has been formulated and implemented for thermo-mechanical problems in [15,17-191, for which the coupling is weak. However, also in this case the possibility of applying such algorithms relies on heuristic attempts only [19], without a correct analysis of all the numerical properties of the solution algorithm, not only stability [5-7). </p><p>In the paper we compare the iteration properties of both algorithms with the aim of checking the applicability and the convenience of the staggered Newton method to consolidation problems, which are usually strongly coupled and characterized by non-symmetric matrices appearing in the semi- discretized equations. From this point of view we also investigate the consequences of a common practice in simulating slow phenomena, which consists in increasing the time step length according to the slowing down of the process. </p><p>2. The model problem </p><p>For the sake of completeness, we summarize here the equations governing the thermo-hydro- mechanical problem of consolidation for a multicomponent medium with a solid phase and two fluid phases (water and gas) [20]. In the following, the index w refers to the generic phase, whereas s, 1 and g represent solid, liquid water and gas, respectively: </p><p>- Linear momentum balance equation for the volume fraction mixture (neglecting inertial effects) in terms of total stresses (T: </p><p>V.a+pb=O (I) </p><p>where p is the average density of the mixture. </p><p>P = (I- 4)P, + 4% + WPg </p><p>C$ is the porosity, p,, is the density force; </p><p>- Dry air conservation equation: </p><p>(2) </p><p>of the m-phase, S the degree of saturation and b the specific body </p><p>(3) </p><p>where ps,, is the mass concentration of dry air in gas phase, u the displacement vector in solid matrix, ug the velocity of the gas phase and u&amp;, the relative average diffusion velocity of dry air species. </p></li><li><p>B.A. Schrefler it al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 95 </p><p>- Water species (liquid and vapor) conservation equation: </p><p>4 $ [(I - W,:,l + (1 - ~)P*w $ (V*u) + v* (P&amp;) + v* (Pg&amp;) </p><p>=-~~l~-sP,$(v.u)+v.(P,u,) (4) </p><p>where pgw is the mass concentration of water vapour in gas phase, ut the velocity of the liquid phase; </p><p>- Energy conservfation (enthalpy balance): </p><p>pep at </p></li><li><p>96 B.A. Schrefler et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>(7) </p><p>, Y ntl -Y, </p><p>n+l = At </p><p>This is a particular case of the implicit mid-point rule for 8 = 1. The use of this temporal integration is particularly useful in convective problems, however the generalisation is straightforward and used in the applications. </p><p>Thus, from Eq. (6) in conjunction with Eqs. (7), (8) we find </p><p>(I+ AtC,+i)%+i = AtF,,+i +x,, (9) </p><p>Now we shall focus our attention on solving Eq. (9) for a fixed time step n + 1. Let C,,,, = [c,], then Eq. (9) in matrix form is written as </p><p>1+ Ate,, </p><p>At c21 AtC12 ][xn+] = At[;+] + [;;I </p><p>I + Ate,, y,+l (10) </p><p>To perform the staggered iteration scheme we first split the matrix on the 1.h.s. of Eq. (10) </p><p>[ </p><p>1 + At cl1 At% l+Atc,, 0 </p><p>At ~21 l+Atc,, = 0 1 [ 1 + At c22 ] +At[c:l c;2] (11) then transfer part of it to the r.h.s. The iteration procedure starts with a predicted value of y,, 1,0, x,+ 1,0 (only with yn+l,o in the linear case) and provides new estimates x,+ l,i and y,, l,i by the following equations: </p><p>X n+l,i = (I+ Atcll)-[AtFZi+l +Xn,K -ClZYn+l,i-11 (12) </p><p>where </p><p>Cl1 =cll(xn+l,i-l~ Yn+l,i-1 )9 Cl2 = dXn+l,i-1, Yn+l,i-1 19 </p><p>Y n+l,i = (1 + AtC22)-[AtY,+l +Yn,k -CZIX~+I,II (13) </p><p>where </p><p>22 =c22(xn+l,iT Yn+l,i-1) 7 c21 = C21(Xn+l,i~ Yn+l,i-1) . </p><p>K denotes the number of performed iterations and depends on n, K = K(n). The predictor usually is chosen as a linear extrapolation of the last obtained values, although to </p><p>obtain better convergence properties it is recommended to use an explicit one or multi-step formula [9]. In concise form the above equations (12), (13) can be written as </p><p>X n+l,i =p(xn+l,i-l~ Yn+l,i-1) (14) </p><p>Y n+l,i = Q</p></li><li><p>If P, </p><p>B.A. Schrefler et al. I Comput. Methods Appl. Mech. Engrg. 14.4 (1997) B-109 </p><p>Q satisfy the necessary conditions then by Taylors theorem </p><p>~(~i_l, yi_l&gt; q = P(x, y) + P,x(x, y)Ef-, + P,y(x, Y)E~-, + O((E~-,, Er-,)2) </p><p>Q(Xi&gt; yi-1) = !Z(X, y) + Q,,(x&gt; Y)E~ + Q,, Y,+1) = 0 </p><p>for the monolithic problem Eq. (9) written as </p><p>(25) </p><p>(26) </p></li><li><p>98 B.A. SchrefIer et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>46 n+l, Yn+l) =o (27) where </p><p>P(X n+l&gt; y,+l) = (I+ Atc,,)~,., + A~GY,+~ --xx - AtC+, (28) </p><p>4(x n+l&gt; Y,+J = Afcx~,z+I + (I+ A~c,,)Y,+, -Y, - AtFY,+I (29) </p><p>As in the previous section, to simplify the subsequent writing, we shall omit subscript IZ + 1: x n+l,i =i, Yn+l,i =Yi The staggered Newton scheme is given by the iterations </p><p>Xi = xi_* -Ptxi-l? Yi-l)P,x(i-19 Yi-1) (30) </p><p>Yi =Yi-1 - 4txi? Yi-1)4,y(xi7 Yi-1) (31) </p><p>where x,,, y, are the starting values. Expanding the quotients in Eqs. (30) and (31) by Taylor series in point (x, y) and taking into account </p><p>that X, y are exact solutions of Eqs. (30), (31) we obtain </p><p>E; = - 2 E;_, - O((E:_,, E;__,)*) (32) ,x </p><p>ET = -p E; - O((El, E;__,)*) (33) .Y </p><p>where all the derivatives are taken in point (x, y). In matrix form </p><p>Ei =ANEi_, </p><p>where </p><p>(34) </p><p>is the amplification matrix. Condition </p><p>p(AN) &lt; I </p><p>must be satisfied to ensure convergence. </p><p>(35) </p><p>As we see, the staggered Newton scheme is only of first order-it has the same order as the standard staggered scheme but requires an additional evaluation of tangent operators. </p><p>For weakly coupled problems, i.e. P,~ = 0, the method is of second order. This was the case considered by [ 15,18,19], where coupling terms were neglected and good numerical properties achieved. </p><p>The comparison of the maximum eigenvalues of the amplification matrices AS and AN tells us which iteration should proceed faster, i.e. the smaller the spectral radius the faster the process. This is concluded from the fact that for any matrix norm I] - (1 compatible with a vector norm: </p><p>&amp; ]]A]] =,1ic ~(4 (36) </p><p>The eigenvalues of matrix AS are </p><p>G.2 = + @,P,Y + Q,y + 9x + 4 (37) </p><p>where </p><p>A = KQ,xf,y + Q.y + &amp;I* - 4f,,Q,ylo.5 </p></li><li><p>B.A. Schrefier et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 99 </p><p>Matrix AN possesses eigenvalues </p><p>(3% </p><p>To be able to compare them we must find a relation between P, Q and p, q. From the definitions of P, Q (found in Eqs. (14), (15)) in conjunction with Eqs. (30) (31) we obtain </p><p>P(X, Y) = (I+ Atc,,)x- (1 +Atc,,)% y) (39) </p><p>dx, Y) = Cl+ At 4~ - Cl+ At 4Qk Y&gt; (40) </p><p>Taking into account that x, y is the exact solution of Eqs. (14), (15) and Eqs. (30), (31) we find that </p><p>P,,(K Y) = (1-t At cii)(l - P&amp;, Y)) (41) </p><p>p,,(x, y) = -Cl + At cn)Sy(x~ Y&gt; (42) </p><p>q,x(-v Y&gt; = -Cl + At 4Q,xk Y) (43) </p><p>q,,b Y) = (1 + At c,,)(l - Q,,b Y&gt;) (9 </p><p>Let us consider some particular cases: l Weakly coulpled problems, i.e. P,~ = 0, q,xO then A,,* = P,x; Q,y and AN = 0. The staggered </p><p>Newton is of second order. l Alternatively coupled problems, i.e. p,, = 0, q,y = 0 then the staggered Newton cannot be used. l Strongly coupled problems in the sense that P,X = 0, Q,Y = 0 then </p><p>hf=O; As = P,,Q,, and </p><p>A;=o; A: =P,yq,J(l + At eii)(l + Ate,,) = 9,Q.x </p><p>thus the spectral radii are equal, both procedures converge at the same rate. l In the case of P,x = 0, the standard staggered procedure may be faster than the Newton one. </p><p>Then Ai;, = 0; P,,Q,, + Q,y and A:z = 0; P,,Q,,/(l - Q,,). A typical example is the system of equations in [2, Appendix C]. </p><p>l The staggered Newton scheme can be divergent, even for starting values very close to the roots. </p><p>5. Application to consolidation problems </p><p>We now apply the standard staggered and staggered Newton procedures to the solution of a partially saturated soil behaviour. We limit the application to physical partitioning. The strength of coupling between the solid and fluid fields is stronger than the coupling with the temperature field [3], hence to analyse the different numerical procedures, we do not take into account the last one. In particular, the analysis deals with the problem of desaturation of a sand column due to gravity, solved for instance in [21]. This problem is certainly a strongly coupled one during the initial part of the desaturation, i.e. the physically observed behaviour can only be obtained by taking the full coupling into account [22]. Furthermore, the strength of coupling results also from the iteration count [23]. Purely thermo- mechanical, such as in [ll], exhibit generally weak coupling. We consider that air remains always at atmospheric pressure, hence the model consists of an equilibrium equation for the multi-phase medium and of a continuity equation for water (Eqs. (1) and (3)). The equations are defined in a closed domain R = [0, h], with h variable so as to be able to analyse the effect of spatial discretization (h was set equal to 50 an.d 100 cm). </p></li><li><p>100 B.A. Schrefer et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>8.M3E-01 - i \ </p><p>7.WE-01 -- I I </p><p>&amp;WE-01 </p><p>.s SXlE-01 </p><p>B </p><p>P 4LoE-01 </p><p>ii 3.00E-01 u) </p><p>-.-**-. Suc.Subst.Ston. </p><p>- N.- R. Stag. </p><p>---- suc.suixt.Stag. </p><p>2.CKlE-01 </p><p>1.00E-01 </p><p>0.coE+cxl I </p><p>1 </p><p>.:.\, </p><p>20 70 300 800 4cra pooo 14oxl 19wl 240Yl </p><p>time (sect) </p><p>Fig. 1. Comparison between spectral radii: Case 1. </p><p>7.aIE-01 </p><p>6.ooE-01 - N.- R.Stag. </p><p>5.ME-O1 </p><p>4.0JE-01 </p><p>3.COE-01 </p><p>2.CQE-01 </p><p>O.ooE+OO </p><p>1 6 20 70 300 8Ou 4000 QXO 14cca 19coO 24Oa </p><p>Ime @ed </p><p>Fig. 2. Comparison between spectral radii: Case 2. </p><p>35 </p><p>--.-- N.-R.Stan. </p><p>---- sw.suw.stag. </p><p>- N. -R. stacl. </p><p>5 t </p><p>04 </p><p>1 6 20 70 3cn eal 4ooo wxx) 14lX0 19Wl 24cal </p><p>Ime Met) </p><p>Fig. 3. Comparison between number of iterations per time step in different solution procedures: Case 1 </p></li><li><p>B.A. Schrejkr et al. I Comput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 101 </p><p>xl -- </p><p>1! i 40-- </p><p>3 30 -- </p><p>% __ 2 </p><p>10 __ f </p><p>----- N.-R.Stan. </p><p>---- suc.subst.stag. </p><p>- N.- R Stag. </p><p>-----__-..- </p><p>01 I </p><p>1 6 20 70 303 800 4ax Km 14CKO 19fYXl 24oml </p><p>lime (see) </p><p>Fig. 4. Comparison between number of iterations per time step in different solution procedures: Case 2. </p><p>3.COE-02 </p><p>2.5OE-02 </p><p>! 200E-02 3 1.5OE-02 </p><p>i l.DlE-02 </p><p>5.oOE-03 </p><p>O.COE+M) </p><p>..---.. succ.subst.stan. </p><p>- - - - succ. sub&amp; stag. </p><p>----N.-R.Stan. </p><p>- N.- R. Stag. </p><p>6 20 70 3al 800 4am woo 14Wl 19KO 24Krl tlme (SW) </p><p>Fig. 5. Maximum error at the beginning of each time step: Case 1. </p><p>1.4OE-02 I I </p><p>1.2OE-02 </p><p>l.KW2 </p><p>0 8 8.C.QE-03 E </p><p>i </p><p>1 mE-03 4.IXE-03 </p><p>2.mE-03 </p><p>O.WE+W </p><p>---- suc.suhlt.stag. - -. -. N.- R. Stan. </p><p>- N.-R. Stag. </p><p>6 20 70 3al 800 4aYl m 14am 190x 24wO </p><p>time (ted Fig. 6. Maximum error at the beginning of each time step: Case 2. </p></li><li><p>102 B.A. Schrefler et al. I Cornput. Methods Appl. Mech. Engrg. 144 (1997) 93-109 </p><p>The test consists in solving a 1-D model problem with only 2 degrees of freedom. Dirichlets boundary conditions at x = h are associated with zero values of variables, while at h = 0 natural boundary conditions with zero force and flow are assigned. The forcing function is the gravitational ef...</p></li></ul>


View more >