g/l − 1 via incremental setting augmented with relaxation/ acceleration

acceleration techniques

For later developments, it proves convenient to reformulate the global equation in incremental form. It is straightforward to see that for a given Global-Local iteration k ≥ 1, this reads: given the triple (u^k−1_G ,λ^k−1_F ,λ^k−1_C ) known from the iteration k−1, as well as (λ^k_F,λ^k_C) ’recovered’ at the iteration k, we solve

Z

B_G

σ(∆u_G) :ε(δu_G) dx− Z

Γ

(λ^k_F−λ^k−1_F )·δu_G da− Z

Γ

(λ^k_C−λ^k−1_C )·δu_Gda= 0, (G_incr) for ∆u_G =: (∆u_G)^k and set

u^k_G:=u^k−1_G + (∆uG)^k. (4.25) We term equation (4.25) a ’direct update’ within the Global-Local iterative procedure.

This is in contrast to the notion of a ’relaxed/accelerated update’ to be considered in the following sections.

Following [36] and [79], we will consider and incorporate two types of relaxation/

acceleration techniques into our approach: Aitken’s ∆²-method (also known as dynamic relaxation, whose efficient implementation in fluid-structure interaction computations has already been reported [75, 37]) and Newton correction. Within the family of Quasi-Newton correction formula, we restrict ourselves to the Symmetric Rank One (SR1) and the Broyden update versions. Technically, both types deal with the global solution update equation (4.25) and modify it specifically. Let us consider (4.25) written in terms of the nodal displacements

uˆ^k_G:= ˆu^k−1_G + (∆ ˆuG)^k s.t. (∆ ˆuG)^k =K⁻¹_G r^k_Γ, (4.26) where, owing to (G_incr), one has

KG :=

Z

B

(B_u^G)^TCeB_u^G dx, with Ce as a 3×3-matrix representation of C, and

r^k_Γ:=

Z

Γ

(N_u^G)^T(λ^k_F −λⁿ⁻¹_F ) da+ Z

Γ

(N_u^G)^T(λ^k_C −λⁿ⁻¹_C ) da.

4.6.1. Necessity of using relaxation/acceleration techniques

This section evaluates the convergence performance of the Global-Local iterative for-mulation. Let us now investigate the numerical example depicted in Fig. 4.5, to evaluate

the convergence performance of the Global-Local iterative scheme. The left edge of a ref-erence problem as well as the global domain, is fixed along the x−and y−directions, and uniform displacement (u_x = 2 mm) is applied monotonically in one increment over the right edge. A quarter of the global scale is considered to be the fictitious domain, B_F, for which re-localization needs to be performed. The local scale, B_L, by refining and keeping the outer boundary of the B_F, is designed by additionally including two soft inclusions.

The elastic properties are set asλ= 121.15 kN/mm², and the shear modulus isµ= 80.77 kN/mm² for both global and local levels.

a b c

Figure 4.5: Geometric setup, and boundary conditions. (a) Reference problem, (b) local domain, and (c) global domain.

By χ = _E^EL

G we mean the mismatch ratio between the local and global Young’s mod-ulus. The efficiency of the Global-Local coupling scheme for three different values of χ is investigated. As soon as χ becomes greater than one, the local response will be stiffer than the global level. Figure 4.6a presents the Global-Local iteration procedure based on different values of χ (represented in a logarithmic base).

a b

−

number of / iterations number of / iterations

−

Figure 4.6: Iterative convergence performance of the Global-Local coupling approach.

(a) Different χ= (1,4,10), and (b) at fixed χ= 4 by imposing differentω.

Here, we only plot the error indicatorηand not its ingredients for the iterative Global-Local formulation; see Section 4.5. We observe that the convergence performance of the Global-Local approach is affected by the stiffer response transferred from the local level.

Hence, if the local response is stiffer than the global one (e.g.,χ= 4 or 10), then divergence behavior is observed.

To overcome this issue, one could use an explicit damping factor (i.e., a constant parameter ω) to modify the incremental solution used in (4.26) (i.e., ω∆ ˆu_G). Now, let χ= 4 to investigate the effect of the different damping factors (and possibly improve the convergence behavior). Figure 4.6b compares the impact of the explicit damping factor on the convergence performance during the Global-Local iterations.

It is clear that by imposing ”an appropriate” damping factor, the stiff local response is relaxed, which for some values convergence is achieved (e.g.,ω = 0.4) while for some others divergence is obtained (e.g., ω = 1.5). Thus, an appropriate damping factor cannot be a priori known; therefore, specific treatment needs to be considered to identify a suitable damping factor for relaxing the Global-Local iterations; see Section 4.6.4. To this end, we discuss in detail relaxation/acceleration techniques in the following sections.

