PARTITIONED SOLUTION PROCEDURE F

Document Sample

```					COMMUNICATIONS IN NUMERICAL METHODS IN ENGINEERING,                                  Vol. 13: 239±247 (1997)

PARTITIONED SOLUTION PROCEDURE FOR
SIMULTANEOUS INTEGRATION OF
COUPLED-FIELD PROBLEMS

JEAN H. PREVOST
Department of Civil Engineering and Operations Research, Princeton University, Princeton, NJ 08544, U.S.A.

SUMMARY
An implicit unconditionally stable partitioned solution procedure for the simultaneous integration of
transient coupled ®eld problems is presented. The procedure does not require that the full system of coupled
equations be assembled, and allows use of existing single-®eld analysis software modules to solve the
coupled ®eld problem. An iterative partitioned conjugate gradient procedure is used to avoid having to form
and assemble the Schur complement matrix. The coupling matrices never need be formed, thus resulting in
substantial computational savings. # 1997 by John Wiley & Sons, Ltd. Commun. numer. meth. engng.,
vol. 13, 239±247 (1997).
(No. of Figures: 2 No. of Tables: 1 No. of References: 20)

1. INTRODUCTION
The transient response analysis of coupled-®eld problems is becoming a major challenge in
many engineering disciplines. In coupled problems two or more physical entities interact with
each other, with the independent solution of any one entity being impossible without simul-
taneous solution of the others. Some problems, such as dynamic ¯uid-structure interaction,
may be characterized by distinct spatial subdomains (e.g. ¯uid and structure) with distinct
sets of dependent variables (e.g. ¯uid pressure and solid displacements). In that case, neither
subsystem can be solved independently of the other due to the unknown interface inter-
action forces. Other problems are characterized by the interaction of essentially dierent
physical phenomena within the same material domain, e.g. thermal and mechanical interaction.
Although it is believed that the proposed procedure could also be useful for the ®rst class of
problems, it is this last class of problems which is of primary concern in this paper. In this
class of problems the coupling occurs through the dierential governing equations describing
the dierent physical phenomena. Each component ®eld is then spatially discretized into a
®nite number of degrees of freedom through a ®nite element method and called a subsystem.
The goal is then to integrate in time the resulting coupled semi-discrete matrix equations of
motion.
A typical system of equations which arises in linear thermo-elasticity is given as
!       !                   !       !                !       !            !
M      0       
u             C     0       
u           K   ÀL       u           fu
                                                                            1
0      0       y           y0 LT   S       y           0   H        y           fy

CCC 1069±8299/97/040239±09\$17.50                                                                     Received 1 May 1996
# 1997 by John Wiley & Sons, Ltd.                                                                Revised 1 December 1996
240                                             J. H. PREVOST

where u  nodal displacement vector ®eld, y  nodal temperature ®eld; M, C and K are the
mass, damping and stiness matrices; f u is the vector of prescribed structural forces; S, H and L
are the heat capacitance, heat diusion and thermal expansion coupling matrices; f y is the vector
of external heat sources; and y0  reference temperature. All the matrices on the diagonals are
symmetric but the overall matrices are not. A similar system of equations arises in saturated
porous media analysis when describing the interaction of the pore-¯uid pressure with the porous
solid soil skeleton. Direct simultaneous time integrations of the coupled semi-discrete equations
have been used,1±4 and the resulting algorithms have been shown to be unconditionally stable.
However, such implementations have several limitations:5 (i) there is a need to develop special
software modules which combine the coupled ®eld equations; (ii) large matrix problems result,
especially for three-dimensional cases. For these reasons, staggered implementations have been
suggested,5±11 by which each ®eld equations are solved separately assuming that the ®eld
variables of the other subsystems are known (via a predictor) and temporarily frozen. There are
many advantages to such staggered procedures: (i) modularity features which allow the coupled
equations to be processed by separate program modules taking full advantage of specialized
features and disciplinary expertises built into independently developed single-®eld analysers;
(ii) the resulting algorithmic structure, which allows the set of analysers to be synchronized to
operate in sequential or parallel fashion. However, simple, straightforward implementations
of staggered procedures are known to be at best only conditionally stable even though
unconditionally stable implicit integrators are used. Various stabilization procedures have
been proposed,5Y6Y9Y11 but they entail costly restructuring and augmentation of the matrix
problem, e.g. by requiring that a stabilization matrix of the form LT M À 1 L be formed and
operated with.6
In the present work, we propose an unconditionally stable, robust and computationally
ecient partitioned procedure for the simultaneous integration of transient coupled ®eld
problems. The procedure allows the ®eld-state vectors of the coupled equations to be processed
by separate discipline-oriented program modules. In Section 2 we brie¯y review the basic
equations for coupled-®eld analysis. In Section 3 we introduce the partitioned solution proce-
dure. Numerical experiments in Section 4 illustrate the proposed agorithm. Finally, conclusions
are drawn in Section 5.

2. COUPLED SYSTEMS
2.1. Governing equations
For simplicity in the presentation we consider a two-®eld coupled system modelled after thermo-
elasticity and governed by the following semi-discrete equations of motion:


M 1 À Lu2  n1 u1 Y u1 Y u2   f 1
u
S 2  y0 LT u1  n2 u2 Y u2 Y u1   f 2
u                                                                2

where u1 and u2 denote the two coupled ®eld-state vectors (e.g. u1  solid displacement;
u2  temperature; a superposed dot is used to indicate time dierentiation; n1 and n2  non-
linear vector-valued functions of the ®eld-state vector and time derivatives; f 1 and f 2  vectors
representing the external loads. The global quantities M, L, etc. are assembled from the
corresponding ®nite element matrices and vectors.

# 1997 by John Wiley & Sons, Ltd.                      COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
SIMULTANEOUS INTEGRATION OF COUPLED-FIELD PROBLEMS                                   241

2.2. Time discretization
The semi-discrete ®nite-element equations are a system of coupled ordinary dierential equations
which are to be integrated in time. The system is solved by applying a step-by-step integration
procedure, and a direct time integration scheme is constructed to advance both ®eld variables
simultaneously in time. Various algorithms are available for that purpose (see, for example,
References 12 and 13). In the following, simple one-step methods are used to carry out the time
integration of both the second-order and ®rst-order systems by assuming auxiliary relationships
between the solution vectors at time steps tn and tn  1 . For the second-order system, the Newmark
update formulas14 are used, namely
Dt 2

u1Yn  1  u1Yn  Dt u1Yn                                 
1 À 2b1  1Yn  2b1 u1Yn 1 
u
2
          
u1Yn  1  u1Yn                            
 Dt1 À a1  1Yn  a1 u1Yn  1 
u                                                 3

                                                     
where Dt  tn 1 À tn ; and u1Yn , u1Yn and u1Yn are approximations for u1 t, u1 t and u1 t,
respectively at time tn . The scalars a1 and b1 are algorithmic parameters which can be chosen to
ensure unconditional stability (namely 1a2 4 a1 4 3a2Y aa2 4 b1 4 1 and second-order accuracy
(namely a1  1a2 in the linear case. In the following, a1 5 1a2 and b1  a1  12 a4,15 to
2

maximize high-frequency numerical dissipation.
Similarly, the generalized trapezoidal family of ®nite dierence time stepping algorithms12Y13 is
used for the ®rst-order system, namely

u2Yn  1  u2Yn  Dt u2Yn  a
                              
u2Yn  a  1 À a2  2Yn  a2 u2Yn  1
u                                               4

where a2 is an algorithmic parameter which can be chosen to ensure unconditional stability
(i.e. 1 4 a2 4 1 and second-order accuracy (i.e. a2  1 in the linear case. Maximal high-
2                                                2

frequency dissipation is provided by selecting a2  1.

2.3. Non-linear iterations
The resulting system (2±4) is a non-linear system of equations in general. The solution to this
system is obtained by using an iterative strategy which is implemented by means of a predictor±
multicorrector scheme applied at each time step. In this method, a series of corrected solutions
are computed starting from an initial or predicted solution for the time step. Each corrected
solution is used in the following iteration to compute the next corrected solution. The procedure
continues until a speci®ed solution convergence has been obtained. There are several ways of
implementing the recursive relationships that take the solution from step n to step n  1. In our
implementation we always chose to use the vector of nodal quantities with the highest time
derivatives as the vector of primary unknowns, namely for the model problem we use the vectors
            
u1Yn  1 and u2Yn  1 as the primary unknown variables. The resulting system of equations is solved
by performing a linearization via a truncated Taylor's series expansion, and results in the
following coupled linear system of equations to be solved at each iteration:
!4 i        5 4           5
M  a1 DtC  b1 Dt 2 K   Àa2 Dt L       u
D 1Yn 1      ri  1
1Yn
                        5
y0 LT         S  a2 Dt H D i 1
u2Yn          ri  1
2Yn

# 1997 by John Wiley & Sons, Ltd.                         COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
242                                                 J. H. PREVOST

where the superscript index i is used to index the non-linear iterations; r1 and r2 denote the
residuals, namely

ri  1  f 1Yn  1 À M i  1  Lui  1 À n1 ui 1 Y ui  1 Y ui  1 
1Yn                    u1Yn          2Yn            1Yn       1Yn       2Yn

ri  1  f 2Yn  1 À S i  1 À y0 LT ui 1 À n2 ui  1 Y ui 1 Y ui 1 
2Yn                    u2Yn               1Yn          2Yn        2Yn 1Yn             6

which can be evaluated from initial and known values from the previous iterate. The corrector
formulas are written as (from (3) and (4))

ui  11  ui  1  D i 1
 1Yn       1Yn       u1Yn
ui  11  uÁ  1  a1 Dt D i  1
 1Yn       1Yn             u1Yn
ui  11  ui  1  b1 Dt 2 D i  1
1Yn        1Yn                u1Yn                                 7

and

ui 11  ui 1  D i 1
 2Yn      2Yn      u2Yn
ui 11  ui 1  a2 Dt D i 1
2Yn       2Yn             u2Yn                                    8

The correctors, (7) and (8), ensure that every set of iterates adheres to the recursive formulas, (3)
and (4). The matrices appearing in (5) are Jacobian-like matrices, namely

2            dr1 
M a1 DtC  b1 Dt K  À
d 1  i Yun YDt
u u
1Yn  1

dr2 
S a2 Dt H  À                                                       9
d 2  i Yun YDt
u u
1Yn  1

     
where un  fu1Yn Y u2Yn Y u1Yn Y u2Yn , etc.} denotes the collection of known parameters. Note that the
coupling matrices appearing in (5) are `incomplete' Jacobian matrices, since in general
                             
