Newton Scheme and Monolithic Approach

φµ^f

k δv^f_h,−τMR^f,h_M

Ωe

. (4.138)

• Biot stabilization

Further, a special stabilization term for the porous medium system (also called Biot system in linear theory) is included. It becomes relevant, when low permeabilities and dynamic effects are considered. More precisely, when the time derivative of the porosity in the continuity equation is non-zero, unphysical stress and pressure oscillations occur for low-permeability problems. From an engineering perspective, this effect can be interpreted as volumetric locking of the porous system [189]. One possible remedy is the introduc-tion of another stabilizaintroduc-tion term to the balance of mass of the fluid phase [18]. The fol-lowing contribution is added to the continuity equation of the fluid phase (first equation of (4.66), (4.80) and (4.95))

nele

X

e=1

∇δp^f_h, τBR^s,h_M

Ωe

. (4.139)

PSPG and the stabilization due to the reactive term are used within all numerical simulations involving a porous medium. As stated before, the Biot stabilization is only necessary in certain kinds of settings. The Biot stabilization is used and analyzed in the example in Section 4.4.1.

4.3.4. Newton Scheme and Monolithic Approach

After discretization of (4.66), (4.80) or (4.95) in time and space the discrete residuals R^s = R^s(d^s,v^f)andR^f =R^f(d^s,v^f)of the skeleton and fluid equations are obtained. In this context,d^s andv^f denote the vectors of nodal displacements (and porosities in case of the mixed approach) and fluid unknowns (i.e., fluid velocities and pressures), respectively. The residuals consist of the discrete form (4.116), (4.80) or (4.95) and the stabilization terms presented in Section 4.3.3.5.

For the solution of this system of non-linear equations at a time stepn+ 1

a Newton-Raphson method is utilized. This fully coupled non-linear problem will be solved within one Newton loop. A so-called monolithic solution scheme is applied, wherein the pri-mary variables of the respective fields will be solved simultaneously. In comparison to parti-tioned schemes, where an iteration between single field solvers is utilized, this approach has been successfully applied to other multi-field applications and proven to be superior especially for complex biological problems (see, [107, 151]). The corresponding linearized equations at a iteration stepiare given as

LinR=R(d^s,i_n+1,v^f,i_n+1) + ∂R(d^s,v^f) The most important linearizations can be found in Appendix A.2. This system is solved for the unknown increments∆d^s,i+1_n+1 and∆v^f,i+1_n+1 . Hence, in compact matrix notation the following linear system is solved

Then, the guess of the nodal values for the next iteration stepi+ 1are obtained as



The iterative procedure is aborted when a convergence criterion is met, e.g.

where^sand^f are problem specific tolerances. In the following Section 4.3.4.1 comments con-cerning the solution of linear system are made. Subsequently, some special cases are presented, where the linear system is manipulated before solving. In Section 4.3.4.2 a method for the ap-plication of impermeability constraints is proposed. There, the constraints are enforce with the Lagrange multiplier method. The resulting additional degrees of freedom are condensed out of the linear system. In Section 4.3.4.3 non-matching fluid and skeleton meshes are treated. The nodal values are projected from one mesh to the other. This leads to manipulation of the cou-pling blocks of the linear system.

4.3.4.1. Linear Solver

The main drawback of the monolithic solution approach compared to partitioned schemes is the increased size of the linear system. As both fields are solved simultaneously, the linear sys-tem includes the tangential matrices from both residuals as well as the respective linearizations of the coupling terms. This leads to a complex matrix pattern, as exemplarily depicted in Fig-ure 4.4. For large problems the most expensive part of the numerical simulation is the linear

Figure 4.4: Exemplary sparsity pattern of tangent matrix of monolithic porous medium system.

Black dot denotes non-zero entry.

solver. In this work, an iterative solution technique known as Generalized Minimal RESidual (GMRES) [211] is applied. In addition, for efficient calculations a block-preconditioner is used.

It has been specifically designed for monolithic coupled problems, see [240].

4.3.4.2. Impermeability Constraints

In the following a methodology to apply no-penetration conditions for Darcy flow and no-slip conditions for Darcy-Brinkman flow is presented. This is realized mainly by manipulating the linear system within the Newton scheme. After consistent linearization and before application of those mentioned conditions, the system to be solved in every iteration stepi+ 1of the non-linear

algorithm reads

for time step n + 1. Therein, the degrees of freedom are split into inner degrees of freedom (marked with the subscript I) and conditioned degrees of freedom (marked with subscript C).

