6. Global-Local Approach for Non-conforming Finite Element Discretiza-
6.2. Benchmark numerical problems
6.3.2. Example 2: I-shaped specimen at finite strain
To gain deeper insights into the performance of the Global-Local formulation equipped with non-conforming discretization techniques, a tensile test on the I-shaped specimen undergoing large deformation is performed. Here, the dual mortar method is used, which its computational efficiency has been demonstrated in the previous example. Following [6, 97], a BVP is applied to the specimen, which is depicted in Fig. 6.14. The main challenge here is to demonstrate the efficiency of the Global-Local formulation from a computational standpoint. In particular, (i) we aim to fix the large global domain in terms of the geometry and modeling, whiledifferent local domains are introduced at the region of interest. Thus, for the simulation standpoint, the large structure at BG (which typically comes from industry) is kept fixed, and only the geometry ofBLwill be changed;
(ii) Furthermore, we aim to use non-conforming finite element discretization, and hybrid finite element discretization. Specifically, this means that the local discretization element’s type isT1 (triangular isoparametric element) while the global mesh is set with Q1 (i.e., quadrilateral isoparametric element); see Fig. 6.15.
Here, for the global domain, we set l1 = 25 mm, l2 = 18.8 mm, l3 = 22.4 mm,
H1 = 110mm, and w1 = 22 mm with a quarter circular notched shape with a radius of r1 = 3.625 mm. The discretizations of the global domain contains 532 Q1 elements. For the local BVP in BL, we consider three different case studies:
• Case 1. The local domain includes double notched at the left and right edges, as shown in Fig. 6.14b. Here, we set H2 = 22.4 mm, w2 = 14.8 mm, l4 = 13 mm, l5 = 3.7 mm, and r2 = 2.5 mm. The discretizations of the local domain contain 16,484 T1 elements.
• Case 2. The local domain includes voids to weaken the specimen and facilitate the material cracking. Here, dimensions of the domain are identical to Case 1. Three types of voids are considered: small, medium, and big voids, which are distributed randomly with radius r= (0.94,1.58,3.15)mm, respectively. The discretizations of the local domain contain 12,155 T1 elements.
• Case 3. The local domain includes hard inclusions to stiffen the specimen, and thus more time is required to initiate the cracking. Here, dimensions of the domain are identical to Case 1. Three types of hard inclusion are considered: small, medium, and big inclusions which are distributed randomly with radiusr= (1.28,1.98,2.8)mm, respectively. The local domain partition contains 18,074 T1 elements.
Here, the elastic material property is set for the hyperelastic Neo-Hookean model identical to the numerical example given in Section 5.8.2. For Case 3, the mismatch ratio
a b
B𝐺
𝑦 𝑥 𝜕𝐷B
𝜕𝐷B 𝒖 = (0, ത𝑢)
Γ1 Γ2
Case 1 Case 2 Case 3
B𝐿 B𝐿 B𝐿
𝑤1 𝑙1
𝑙2 𝑙3 𝐻1
𝑟1
𝑤2 𝑤2 𝑤2
𝑟2 𝐻2
𝑙4
𝑙5
Figure 6.14: Example 2. I-shaped fracture test at finite strain. (a) Specimen geometry and loading conditions for the global domainBG, and (b) the local domainsBLfor different case studies.
χ= EEf iber
matrix is set to χ = 10, which is a ratio between the inclusions and matrix Young’s modules.
The overall response of the Global-Local approach in terms of accuracy/robustness and efficiency for the two nested finite element models is verified using a reference problem.
The reference domain is made from simultaneously considering BG, including BL within a single geometry. In the reference domain, a single type of finite element is used by the Q1 quadrilateral isoparametric element.
The load-displacement curves that correspond to Cases 1, 2, and 3 are depicted in Figs. 6.16a, 6.17a, and 6.18a, respectively. Accordingly, a comparison of the reference solution and the Global-Local formulationsg/l−1 andg/l−2 are presented. The resulting Global-Local curve for bothg/l−1 andg/l−2 in all examples are in very good agreement with the reference solution. Hence, the transition of local non-linearity and heterogeneity responses to the global level (through the up-scaling procedure) for all the numerical case studies is consistently performed. Notably, the fracturing process in Case 3 takes longer to initiate due to the hard inclusions (see Fig. 6.18a), while in Case 2, crack initiates earlier compared with the others because of voids being used in the BL to weaken the material stiffness; see Fig. 6.17a.
Next, we evaluate the corresponding accumulative computational time, that is the CPU simulation time per the prescribed displacement. The corresponding accumulative computational time for Cases 1, 2, and 3 are depicted in 6.16b, 6.17b and 6.18b, respec-tively. It can be grasped that the resulting accumulative computational times through g/l−2 is required 13%, 13.3% and 16% less time compared to the reference computational time for the Cases 1, 2, and 3, respectively. Additionally, the accumulative computational time through g/l−1 shows higher computational time compared with g/l−2. But, as
a b
Case 1 Case 2 Case 3
Figure 6.15: Example 2. Finite element discretization of (a) the global domain BG, and (b) the local domains BL for different case studies with the non-matching discrete interface, i.e., ΓG6= ΓL.
,KN ,sec
a b
, mm , mm
Figure 6.16: Example 2 (Case 1). Comparison of the reference solution and the Global-Local formulation g/l −1, and g/l −2. (a) Load-displacement curves, and (b) time-displacement curves in terms of ’accumulated time’.
,KN ,sec
a b
, mm , mm
Figure 6.17: Example 2 (Case 2). Comparison of the reference solution and the Global-Local formulation g/l −1, and g/l −2. (a) Load-displacement curves, and (b) time-displacement curves in terms of ’accumulated time’.
a b
, mm
,KN
, mm
,sec
Figure 6.18: Example 2 (Case 3). Comparison of the reference solution and the Global-Local formulation g/l −1, and g/l −2. (a) Load-displacement curves, and (b) time-displacement curves in terms of ’accumulated time’.
already mentioned, the advantage of using g/l−1 is being strongly non-intrusive for the coupling BG and BL and not computational efficiency (as observed for g/l−2).
Furthermore, we investigate the energy response when solving a reference problem
and a Global-Local approximation. Recall, the consistency of the energy functional (the departing point of the Global-Local approximation; see Section 5.2)
E(ϕ, s)≡Ee(ϕG,ϕL, sL,ϕΓ,λC,λL),
between referenceE and the Global-LocalEefunctional which are indicated in Formulation 3.2.2 and 5.2.1, respectively. We investigate this approximation by means of the evolution of the total stored elastic strain energy plotted in Fig. 6.19 for different case studies, and zooming into the framed region of the left plot. An important observation is that the g/−1 gives an overestimate energy response compared with the reference energy for all the numerical examples; see. Fig 6.19. This observation holds when Dirichlet boundary condition being used for the coupling of the global and local domains. This is similar to the computational homogenization when linear displacement boundary conditions are used; see [42]. In summary, the Global-Local simulation results for both g/l −1 and g/l−2 show very good agreement with the reference solution, for all case studies.
The computed crack phase-field is a good indicator for evaluating the down-scaling procedure (i.e., the transition of external loading increments from the global level to the local one). To this end, the crack phase-field profiles are presented in Fig. 6.20, Fig. 6.21, and Fig. 6.22 on the deformed configuration for Case 1, Case 2, and Case 3, respectively. A qualitative crack profile from the Global-Local formulation are compared with the reference solution. For visualization purposes, the elements wheresL<0.05 are removed from the contour plots. Notably, the deformed shapes are not magnified. Due to the different local structures used forBL, different crack patterns are observed. In Case 1, the crack initiates at the tip of the circular notched at the both left and right edges and continues to propagate straight till the end of the computation; see Fig. 6.20. In Case 2, the crack initiates from the void near the left circular notch, and after merging to that, it continues to propagate to the right edge till the end of the computation, see Fig. 6.21.
In the last case, the crack initiates between two fibers, as depicted in Fig. 6.22, and after merging to the left circular notch continues to propagate to the right edge until the end of the computation.
As a result, the crack phase-field profiles, for different time steps, are in very good agreement with the reference solution, which demonstrates the consistent transition be-tween the global and local BVPs. Thus, the feasibility/efficiency of using the proposed Global-Local approach equipped with non-matching discretization techniques are observed for these complicated local structures.
ഥ𝒖, mm
ഥ 𝒖, mm
ഥ 𝒖, mm ℰ−ℰ~ ℰ−ℰ~ ℰ−ℰ~
c b a
Figure 6.19: Example 2. Comparison of the total free-energy functional between the reference solution E and the Global-Local approach Eein (a) Case 1, (b) Case 2, and (c) Case 3 (left); zooming into the framed region of the left plot (right).
/ :
/ :
:
Figure 6.20: Example 2 (Case 1). Comparison of the crack phase-field profiles. (a) The reference solution s. The local solution sL through (b) g/l−1, and (c) g/l −2 on the deformed configuration at ¯u= [4.61,4.65,4.7,4.84] mm.
/ :
/ :
:
Figure 6.21: Example 2 (Case 2). Comparison of the crack phase-field profiles. (a) The reference solution s. The local solution sL through (b) g/l−1, and (c) g/l−2 on the deformed configuration at ¯u= [3.17,3.18,3.19,3.21] mm.
/ :
/ :
:
Figure 6.22: Example 2 (Case 3). Comparison of the crack phase-field profiles. (a) The reference solution s. The local solution sL through (b) g/l−1, and (c) g/l −2 on the deformed configuration at ¯u= [4.66,4.67,4.68,4.69] mm.
Conclusion and Future Research
7.1. Conclusion
Variational phase-field modeling is the regularized fractured formulation with a strong capability to simulate complicated failure processes. These include crack initiation (also in the absence of a crack tip singularity), propagation, coalescence, and branching without additional ad-hoc criteria. This feature is particularly attractive for industrial applica-tions, as it minimizes the need for time-consuming and expensive calibration tests. In contrast to these advantages, the finite element treatment of the phase-field formulation is known to be computationally demanding, mainly due to the non-convexity of the en-ergy functional to be minimized with respect to the displacement and the phase-field [46, 58, 128]. Other challenges for the phase-field fracturing formulation is two-fold.
First, that is a regularized-based formulation which is strongly linked to the element discretization sizehdue to the principal parameters: a small residual scalarκand charac-teristic length-scalel. Specifically,κ:=κ(h) andl :=l(h) hold such thathlandhκ through discretization error estimates [85]. Hence, the equations to be minimized for the variational phase-field formulation are strongly related to the element size h. Thus, for resolving the crack phase-field, a sufficiently small his chosen to obtain the experimental resolution.
The second challenge is to use the phase-field fracture approach for structures of indus-trial complexity. This has been the subject of limited investigations, and further studies in this direction will pave the way for the wide adoption of phase-field modeling within legacy codes for industrial applications.
In fact, when dealing with large structures, the failure behavior is solely analyzed in a (small) local region, whereas in the surrounding medium, a simplified and linearized system of equations can be solved. Thus, the idea of a two-scale formulation, in which the full displacement/phase-field problem is solved on a lower(local) scale while dealing with a purely elastic problem on an upper(global) level, is particularly appealing. These features
145
lead us to use the Global-Local approaches as they make it possible first to compute the global model elastically, and then to determine the critical areas to be re-analyzed, while storing the factorization of the structural stiffness decomposition. The local model is then iteratively substituted within the unchanged/fixed global one, which avoids the reconstruction of the global mesh. Here, we proposed an efficientvariational-based Global-Local approach for a phase-field formulation of a fracturing material at small and large deformations.
In this thesis, two types of Global-Local formulation were adopted for the phase-field brittle fracture. These Global-Local formulations are based on Dirichlet-Neumann-type boundary conditions and Robin-type boundary conditions, namely g/l−1 and g/l−2, respectively.
The first type, g/l−1, is strongly non-intrusive in the computational aspect enabling the formulation to be performed within legacy codes. Due to the extreme difference in stiffness between the global counterpart of the zone to be analyzed locally and its actual response when undergoing extensive cracking, relaxation/acceleration techniques are used.
These include the Aitken’s ∆2-method, the SR1, and Broyden’s methods. Our findings showed that the iterative convergence could be improved significantly, and to a similar extent for all investigated methods. We observed that Aitken’s ∆2-method is the most convenient choice for the implementation of the approach within legacy codes, as this method needs only tools already available for the so-called sub-modeling approach, which is well-known and widely used in industrial contexts.
The second type, g/l−2, has the advantage of computational efficiency. In contrast to the Dirichlet boundary conditions (being used ing/l−1), Robin-type boundary conditions (is used ing/l−2) did not lead to a stiff local response, particularly in a softening regime;
thus no ”extra efforts” were required (e.g., relaxation/acceleration procedures).
Here, the proposed framework for the Global-Local approach was first used for the material undergoing small deformation. A successful extension of the Global-Local ap-proach was further extended towards large defamation for bothg/l−1 and g/l−2. Here, the performance of the Global-Local formulation was demonstrated in a quantitatively and qualitatively manner, as follows:
• Crack patterns at the local scale at the complete failure state to evaluate the down-scaling procedure (i.e., the transition of external loading increments from the global level to the local one).
• Load-displacement curve for evaluating the up-scaling procedure during the Global-Local coupling approach (i.e., the transition of local non-linearity and heterogeneity responses to the global level).
• Investigations of the thermodynamical consistency between the reference energy and its Global-Local energy functional.
Accordingly, the overall response of the Global-Local approach in terms of accuracy and efficiency for the two nested finite element models was verified using a reference problem. The resulting crack phase-field patterns and load-displacement curves for both g/l−1 and g/l −2 showed very good agreement with the reference solution. This is
observed for all case studies in the entire range of loading, including the pre- and post-peak behavior. Furthermore, the consistency of the energy functional (the departing point of the Global-Local approximation) between the reference and the approximated Global-Local energy functional were investigated. An important observation was that the g/−1 gave an overestimated energy response compared with the reference energy due to the Dirichlet boundary condition, whereas g/l−2 showed very good agreement with the reference solution.
For deeper insight into the Global-Local formulation, we quantitatively investigated convergence performance based on the newly proposed indicator for the stopping criterion of the iterative process. Another evaluated factor is the corresponding accumulative com-putational time (i.e., CPU simulation time) per prescribed displacement. An important observation was that g/l−2 required fewer iterations for the convergence of the solution.
The required accumulative computational time through g/l−2 was approximately 14%
less than the reference computational time which underlines the efficiency of using Robin-type boundary conditions. Additionally, the accumulative computational time through g/l−1 shows the high computational time versusg/l−2. Yet, the advantage of using the g/l−1 is being strongly non-intrusive way of coupling (and not computational efficiency).
High efforts within g/l −1 are not surprising because the Global-Local problem has a larger discretization size than the reference problem and three nested iterative processes versus two for the reference problem.
Because the strong displacement continuity through the matching discretized inter-face is too restrictive (from the computational standpoint), the Global-Local formulation was further extended toward the non-conforming finite element discretization. The main advantage here is to achieve more regularity at the interface. Another impacting factor is to discretize the local domain independently from the global mesh, thereby gaining more freedom in terms of implementation. In this study, we investigated the mortar method, the dual mortar method, and the localized mortar method. Our findings showed that the mortar method and dual mortar method gave exact results, while the localized mor-tar method was not precise as other methods. On the other hand, the localized mormor-tar method was computationally faster, yet preserved computational accuracy with some or-der of approximation. The reason is that a Lagrange multiplier is assigned with a dual basis function in a dual mortar method, whereas in the localized mortar method is set with Dirac delta function. For all numerical examples, the excellent performance of the proposed framework in terms of accuracy for both g/l−1 and g/l−2 were observed.
7.2. Future research
Possible extensions of the Global-Local formulation for the fracturing material are as follows:
• In practical applications, when the evolving localization areas are not known ´a-priori, the Global-Local approach must be supplied with the possibility of the adaptive choice of the local domain. The adaptive procedure has two goals: (i) to adjust the local domain when fractures are propagating dynamically; and (ii) to reduce the total computational cost because the local domains are tailored to ´a-priori unknown fracture path.
• To approximate the Robin-type boundary conditions rather than using the exact representation. For instance, a two-scale approximation of the Schur complement within the Global-Local formulation is a possibility to reduce the computational effort, yet preserving accuracy; see [44].
• The careful investigation of the error evolution during an iterative Global-Local process to evaluate the decay of the iterative coupling error.
• A further application of the proposed framework which needs to be investigated is the phase-field modeling of ductile fracture undergoing small/large strains toward the efficient Global-Local approach.
• An extension of phase-field modeling within a multi-physics framework towards the Global-Local approach. For instance, hydraulic phase-field fracture, and thermo-elastic solids.
• Global-Local approaches easily allow different numerical techniques to be used for the global and local domains. In this regard, a flexible choice of the discretization scheme can be employed on each domain, individually, such as FEM [135], IGA [66], and VEM [3], among other numerical treatments.
• To extend the Global-Local formulation toward a three-dimensional setting. Thus, a rigorous numerical analysis must be left for future work.
Derivation of Robin-Type Boundary Conditions
In this section, we investigate the relationship between ∆ ˆu and ∆ ˆλ (in the incremental sense) for the complementary, fictitious and local domains at the converged solution state.
We aim to derive the Robin-type boundary conditions such that all coupling terms given in (C1), (C2), and (C3) are satisfied, simultaneously, at the Global-Local iteration k.
Recall, the complementary term used in (5.11), and letuC and λC be the stationary of the following functional:
L=L(uC,λC) :=
Z
BC
w(εC,1,1) dx+ Z
Γ
λC ·(uΓ−uC) ds− Z
ΓN,C
¯
τ ·uCds. (A.1) Here, Γ ∈ Rδ−1 ⊂ BC is denoted as an interface, and uΓ := tr uC ∈ H1/2(Γ) can be given implicitly, i.e., (Ce1)+(Ce2), or explicitly, i.e., (eC4). Recall that (A.1) lives in BC (the following description holds true for BF except τ¯ = 0). In the discrized setting, the stationary points of the energy functional for the L is characterized by the first-order necessary conditions through L1 = Lu(uC,λC;δu) =0 and L2 =Lλ(uC,λC;δλ) = 0.
Here, δu and δλ are test functions. We split L1 into inner nodes and interface nodes denoted as, {a, b}, respectively, by
La1(u) =fa−F¯ =! 0 x∈ B\Γ, Lb1(u)=fb−LTGλˆC =! 0 x∈Γ, L2 =LGuˆΓ−JGuˆC =! 0 x∈Γ.
(A.2)
Here, f= Z
BC
(BuG)Tσ(uC) dxis an internal nodal force vector and ¯F= Z
ΓN
(NuG)Tτ¯ds stands for the external force vector. It is trivial that the Lagrange multiplier acts as an
149
external force at the interface. A Newton-type solution for the residual-based system of equations for (uC,λC) is provided by the following linearization
(fa−F) +¯ Kaa∆ ˆuC,a+Kab∆ ˆuC,b =0,
(fb−LTGλˆC) +Kba∆ ˆuC,a+Kbb∆ ˆuC,b−LTG∆ ˆλC =0, (LGuˆΓ−LG∆ ˆub) +LG∆ ˆuΓ−LG∆ ˆuC,b =0,
(A.3)
whereK:=∂f/∂uˆC is the standard tangent stiffness matrix. Thus yielding the following iterative update
uˆC,a ←uˆC,a+ ∆ ˆuC,a, uˆC,b ←uˆC,b+ ∆ ˆuC,b, and ˆλC ←λˆC + ∆ ˆλC. (A.4) Let us now assume that the equilibrium state is achieved such thatLa1 =0,Lb1 =0, and L2 =0. Through further algebraic analysis, (A.3) is reduced to the interface which takes the following form
SC∆ ˆuC,b =SC∆ ˆuΓ =LTG∆ ˆλC with SC :=S(KC) =Kbb−KbbaK−1aaKab, (A.5) where S refers to the Steklov-Poincar´e mapping [70]. By means of (A.5), a displace-ment field ˆuC,b is extracted from the interface Γ and then through the Poincar´e-Steklov mapping S returns the outward normal stress derivative with respect to the trace of the displacement. This is called Dirichlet-to-Neumann mapping [33, 52]. Notably, because an identity uG =uC holds on Γ in the sense of a trace (see Section 4.2), then (A.5) can be written as,
SC∆ ˆuG,b =LTG∆ ˆλC. (A.6) Similarly, we have the following identities:
SL∆ ˆuL,b=TTL∆ ˆλL and SF∆ ˆuF,b =LTG∆ ˆλF. (A.7) Here, TL :=JL|ΓL is the restriction ofJL from BL to ΓL. Furthermore, we define SL:=
S(KL), and SF :=S(KF) in (A.7).
Proposition 1. Let the global solutions be at the converged state and let the following identity holds true:
uk,
1 2
Γ =ukΓ ∈ Γ, (A.8)
then, the Global-Local formulation is converged. In addition, (A.8) holds true if and only if
∆ΛL =ΛkL−Λk−1L = 0. (A.9)
Proof. Let we are at the converged solution of the g/l−2 at iteration k. The proof comprises two parts:
a. The Global-Local procedure is in the convergence stateif all the coupling terms, (C1),
a. The Global-Local procedure is in the convergence stateif all the coupling terms, (C1),