n1  n1 u1 Y u1 Y u2  and n2  n2 u2 Y u2 Y u1 .
Remark 1: The resulting coupled matrix equation in (5) may be symmetrized by multiplying
the second equation by Àa2 Dtay0.
Remark 2: In the linear case, (2) has the form of (1), and (5) needs to be solved only once at
every time step (no iterations are required), and only needs to be re-evaluated when the time
step value changes.
Remark 3: It is easily shown by standard procedures13Y16 that unconditional stability of the
algorithm (in the linear case) is assured when unconditionally stable implicit integrators are
used to integrate each ®eld equation.2Y4
Remark 4: Consistency and accuracy of the combined solution is inherited directly from the
individual integrators used for each ®eld equation.

3. PARTITIONED SOLUTION PROCEDURE
At the core of the combined simultaneous solution procedure is the solution of the set of linear
equations (5). If the full matrix system is to be used, the assembly and factorization of the implicit

# 1997 by John Wiley & Sons, Ltd.                          COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
SIMULTANEOUS INTEGRATION OF COUPLED-FIELD PROBLEMS                                 243

coecient matrix is found to pose enormous computational demands for three-dimensional
problems. It also requires that special software modules be developed to combine the ®eld
equations, hampering the direct use of discipline-oriented software modules for each ®eld
equation. To preserve software modularity and to avoid having to operate with the full matrix
equation, we solve (5) by using a partitioned iterative conjugate gradient procedure17 as follows.
For ease of presentation we rewrite (5) as
Ax  b                                           10
where A is the coecient matrix appearing in the left-hand side of (5); x is the vector of
unknowns; b is the right-hand side. Further, we partition the matrix problem as follows:
!    !         !
A11 A12 x1           b1
                                   11