Note that in this notation all fluid pressure degrees of freedom are counted as inner degrees of freedom, i.e. included in ∆v^f_I. The individual sub-matrices of the system matrix in equa-tion (4.145) denote the derivatives of the residuals with respect to the degrees of freedom, split into inner and conditioned degrees of freedom, i.e. for instance

S_CI = ∂R^s_C

In case of impermeable boundaries this system is extended by the Lagrange multiplier field and the constraint itself. The most important steps in order to condense them out of the system are sketched in the following.

4.3.4.2.1. No-Penetration Condition for Darcy Equation

In case of Darcy flow, the no-penetration constraints (4.102), (4.103) and (4.106) are de-manded in a weak sense. This leads to the discrete residuals R^f_C⊥ = R^f_C⊥(d^s,v^f), R^f_Ck = R^f_Ck(d^s,λ)andR^f,λ_C =R^f,λ_C (d^s,λ)of the normal, tangential condition and the additional bound-ary integrals, respectively. Here λ denotes the vector of nodal unknown Lagrange multipliers.

Including

R^f_C⊥ = 0, (4.150)

R^f_Ck = 0 (4.151)

into the monolithic system (4.145) and adding the boundary integral residual to the fluid equation (forth line in (4.145))

R^f_C+R^f,λ_C = 0, (4.152)

leads after consistent linearization to



n. Thus, the system size is increased by the number of nodal Lagrange multipliers. A noteworthy technical detail, is that in contrast to the other unknown quantities, the system is solved for the discrete Lagrange multiplier λⁱ⁺¹_n+1, instead of an increment. Also, the time discretized tangential condition is demanded at the time stepn+ 1.

The additional sub matrices denote the derivatives of the following residuals with respect to the degrees of freedom, i.e. for instance

N^FF_CC = ∂R^f_C⊥

Solving for the Lagrange multiplier in the Darcy flow equation (forth line in (4.153)) gives λⁱ⁺¹_n+1=−1

θD^−T·(−R^f_C−F_CI·∆v^f_I−F_CC·∆v^f_C−C^FS_CI ·∆d^s_I (4.159)

−C^FS_CC·∆d^s_C+ (1−θ)D^T_n ·λn), (4.160)

where all indices referring to the actual time stepn+1on the right hand side are omitted for clar-ity of notation. The Lagrange multipliers are discretized by dual shape functions. This idea was first presented in [255]. Details on the construction approach used here can be found in [193].

Utilizing these shape functions leads to a diagonal matrixD^T and therefore enables straightfor-ward and efficient evaluation ofD^−T, and hence of equation (4.160). By insertingλⁱ⁺¹_n+1into the tangential condition (sixth line in (4.153)), the Lagrange multiplier can be condensed out of the system and the original system size is restored.

The final system to be solved by the non-linear algorithm then reads

 with A = ¹_θT^λ_C ·D^−T. After each iteration, the Lagrange multiplier needs to be recovered by evaluating equation (4.160), as it is required for the calculation ofT^λ_C.

Remark 4.7 It has to be noted, that the construction of dual shape functions for NURBS is topic of current research [44, 225] and not considered here. Fortunately, the higher continuity of NURBS is not necessary for Darcy flow, as elaborated in Section 4.3.1.

4.3.4.2.2. No-Slip Condition for Darcy-Brinkman Equation

The methodology to include no-slip constraints into the monolithic solution scheme is very similar to modifications performed in FSI problems, see [148, 168]. After including the no-slip condition (4.104), its consistent linearization, and after adding the additional boundary integrals, the system reads

withRf^s_C = R^s_C −(1−θ) D^T·λ

n As done before,λⁱ⁺¹_n+1 can be condensed out of the system by solving for the Lagrange multiplier λⁱ⁺¹_n+1 in the Darcy-Brinkman flow equation (forth line in (4.162)). Additionally, the last equation from (4.162), i.e. the identity ∆v^f_C = _θ∆t¹ ∆d^s_C, is used to condense the fluid velocities on the constraint boundary, further decreasing the number of degrees of freedom. Note, that this simple relation only holds for conforming structure and fluid meshes. In case of non-conforming meshes a projection operatorPis needed to relate velocities and displacements ∆v^f_C = _θ∆t¹ P∆d^s_C. Such an operator can, for instance, be calculated from a mortar approach, as done for FSI with non-matching interfaces [148]. The final condensed system then reads