4.6.2. Symmetric Rank One method

The quasi-Newton method is used at the global level to obtain a new approximation of the global tangent stiffness matrix in the recursive manner. The main ingredient of this method is to replace the tangential stiffness matrix with its secant approximation. More precisely, the quasi-Newton method modifies (4.26) at any iteration k ≥ 2 by replacing the matrix KG with Ke^k =f(Ke^k−1,r^k_Γ,(∆ ˆuG)^k−1) s.t. Ke¹ =KG, thus resulting in

uˆ^k_G := ˆu^k−1_G + (Ke^k)⁻¹r^k_Γ . (4.27)

Let us now consider symmetric rank one (SR1) update. The classical recursive update formulation for the SR1 method renders the following form:

Ke_i+1=Ke_i+ rⁱ⁺¹_Γ (rⁱ⁺¹_Γ )^T

∆ ˆuⁱ_G(rⁱ⁺¹_Γ )^T . (4.28) Herein, the subscripted 0≤i < k designates the iteration number of the SR1 formulation that is recursively computed up to the current Global-Local iteration, k. Through SR1 formulation in (4.28) symmetric property of the tangent stiffness matrix still holds, which means that ifKe_i is symmetric then Ke_i+1 is also symmetric [32].

We apply the Sherman-Morrison formula [119]; that is, (A+∆A)⁻¹ =A⁻¹+A⁻¹B(A,

∆A)A⁻¹for any given square and invertible matrixAon (4.28). This results in the update formulation for the inverse of tangent stiffness matrix, which yields

Ke⁻¹_i+1=Ke⁻¹_i +Ke⁻¹_i B(Ke_i,∆Ke_i)Ke⁻¹_i , (4.29) such that,

B(Ke_i,∆Ke_i) =−rⁱ⁺¹_Γ (rⁱ⁺¹_Γ )^T

(rⁱ⁺¹_Γ )^T(∆ˆuⁱ_G+Ke⁻¹_i rⁱ⁺¹_Γ ) . (4.30) Multiplying the left- and right- hand sides of (4.29) by the current residual vector (i.e., r^k_Γ) results in

Ke⁻¹_i+1r^k_Γ=Ke⁻¹_i r^k_Γ−Ke⁻¹_i rⁱ⁺¹_Γ (rⁱ⁺¹_Γ )^T(Ke⁻¹_i r^k_Γ)

(rⁱ⁺¹_Γ )^T(∆ˆuⁱ_G+Ke⁻¹_i rⁱ⁺¹_Γ ) . (4.31) In the last iteration (i.e., i = k − 1), the modified SR1 recursive formulation (4.29) resembles ∆ˆu^k_G =Ke⁻¹_k r^k_Γ, which turns out to be the current global incremental solution.

For the implementation standpoint, we reformulate (4.31) in a compact form. Let us define w_i+1 =Ke⁻¹_i+1r^k_Γ, and v_i+1 = Ke⁻¹_i rⁱ⁺¹_Γ , for which the following re-formulation of

(4.31) is achieved:

wi+1 =wi+∆wi s.t. ∆wi =−vi+1

(rⁱ⁺¹_Γ )^Tw_i