A21 A22 x2           b2
where x1  D i and x2  D i . A direct partitioned solution of (11) in general requires that the
u1               u2
Schur complement matrix be formed and assembled. In order to circumvent that diculty, we use
an iterative partitioned conjugate gradient procedure as shown in Table I.
In Table I, dtol  convergence tolerance, B a symmetric positive-de®nite conditioner matrix
and x0  initial guess to the solution. In our implementation we have used B  A22 and
2
x0  0 successfully. The computations of the matrices A11 and A22 pose no diculty and are
2
handled by their respective discipline-oriented software modules developed for the single-®eld
problems. The coupling matrices A12 and A21 are never to be formed and assembled for the
procedure to retain its versatility and modularity. The matrices are only needed to compute
matrix±vector products (Ap-type products) in the algorithm. Since the matrix A is a Jacobian-
like matrix, it is possible to de®ne residual-like vectors such that:
A12 p2  À lim ~1 u1  ep2  À ~1 u1 ae
r                r
e30
A21 p1  À lim ~2 u2  ep1  À ~2 u2 ae
r                r                                  12
e30

for given solution vectors u1 and u2 . In this implementation, evaluation and storage of the
components of A12 and A21 are avoided altogether, and an explicit ®rst-order-accurate forward
central dierencing scheme is adopted to approximate (12) as
A12 p2  À ~1 u1  ep2  À ~1 u1 a e
r        "       r          "
A21 p1  À ~2 u2  ep1  À ~2 u2 a e
r        "       r          "                          13
where e  e1a2 a k p k and em  machine precision. We have also experimented with a second-
"   m
order-accurate central dierence scheme (twice as expensive), but noted no signi®cant
improvement in the performance of the CG iterations.
Remark 1: The iterative partitioned conjugate gradient algorithm shown in Table I is used
to solve for the unknown x2 , namely if one eliminates x1 in (11), then
~
Sx2  ~2
b                                         14