It has to be pointed out, that in this case the Lagrange multiplier does not influence the single non-linear iteration steps directly. It rather plays the role of a boundary traction which can be post-processed after solving one time step. This is because the momentum of the whole mixture already includes the coupling with the fluid stresses, which are equal to the Lagrange multiplier on the constraint boundary and therefore cancels out as described in Section 4.3.2.3.

4.3.4.3. Porous Medium Problem on Non-Matching Meshes

From a numerical point of view an interesting aspect is the consideration of non-matching fluid and solid meshes. If, for instance, the flow is complex, but the deformations are easily re-solved, it would be computationally advantageous to only use a very small number of elements to resolve the displacement field while the fluid velocity is approximated with far more elements.

Independent meshing of the fluid and the solid also implies the possibility of using different interpolation functions, even if the element sizes do not differ significantly. A general method-ology for solving volumetrically coupled multi-field problems is presented in [91]. Therein, two coupling methods are presented: one is based on straightforward interpolation of nodal values on the respective mesh and the other uses a mortar-like coupling of the primary fields. The former will be called collocation approach, while the latter will be referred to as mortar-based approach in the following. In the end, both methods provide the global algorithm with two discrete projec-tion operatorsP₁₂and P₂₁, which are used to transfer nodal information from one mesh to the other and vice versa. In this thesis, no further details on the construction of these operators are given. Instead, it is referred to the above contribution. It will be shown here, how this projection is included into the monolithic porous medium system.

The two different discrete domainsΩ^h₁ andΩ^h₂ are considered. Both domains approximate the same continuous porous volume. After spatial discretization, the solid problem is going to be solved on mesh Ω^h₁ and the fluid problem on Ω^h₂. The two discrete sets of primary unknowns are the skeleton unknowns d^s₁ on Ω^h₁ and the fluid unknowns v^f₂ on Ω^h₂. The auxiliary coun-terparts of the primary unknowns are d^s₂ on Ω^h₂ and v^f₁ on Ω^h₁, respectively. Simply speaking,

both approaches demand the equality of the primary and the corresponding auxiliary variable in a certain sense. The collocation method enforces the coupling in a strong, point-wise sense at each node. The mortar-based approach demands an equality of the two fields in a weak, in-tegral sense. Thereby, so-called dual shape functions [255] are used for efficient calculations.

Discretizatuion with standard finite elements then gives the desired projection of the discrete nodel valuesd^s₂ =P21d^s₁. Again, it is referred to [91] for details.

Now, the form of the linearized monolithic system will be sketched. The discrete residuals read

R^s₁(d^s₁,v^f₁) =0 and R^f₂(d^s₂,v^f₂) =0. (4.164) As before, the arising system of non-linear discrete algebraic equations is going to be solved in a fully monolithic manner. As iterative non-linear solution technique, a Newton-Raphson method is employed. This requires linearization of the residuals with respect to the primary unknowns, which is obtained from a truncated Taylor expansion:

LinR^s₁(d^s,i₁ ,v^f,i₁ ) = R^s₁(d^s,i₁ ,v^f,i₁ ) + ∂R^s₁(d^s₁,v^f₁) In order to eliminate the auxiliary quantities ∆d^s,i+1₂ and ∆v^f,i+1₁ from the resulting system of equations, the nodal transfer operatorsP12andP21are employed:

∆d^s,i+1₂ =P₂₁∆d^s,i+1₁ and ∆v^f,i+1₁ =P₁₂∆v^f,i+1₂ . (4.167) Hence,P₁₂is a projection of nodal values from meshΩ2 to meshΩ1 andP₂₁ a projection from meshΩ1 to meshΩ2. Written in matrix form, the off-diagonal blocks of the resulting system of equations are multiplied with the projection operators:



Equations (4.165) and (4.166) have to be fulfilled at each Newton step. After convergence of the Newton scheme, the incremental solutions∆d^s,i+1₁ and∆v^f,i+1₂ are used to perform the updates

d^s,i+1₁ =d^s,i₁ + ∆d^s,i+1₁ , v^f,i+1₂ =v^f,i₂ + ∆v^f,i+1₂ , (4.169) and the auxiliary quantities are recovered from the primary unknowns employing equation (4.167).

The procedure is repeated until a user-defined convergence criterion is fulfilled, as described in Section 4.3.4.

Im Dokument A Computational Approach to Coupled Poroelastic Media Problems  (Seite 85-93)