(rⁱ⁺¹_Γ )^T(∆ˆuⁱ_G+v_i+1 . (4.32) Here, the recursive formulation is initialized through a given initial global solution, namely w₀ =u⁰_G. Furthermore, for the next Global-Local iteration, three vectorsv_k =w_k−1, r^k_Γ and ∆ˆu^k_G have to be retained. The corrected incremental global solution is obtained by the last recursive step (i.e., w_k=: ∆ˆu^k_G).

Remark 4.6.1. The SR1 algorithm is in convergence state if the updated term of Ke_i+1 (second term of the right-hand side) in (4.28) is well-defined [32]. Moreover, the denom-inator of the update term within (4.28) should not to be so small. Therefore, the update SR1 formulation given in (4.29) is used only if the following criterion is satisfied:

∆ˆuⁱ_G(rⁱ⁺¹_Γ )^T

2 > c₁

∆ˆuⁱ_G 2

rⁱ⁺¹_Γ 2.

Herein, c₁(0,1) is a constant variable that is usually taken to be square root of our floating-point computational environment,

√

10⁻¹⁶ [32]. If this criterion does not hold, then the SR1 update formulation is not applied (i.e., w_i+1 =w_i).

4.6.3. Broyden’s method

One of the drawbacks of the SR1 method in the context of the Global-Local cou-pling approach is storing three vectors at each Global-Local iteration. Additionally, the denominator of the SR1 update formulation can be very small quantity. Therefore, an alternative quasi-Newton method is also considered. Broyden’s method, as one of the standard quasi-Newton methods, can be used to accelerate the convergence rate of the Global-Local iterations. Analogous to the SR1 method, Broyden’s method is used at the global level for determining stable and efficient updates of the tangent stiffness matrix, but with less computational cost compared to SR1. The recursive update formulation for Broyden’s method through basic algebraic manipulation yields the following form:

Ke_i+1=Ke_i+rⁱ⁺¹_Γ (∆ˆuⁱ_G)^T

(∆ˆuⁱ_G)^T∆ˆuⁱ_G . (4.33) It can be observed that the update term of Ke_i+1 in the SR1 method, namely the second term of the right side of (4.28), is more influenced by the global residual vector, namely rⁱ_Γ, whereas the update term of Kei+1 in Broyden’s method in (4.33) is more sensitive to the global incremental displacement solution, namely ∆ˆuⁱ_G. Applying Sherman-Morrison [119] to (4.33) (in a similar manner to the SR1 method) and then multiplying left- and right-hand sides of equality by the current residual vector (i.e.,r^k_Γ) results in the following equation:

Ke⁻¹_i+1r^k_Γ=Ke⁻¹_i r^k_Γ−Ke⁻¹_i rⁱ⁺¹_Γ (∆ˆuⁱ_G)^T(Ke⁻¹_i r^k_Γ)

(∆ˆuⁱ_G)^T(∆ˆuⁱ_G+Ke⁻¹_i rⁱ⁺¹_Γ ) . (4.34) In contrast to the SR1 method, a significant feature of Broyden’s update formulation in (4.34), is its capability to approximate Ke⁻¹_i+1r^k_Γ with the knowledge of two quantities,

namely ∆ˆuⁱ_G and Ke⁻¹_i rⁱ⁺¹_Γ , through previous Global-Local iterations and not necessarily rⁱ⁺¹_Γ .

For the implementation standpoint, we need to reformulate (4.34) in a compact form.

Let us define w_i+1 = Ke⁻¹_i+1r^k_Γ, v_i+1 = Ke⁻¹_i rⁱ⁺¹_Γ , for which, in turn, the following re-formulation of (4.29) is achieved:

w_i+1 =w_i+∆w_i s.t. ∆w_i =−v_i+1 (∆ˆuⁱ_G)^Tw_i

(∆ˆuⁱ_G)^T(∆ˆuⁱ_G+v_i+1) . (4.35) Similar to the SR1 method, the recursive formulation is initialized through the given initial global solution, namely w₀ = u⁰_G . Moreover, for the next Global-Local iteration, the two vectors vk = wk−1 and ∆ˆu^k_G have to be retained. The corrected incremental global solution is obtained by the last recursive step; that isw_k =: ∆ˆu^k_G.

4.6.4. Aitken’s ∆² method

Aitken’s ∆² method [75, 100] is used to relax the interface displacement quantity by imposing a damping factor on the global incremental solution. This method modifies (4.26) at any iteration ofk ≥2 introducing the damping factorω_k−1 =f ω_k−2,(∆ ˆu_G)^k−1, (∆ ˆu_G)^k

s.t. ω₀ = 1 as follows:

uˆ^k_G := ˆu^k−1_G +ωk−1(∆ ˆuG)^k, (4.36) Considering a well-chosen damping factor, stability is enforced on the computational scheme which enhances the convergence speed of the Global-Local approach [48]. Within Aitken’s ∆² method, a damping factor has the following explicit form:

ωk−1 :=ωk−2

((∆ ˆu_G)^k−1)^T (∆ ˆu_G)^k−1−(∆ ˆu_G)^k

|(∆ ˆuG)^k−1−(∆ ˆuG)^k|² , k ≥2, (4.37) with ω₀ = 1. It is clear that in (4.37), two previous global information, namely ˆu^k_G and ˆ

u^k−1_G , are retained to establish ω and hence modify the current global solution.

Remark 4.6.2. It is known that by means of Aitken’s ∆² method for scalar sequence, the converged solution is achieved if ω ∈(0,2) [37]. In particular, the stabilization of the iterative solution method is achieved while ω <1, and ω >1 implies that the acceleration of the iterative process is achieved. But, for the vector sequence this kind of argument is not guaranteed.

Im Dokument Variational-based Global-Local approach for a phase-field formulation of a fracturing material (Seite 80-84)