with
~2  b2 À A21 A À 1 b1
b                                                      15
11

# 1997 by John Wiley & Sons, Ltd.                       COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
244                                              J. H. PREVOST

Table I. Flow chart for conjugate gradient

Initialize:
x0
2

x0  A11 b1 À A12 x0 
1
À1
2

r0  b2 À A21 x0 À A22 x0
2               1          2

z0  B À 1 r0
2            2

p0  z0
2      2

Iterate:
do for k  0Y 1Y F F F Y
pk  ÀA11 A12 pk
1
À1
2

Apk  A21 pk  A22 pk
2          1          2

ak  rkT zk apkT Apk
2     2     2        2
k
x1 1  xk  ak pk
1           1
k
x2 1  xk  ak pk
2           2
k
r2 1  rk À ak Apk
2              2
k
Convergence check: if  k r2 1 k 4 dtol k r0 k  end doloop
2
k            k
z2 1  B À1 r2 1
k       k
bk  r2  1T z2  1 arkT zk
2     2
k       k
p2 1  z2 1  bk pk
2

end do

and where
~              À
S  A22 À A21 A111 A12                                  16

is the Schur complement. Note that the Schur complement matrix is never formed or
assembled in the implementation shown in Table I. The success of the partitioned iterative
À
procedure used to solve for the Schur complement hinges upon A22 and A21 A111 A12 being
18
symmetric positive-de®nite matrices. This is ensured since both A11 and A22 are symmetric
positive-de®nite matrices (5), and A12 and A21 are transposes of each other (up to a scalar
multiplier).
Remark 2: Computing the matrix±vector product using ®nite dierences requires working
with at least 64-bit word lengths em % O10 À 16 . Consequently, working in double
precision on a 32 bit machine is mandatory.
Remark 3: Clearly, the ordering 1±2 of the ®eld-state vectors is arbitrary. In many cases, the
dependent variables are not of the same order (for example, displacement vector and scalar
temperature). In that case, computational savings can be achieved by selecting for partition 2
the set of variables of the lowest order (namely temperature).

# 1997 by John Wiley & Sons, Ltd.                      COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
SIMULTANEOUS INTEGRATION OF COUPLED-FIELD PROBLEMS                         245

Figure 1.

4. NUMERICAL EXAMPLES
To illustrate the proposed procedure, we consider the consolidation of a ¯uid saturated poro-
elastic sphere subjected to a pressure step loading with uniform intensity p on its outer surface at
time t  0  . The ®nite element mesh and problem geometry are shown in Figure 1. In this
calculation, 4-node axisymmetric bilinear isoparametric elements are used. In order to
accommodate the near incompressible isochoric response of the solid porous matrix in the
early stage of the consolidation process, the mean-dilatation formulation12 is used to integrate the
dilatational component of the solid response.
The ®nite element mesh consists of 73 nodes and 60 bilinear axisymmetric quads. The
following material properties are used: rs  solid mass density  2Á0; r2  fluid mass
density  1Á0; nw  porosity  0Á3; E s  solid Young's modulus  104 ; ns  Poisson's
ratio  0Á0; k  permeability  1Á6 Â 10 À 3 . Figure 2 shows the computed radial contraction
of the sphere d and pore-¯uid pressure pw at the centre of the sphere as a function of time
T  C v taR2 , where Cv  karw ls  2ms   consolidation coecient; ls and ms  solid Lame     Â
constants (in this case ls  2ms  E s ; and T  dimensionless time. The results are presented as
ratios with the ultimate t  I radial contraction dI  pRa3ls  2ms   4Á0 Â 10 À 4 and with
the initial t  0   pore-¯uid pressure pw 0  p. A closed form analytical solution is available
for this problem,19Y20 and numerical agreement with the analytical solution is shown in Fig. 2.
The numerical solution accurately predicts the increase in pore-¯uid pressure which occurs at the
early times of the consolidation process. This phenomenon is characteristic of the coupled
theory and is referred to as the Mandel±Cryer eect. The numerical results shown in Figure 2
were obtained by using 20 time steps with a constant time step Dt  0Á025, and by taking
a1  a2  1Á0. For a convergence tolerance of 10 À 4 , a maximum of seven conjugate gradient
iterations were needed to compute the coupled solution at any time step. A closer match with the
exact solution can be achieved by using a smaller time step (especially in the early stages of the
diusion process).

# 1997 by John Wiley & Sons, Ltd.               COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
246                                            J. H. PREVOST

Figure 2.

5. CONCLUSIONS
An implicit unconditionally stable partitioned solution procedure for the simultaneous
integration of transient coupled ®eld problems is presented. The procedure does not require
that the full system of coupled equations be assembled, and allows use of existing single-®eld
analysis software modules to solve the coupled ®eld problem. An iterative partitioned conjugate
gradient procedure is used to avoid having to form and assemble the Schur complement matrix.
The coupling matrices never need to be formed, thus resulting in substantial computational
savings.

REFERENCES

1. J. H. Prevost, `Consolidation of anelastic porous media', J. Eng. Mech., ASCE, 107, EM1 169±186
(1981).
2. J. H. Prevost, `Implicit-explicit schemes for nonlinear consolidation', Comput. Methods Appl. Mech.
Eng., 39, 225±239 (1983).
3. J. H. Prevost and D. Tao, `Finite element analysis of dynamic coupled thermoelasticity problems with
relaxation times', J. Appl. Mech., ASME, 50, 42, 817±822 (1983).
4. O. C. Zienkiewicz and R. L. Taylor, `Coupled problems Ð A simple time-stepping procedure',
Commun. appl. numer. methods, 1, 233±239 (1985).
5. K. C. Park and C. A. Felippa, `Partitioned analysis of coupled systems', in T. Belytschko and T. J. R.
Hughes (Eds.), Computational Methods for Transient Analysis, North-Holland, Amsterdam, 1983,
pp. 158±219.
6. C. Farhat, K. C. Park and Yves DuBois-Pelerin, `An unconditionally stable staggered algorithm for
transient ®nite element analysis of couple thermoelastic problems', Comput. Methods Appl. Mech. Eng.,
85, 349±365 (1991).
7. C. A. Felippa and K. C. Park, `Staggered transient analysis procedures for coupled mechanics problems:
Formulation', Comput. Methods Appl. Mech. Eng., 24, 61±111 (1980).
8. K. C. Park, `Partitioned transient analysis procedures for coupled ®eld problems: Stability analysis',
J. Appl. Mech., 47, 370±376 (1980).
9. K. C. Park, `Stabilization of partitioned solution procedure for pore ¯uid-soil interaction analysis', Int.
J. Numer. Methods Eng., 19, 1669±1673 (1983).

# 1997 by John Wiley & Sons, Ltd.                    COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)
SIMULTANEOUS INTEGRATION OF COUPLED-FIELD PROBLEMS                                247

10. K. C. Park and C. A. Felippa, `Partitioned transient analysis procedures for coupled ®eld problems:
Accuracy analysis', J. Appl. Mech., 47, 919±926 (1980).
11. O. C. Zienkiewifcz, D. K. Paul and A. H. C. Chan, `Unconditionally stable staggered solution
procedure for soil-pore ¯uid interaction problems', Int. J. Numer. Methods Eng., 26, 1039±1055 (1988).
12. T. J. R. Hughes, The Finite Element Method, Prentice Hall, Englewood Clis, NJ, 1987.
13. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol. 2, McGraw-Hill, London, 1991.
14. N. M. Newmark, `A method of computation for structural dynamics', J. Eng. Mech. Div., ASCE, 85,
EM3, 67±94 (1959).
15. H. M. Hilber, `Analysis and design of numerical integration methods in structural dynamics',
Rep. No. EERC 76-29, Earthquake Engineering Research Center, University of California, Berkeley,
CA, 1976.
16. T. J. R. Hughes, K. S. Pister and R. L. Taylor, `Implicit-explicit ®nite elements in nonlinear transient
analysis, Comput. Methods Appl. Mech. Eng., 17/18, 159±182 (1979).
17. M. R. Hestnes and E. Stiefel, `Method of conjugate gradients for solving linear systems', J. Res. Nat.
Bus. Standards, 49, 409±436 (1952).
18. G. Golub and J. M. Ortega, Scienti®c Computing: An Introduction with Parallel Computing, Academic
Press, 1993.
19. C. W. Cryer, `A comparison of the three-dimensional consolidation theories of Biot and Terzaghi',
Q. J. Mech. Appl. Math., 16, 401±412 (1963).
20 J. Mandel, `Consolidation des sols (etude mathematique)', Geotechnique, 3, 287±299 (1953).

# 1997 by John Wiley & Sons, Ltd.                   COMMUN. NUMER. METH. ENGNG., VOL. 13, 239±247 (1997)

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 5 posted: 12/17/2009 language: English pages: 9
How are you planning on using Docstoc?