Latest revision as of 13:39, 6 July 2021

Abstract

This article is focused on the study of a micro-macro LaTIn based Domain Decomposition Method (LaTIn-DDM) for the prediction of the nonlinear behavior of slender composite structures subjected to bending, buckling and delamination. Previous studies have shown that an adequate selection of the iterative parameters (search directions and macroscopic space) allow to improve the convergence rate and ensure scalability (i.e. number of iterations is independent of the number of subdomains) of the iterative schema. To obtain precise solutions, only the size reduction of the subdomains' discretization has been addressed (${\textstyle h}$-refinement), disregarding the option of increasing the polynomial degree of the finite elements (${\textstyle p}$-refinement) and ignoring their underlying effects on the information's transmission through the interfaces between subdomains. In this work and using linear and quadratic finite elements, ${\textstyle h}$ and ${\textstyle p}$ refinements on the subdomains and local ${\textstyle h}$-refinement only along the edges of the subdomains were investigated. It is conclude that the ${\textstyle p}$-refinement in the whole subdomain not only enables to reach more exact solutions than using global or local ${\textstyle h}$-refinement, but also the convergence rate is improved. These enhancements allow more complex simulations but using less degrees of freedom and less calculation time, even up to 97% faster.

Keywords: Domain decomposition method, quadratic finite elements, composites, delamination, buckling

1. Introduction

Since the middle of the last century, composite materials have been widely used in several industrial applications, showing advantages over materials such as steel and aluminum due to their specific properties. Furthermore, scientists and engineers have made efforts to understand their behavior and to predict them [1,2]. The usage of models which are defined at micro-length scales would be ideal, but numerical complexity and computational limitations (memory and time) appear when simulations are performed [3]. Instead, Domain Decomposition Methods (DDM) [4,5,6] are suitable to face these issues taking advantage of the power of parallel and distributed computations (high performance computing). By partitioning the structure into smaller subdomains connected by interfaces, these algorithms lead with local problems defined in each subdomain and a condensed interface problem. Then computational limitations are overcome because they are numerically cheaper and adapted to the parallel architecture of hardwares. The scalability of these methods (i.e. convergence rate does not depends on the number of subdomains) is often managed using a coarse problem ensuring a global communication between subdomains (i.e. a second calculation's scale) [7,8,9].

The LaTIn method [10] is in principle a non-incremental schema to solve non-linear problems; however its extension as a multiscale mixed DDM has been easily done [11]. Contact [12], crack propagation [13,14,15], buckling-delamination interaction [16,17] and heterogeneous structures [11] are among the different nonlinear and complex problems solved with this method. It distinguishes the no-linear equations defined on the subdomains from the non-linear ones defined on the interfaces, in order to define two groups of partial solutions over which a fixed point is applied to reach the intersection of them at convergence. This algorithm is configured with two search directions linking the force to the displacement distributions on the interfaces (i.e. mixed or Robin conditions). The LaTIn method aims to improve its convergence rate by introducing a second space scale at the interface level (classically called the “macroscopic problem” in contrast to the local problems solved in each subdomain which are called “microscopic problems”). For this reason, the global equilibrium over the structure is enforced in a few interface's unknowns by defining a macroscopic basis of the interface's displacements. Additionally, the Robin conditions need to be optimized (large-wavelength components converge rapidly whereas small-wavelength components converge slowly) as shown in [18,19,20].

Delamination is one of the main degradation mechanisms of laminated composite materials. This phenomenon is generally initiated by large interlaminar stresses and can be accelerated under geometric instabilities as buckling, leading eventually to a structural failure. For simulating the buckling-delamination interaction in composite laminates, the work of [18] uses the mesoscopic scale defining two constituents: the plies (3D elements) and the interfaces (2D elements), which are naturally linked to the domain partitioning. The geometrically nonlinear evolution is handled through a total Lagrangian formulation as proposed by [16], while delamination is modeled using a Cohesive Zone Model (CZM) which are based on damage mechanics [21]. Conditions of unilateral contact are considered to avoid interpenetration between delaminated surfaces.

The simulations previously carried out need a large amount of degrees of freedom (dof) to correctly capture the different local and global phenomena. Until now, the strategy has privileged to use sufficiently refined meshes (${\textstyle h}$-refinement), but this has implyied expensive computations. The other classical technique to reach more exact solutions rather than dividing elements into smaller ones is to increase the polynomial degree of the finite element approximation (${\textstyle p}$-refinement), as proposed in [22,23] for problems using direct solvers (without parallel computations). Therefore, this work studies the influence of the ${\textstyle p}$-refinement not only on the accuracy of the results, but also on the iterative solver (LaTIn-DDM). The numerical implementation is made using the C++ research code MULTI developed at the Laboratoire de Mécanique et Technologies de Cachan¹ using MPI and METIS libraries for the parallel assignments.

This work proposes to use second order finite elements in the LaTIn-DDM and to study their effects on the convergence rate and on the calculation time with respect to the LaTIn-DDM previously defined in [18]. To achieve this goal, Section 2 shows the general aspects of the multiscale strategy, then the reference problem, the domain partitioning, governing equations and the multiscale iterative algorithm are detailed. Subsequently, in Section 3, different ${\textstyle h}$- and ${\textstyle p}$-refinements are compared in 3D beam problems: bending, buckling and delamination. Finally, the conclusions and ongoing work are presented in Section 4.

(¹) LMT-Cachan (ENS Paris-Saclay/CNRS/Université Paris-Saclay)

2. The micro-macro LaTIn based Domain Decomposition Method

Let us consider a laminate composite (Figure 1) occupying the domain ${\textstyle \Omega }$ bounded by ${\textstyle \partial \Omega }$ in the current configuration, and consisting of ${\textstyle N_{P}}$ plies. Each ply ${\textstyle P}$ is connected to an adjacent ply ${\textstyle P'}$ through the interface ${\textstyle {\Gamma _{PP'}}}$. The structure is subjected to an external surface traction field ${\textstyle {{\underline {F}}_{d}}}$ on the part ${\textstyle {\partial \Omega _{F_{d}}}}$ and to a displacement field ${\textstyle {{\underline {U}}_{d}}}$ on the complementary part ${\textstyle {\partial \Omega _{U_{d}}}}$. The body force per unit of mass is written ${\textstyle {{\underline {f}}_{d}}}$. The relevant quantities are described in reference to the undeformed configuration using the index ${\textstyle \cdot _{\mathit {0}}}$. The geometrically nonlinear evolution is handled through a total Lagrangian formulation and delamination (damageable interfaces) is modeled using CZM and unilateral contact conditions. For the sake of simplicity, an extensive description of the CZM is found in [24], while [25] describe contact inequalities.

To propose the partitioned problem, the whole domain ${\textstyle \Omega }$ is split into subdomains which are connected by interfaces with mechanical behaviors. Two possibilities are considered: “material” interfaces between plies with localized non-linearities (damage, contact) that are compatible with the mesomodeling, and “numerical” interfaces (the perfect ones) within the plies to conceive smaller problems that are suited for parallelism, as schematized in Figure 1. A subdomain ${\textstyle {S_{\mathit {0}}^{\quad }}}$ defined in the undeformed domain ${\textstyle {\Omega _{S_{\mathit {0}}}}}$ is connected to an adjacent undeformed subdomain ${\textstyle {{S'}_{\mathit {0}}}}$ through an undeformed interface ${\textstyle {\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}=\partial {\Omega _{S_{\mathit {0}}}}\cap \partial {\Omega _{S_{\mathit {0}}'}}}$. The surface entity ${\textstyle {\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}}$ applies force distributions ${\textstyle {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}}$, ${\textstyle {{\underline {F}}_{S_{\mathit {0}}'}}}$ and displacement distributions ${\textstyle {{\underline {W}}_{S_{\mathit {0}}^{\quad }}}}$, ${\textstyle {{\underline {W}}_{{S'}_{\mathit {0}}}}}$ to ${\textstyle {S_{\mathit {0}}^{\quad }}}$ and ${\textstyle {{S'}_{\mathit {0}}}}$ respectively (Figure 2). Let us define ${\textstyle {\Gamma _{S_{\mathit {0}}}}=\cup _{{{S'}_{\mathit {0}}}\in {\mathbf {E} }}{\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}}$. For a subdomain ${\textstyle {S_{\mathit {0}}^{\quad }}}$ such that ${\textstyle {\Gamma _{S_{\mathit {0}}}}\cap ({\partial \Omega _{F_{d_{\mathit {0}}}}}\cup {\partial \Omega _{U_{d_{\mathit {0}}}}})\neq \emptyset }$, the boundary condition ${\textstyle ({{\underline {F}}_{d_{\mathit {0}}}},{{\underline {U}}_{d_{\mathit {0}}}})}$ is applied through a boundary interface ${\textstyle {\Gamma _{{S}_{d_{\mathit {0}}}}}}$.

Figure 1. The reference problem, the mesomodel and its partitioning [18]

Figure 2. Subdomains and interfaces [16]

The purpose of the method is to find the subdomain fields ${\textstyle {{\underline {u}}_{S_{\mathit {0}}^{\quad }}}}$ (displacement) and ${\textstyle {{\underline {\underline {\pi }}}_{S_{\mathit {0}}^{\quad }}}}$ (second Piola-Kirchhoff stress), and the interface fields ${\textstyle {{\underline {W}}_{S_{\mathit {0}}^{\quad }}}}$ (displacement) and ${\textstyle {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}}$ (inter-forces), where the index ${\textstyle \cdot _{S_{\mathit {0}}}}$ ranges over all subdomains. After partitioning, the governing equations are separated into two groups:

1. Non-linear equations in subdomains and macroscopic admissibility of interfaces, whose solutions are elements of the manifold ${\textstyle \mathbf {A_{d}} }$:

• non-linear kinematic admissibility of the subdomains:

• non-linear static admissibility of the subdomains:

where ${\textstyle {{\underline {\underline {\dot {E}}}}({\underline {u}}_{S_{\mathit {0}}^{\quad }}^{\star })}={\frac {1}{2}}({\nabla }{{\underline {u}}_{S_{\mathit {0}}}^{\star }}+{}^{t}{\nabla }{{\underline {u}}_{S_{\mathit {0}}}^{\star }}+{}^{t}{\nabla }{{\underline {u}}_{S_{\mathit {0}}^{\quad }}}{\nabla }{{\underline {u}}_{S_{\mathit {0}}}^{\star }}+{}^{t}{\nabla }{{\underline {u}}_{S_{\mathit {0}}}^{\star }}{\nabla }{{\underline {u}}_{S_{\mathit {0}}^{\quad }}})}$.

• behavior of the subdomains:

where ${\textstyle \psi }$ is the stored energy function per unit of undeformed volume. For this work, ${\textstyle \psi ={\frac {1}{2}}\mathbf {K} _{S_{\mathit {0}}}{{\underline {\underline {E}}}_{S_{\mathit {0}}}}:{{\underline {\underline {E}}}_{S_{\mathit {0}}}}}$ has been used.

• macroscopic admissibility of the interfaces (after the linearization of the previous equations), which is a partial verification of the action-reaction principle:

where the subspace ${\textstyle {{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{M}}}$ of interface macroscopic admissible displacements is a parameter described in Section 2.1.

2. Local (non-linear) equations in the interfaces whose solutions belong to the manifold ${\textstyle {\mathbf {L} _{\Gamma }}}$:

• constitutive equation ${\textstyle \displaystyle {\mathcal {R}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}([{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}]\,,\,{{\underline {F}}_{S_{\mathit {0}}^{\quad }}}\,,\,{{\underline {F}}_{S_{\mathit {0}}'}})={\underline {0}}\;{\hbox{over}}\;{\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}\in {\Gamma _{S_{\mathit {0}}}}\;,}$ where ${\textstyle [{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}]={{\underline {W}}_{{S'}_{\mathit {0}}}}-{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}}$ is the interface displacement gap.
• boundary behavior ${\textstyle {\mathcal {R}}_{S_{d_{\mathit {0}}}}({{\underline {W}}_{S_{\mathit {0}}^{\quad }}},{{\underline {F}}_{S_{\mathit {0}}^{\quad }}})={\underline {0}}\;{\hbox{over the boundary}}\;{\Gamma _{{S}_{d_{\mathit {0}}}}}\;}$. The relation ${\textstyle {\mathcal {R}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}={\underline {0}}}$ can be made explicit for:

- Perfect interface: ${\textstyle \left\{{\begin{array}{l}{{\underline {F}}_{S_{\mathit {0}}^{\quad }}}+{{\underline {F}}_{S_{\mathit {0}}'}}={\underline {0}}\\[][{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}]={\underline {0}}\end{array}}\right.}$

- Cohesive interface: ${\textstyle \left\{{\begin{array}{l}{{\underline {F}}_{S_{\mathit {0}}^{\quad }}}+{{\underline {F}}_{S_{\mathit {0}}'}}={\underline {0}}\\\displaystyle {\mathcal {A}}_{PP'}([{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}]\,,\,{{\underline {F}}_{S_{\mathit {0}}^{\quad }}})={\underline {0}}\end{array}}\right.}$

- Unilateral contact interface: ${\textstyle \left\{{\begin{array}{l}{{\underline {F}}_{S_{\mathit {0}}^{\quad }}}+{{\underline {F}}_{S_{\mathit {0}}'}}={\underline {0}}\\{\underline {n}}\cdot [{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}]\geq 0\;\;\mathbf {or} \;\;{\underline {n}}\cdot {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}\geq 0\\({\underline {n}}\cdot [{{\underline {W}}_{S_{\mathit {0}}^{\quad }}}])({\underline {n}}\cdot {{\underline {F}}_{S_{\mathit {0}}^{\quad }}})=0\\\mathbf {P} {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}=\mathbf {P} {{\underline {F}}_{S_{\mathit {0}}'}}={\underline {0}}\\\end{array}}\right.}$

where ${\textstyle {\underline {n}}}$ is the unit normal to ${\textstyle \Gamma _{SS'}}$ and ${\textstyle \mathbf {P} }$ is the corresponding tangential projection operator.

The algorithm consists in seeking the interface solution ${\textstyle s_{ref}}$ alternatively in these two spaces: first, one finds a solution ${\textstyle s=(s_{S_{\mathit {0}}^{\quad }})_{{S_{\mathit {0}}^{\quad }}\in {\mathbf {E} }}=({{\underline {W}}_{S_{\mathit {0}}^{\quad }}},{{\underline {F}}_{S_{\mathit {0}}^{\quad }}})_{{S_{\mathit {0}}^{\quad }}\in {\mathbf {E} }}}$ in ${\textstyle \mathbf {A_{d}} }$, then a solution ${\textstyle {\widehat {s}}=({\widehat {s}}_{S_{\mathit {0}}^{\quad }})_{{S_{\mathit {0}}^{\quad }}\in {\mathbf {E} }}=({{\underline {\widehat {W}}}_{S_{\mathit {0}}^{\quad }}},{{\underline {\widehat {F}}}_{S_{\mathit {0}}^{\quad }}})_{{S_{\mathit {0}}^{\quad }}\in {\mathbf {E} }}}$ in ${\textstyle {\mathbf {L} _{\Gamma }}}$. In order for the two problems to be well-posed, two search directions ${\textstyle {k^{+}}}$ and ${\textstyle {k^{-}}}$, which link the solutions ${\textstyle s}$ and ${\textstyle {\widehat {s}}}$ during the iterative process, are introduced. The converged interface solution ${\textstyle s_{ref}}$ is such that ${\textstyle s_{ref}\in \mathbf {A_{d}} \cap {\mathbf {L} _{\Gamma }}}$. A complete description of this iterative method can be found in [16].

In order to evaluate the convergence, a criterion based on the verification of the constitutive laws of the interfaces quantities issued from the admissibility stage has been implemented, as proposed by [26].

2.1 Separation of the micro-macro scales

The definition of the macroscopic quantities is done through an orthogonal projector ${\textstyle \Pi }$ for each interface ${\textstyle {\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}}$, such as the macroscopic fields are ${\textstyle {{\underline {F}}_{S_{\mathit {0}}^{\quad }}^{M}}~=~\Pi {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}}$ and ${\textstyle {{\underline {W}}_{S_{\mathit {0}}^{\quad }}^{M}}~=~\Pi {{\underline {W}}_{S_{\mathit {0}}^{\quad }}}}$. Consequently, the microscopic spaces, ${\textstyle {{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{m}}}$ and ${\textstyle {{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{m}}}$, are orthogonal to the macroscopic spaces, ${\textstyle {{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{M}}}$ and ${\textstyle {{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{M}}}$, with respect to the inner product ${\textstyle L^{2}({\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}})}$, so the scales can be uncoupled with respect to the interface virtual work as follows:

In order to ensure the scalability of the iterative scheme, a global linear coarse grid problem 5 has been introduced. It is fully characterized by the set of interface macroscopic spaces ${\textstyle ({{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{M}})}$ which is a parameter of the method as well as its dual ${\textstyle ({{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}}}^{M}})}$. Usually, a common macroscopic basis for both the traction and displacement macroscopic fields is chosen so that the uniqueness of the micro-macro descomposition is ensured. Classically, the macroscopic space contain at least the affine part of the interface displacements; this corresponds to verify the balance of the first moments of forces at the interface. However, if complex structures are involved, the geometry and its partitioning degrade the convergence rate. To palliate this problem, [18] propose to find an optimized search direction instead to enrich the macroscopic space.

In this paper, simulations have been carried out using the standard linear macroscopic space and the local search direction previously used in [16].

2.2 Discretization

To solve both equations' sets of the multiscale algorithm, interfaces and subdomains are discretized in space using classical finite elements. At an interface ${\textstyle {\Gamma _{S_{\mathit {0}}{S'}_{\mathit {0}}}}}$, the displacements ${\textstyle {{\underline {W}}_{S_{\mathit {0}}^{\quad }}}}$ and forces ${\textstyle {{\underline {F}}_{S_{\mathit {0}}^{\quad }}}}$ belong to the approximation spaces ${\textstyle {{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}},h}}}$ and ${\textstyle {{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}},h}}}$. These spaces are chosen such that the bilinear form (in the sense of the interface mechanical work) is non-degenerate. Additionally, a wrong discretization for ${\textstyle {{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}},h}}}$ could generate spurious oscillating modes leading to numerical instability. These inconveniences can be evaded using a common space for the displacements and forces and including a local refinement of the mesh near the boundary (over-discretization) of the subdomains (approximation space ${\textstyle {{\mathcal {U}}_{S_{\mathit {0}}}}}$). Two manners are possible: to increase the number of elements (${\textstyle h}$-version) or to use a higher degree of approximation (${\textstyle p}$-version) for the field ${\textstyle {{\underline {u}}_{S_{\mathit {0}}^{\quad }}}}$ near to the interface, as illustrated in Figure 3 [27]. Classically, the code MULTI has considered linear elements ${\textstyle {\mathcal {P}}_{1}}$ with local ${\textstyle h}$-refinement along the subdomain's boundary, while the spaces ${\textstyle {{\mathcal {W}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}},h}}}$ and ${\textstyle {{\mathcal {F}}_{{S_{\mathit {0}}^{\quad }}{{S'}_{\mathit {0}}},h}}}$ are piecewise constant functions ${\textstyle {\mathcal {P}}_{0}}$.

In this work, the ${\textstyle p}$-refinement is also explored. Indeed, it is here proposed to use the second order approximation for ${\textstyle {{\underline {u}}_{S_{\mathit {0}}^{\quad }}}}$ not only in the boundary but in the whole subdomain.

Figure 3. Modification of the classical approximations of the inter-force and local displacement along the edge of a subdomain: ${\displaystyle h}$-and ${\displaystyle p}$-versions [13]

3. Influence of the discretization on the LaTIn-DDM

This section analyses the influence of the subdomains' discretization for three different problems: bending, buckling and delamination. The following three meshes ${\textstyle {{\mathcal {U}}_{S_{\hbox{0}}}}}$ are considered:

• ${\textstyle i}$-version: linear six-node wedge elements of equal size in the whole subdomain (“initial” version without local over discretization, (Figure 3);
• ${\textstyle h}$-version: linear six-node wedge elements, where the elements along the subdomain's boundary are divided into three four-node tetrahedral elements as shown in Figure 4 (local ${\textstyle h}$ over discretization);
• ${\textstyle p}$-version: quadratic fifteen-node wedge elements of equal size in the whole subdomain (global ${\textstyle p}$ over discretization).

For all cases, the approximation functions of the interfaces are ${\textstyle {\mathcal {P}}_{0}}$. Displacements, convergence rate and time are used to compare the refinements.

Figure 4. before and after local ${\displaystyle h}$-refinement. (a) Wedge element. (b) Mesh of the subdomain's boundary (the inner mesh remains unchanged)

3.1 Bending

The problem consists of a thin 3D beam whose length is ${\textstyle L=160}$ [mm] (${\textstyle x}$-direction), width is ${\textstyle a=120}$ [mm] (${\textstyle y}$-direction) and thickness is ${\textstyle h=4}$ [mm] (${\textstyle z}$-direction). Regarding the boundary conditions, the cantilever beam is embedded at one end and subjected to a surface load ${\textstyle {\underline {F}}_{d_{0}}=0.1}$ [N/m${\textstyle ^{2}}$] on its upper face. Here, the hypothesis of small displacements is considered. First an isotropic material is studied, then an orthotropic one is analysed to study the influence of material heterogeneities.

3.1.1 Isotropic material

It is considered an elastic modulus ${\textstyle E=210}$ [GPa] and a poisson coefficient ${\textstyle \nu =0.3}$ [-]. The geometry is partitioned into ${\textstyle 64}$ identical subdomains and ${\textstyle 184}$ interfaces (Figure 5); each subdomain has length ${\textstyle L_{sst}=20}$ [mm], width ${\textstyle a_{sst}=15}$ [mm] and thickness ${\textstyle h_{sst}=4}$ [mm]. Five meshes are considered: two different initial discretization with their respective ${\textstyle h}$-refinement while the last one is a ${\textstyle p}$-version. Each subdomain has ${\textstyle n_{x}}$(${\textstyle n_{y}}$,${\textstyle n_{z}}$) elements in the ${\textstyle x}$(${\textstyle y}$,${\textstyle z}$)-direction, respectively, as shown in Table 1 as well as the total number of elements and the total degrees of freedom used for each mesh. In order to estimate the solution's error (Table 1), the maximum vertical displacement is compared with the theoretical elastic curve [28].

Figure 5. Partitioning of the geometry (bending problem with isotropic material)

In Figure 6a, theoretical and numerical displacements of the neutral line are compared. As naturally expected, when increasing the number of linear elements in the thickness, the solution's error decreases, but the dof and the computational cost increase. However, it is important to notice that for the ${\textstyle h}$-refinements (a.2) and (b.2), the displacement's error and calculation time increase (Table 1) while the convergence rate decrease respect to the corresponding ${\textstyle i}$-versions (a.1) and (b.1). Even after 1000 iterations, the iterative LaTIn error for mesh (a.2) is twice the detention criteria. This phenomena could be explained by the fact that ${\textstyle h}$-version has only a localized refinement along the edge of a subdomain, while the element's size inside the subdomain remains the same as the ${\textstyle i}$-version. This choice could induce different stiffness through a subdomain, implying additional difficulties to transfer informationn between subdomains.

Finally, the mesh (d) (${\textstyle p}$-version) is twice more accurate, has 73% less dof and is 19,3% more quickly than mesh (b.1). Differences in the computation time (Table 1) are mainly related to the mesh size, because convergence rates are similar, as shown in Figure 6b.

Figure 6. Bending problem with isotropic material. (a) Deflection. (b) Evolution of the iterative LaTIn error

3.1.2 Orthotropic material

The precedent problem is now studied considering a composite laminate made of four plies ${\textstyle [0^{o},90^{o}]_{S}}$, each 1 [mm] a thick. A ${\textstyle 0^{o}}$-layer is transversely isotropic with the following elastic properties: ${\textstyle E_{1}=165}$ [GPa], ${\textstyle E_{2}=E_{3}=9}$ [GPa], ${\textstyle \nu _{12}=\nu _{13}=0.3}$ [-], ${\textstyle \nu _{23}=0.5}$ [-], ${\textstyle G_{12}=G_{13}=5.6}$ [GPa] and ${\textstyle G_{23}=2.8}$ [GPa]. The geometry is divided into 256 identical subdomains of ${\textstyle L_{sst}=20}$ [mm], ${\textstyle a_{sst}=15}$ [mm] and ${\textstyle h_{sst}=1}$ [mm], generating 736 interfaces. Therefore, for each ply there is one subdomain in the ${\textstyle z}$-direction (four in total through the thickness). Table 2 compares the different discretizations.

In Figure 7a is observed that the vertical displacement of the neutral axis is similar except for simulations (a.1) and (a.2). In these two cases, the iterative LaTIn error (Figure 7b) is twice the detention criteria, even after 1000 iterations. Using the ${\textstyle p}$-version mesh (c), the LaTIn error is less than ${\textstyle 10^{-5}}$ in only 68 iteraciones, less than half of the iterations made by the curve (b.1) to converge. If the upper stresses ${\textstyle \sigma _{xx}}$ are compared respect to the theoretical ones (Figure 8), it is possible to confirm that (a.1), (a.2) and (b.2) do not fit the desired solution.

Figure 7. Bending problem with orthotropic material. (a) Deflection. (b) Evolution of the iterative LaTIn error

Table 2 compares calculation time, total number of elements and total dof. It is verified that using ${\textstyle p}$-version, mesh (c), it is possible to obtain the same results as in (b.1) but consuming only 54.2% of the time, even considering that the number of dof increases in a 50%. This is explained because (c) has the best convergence rate (Figure 7b).

Figure 8. Results for the orthotropic material: normal stresses ${\displaystyle \sigma _{xx}}$

Similar to the precedent example, the local ${\textstyle h}$-version does not ensure better results. Therefore, the meshes considered for the next examples will be ${\textstyle i}$ and ${\textstyle p}$-versions.

3.2 Buckling

The problem to be addressed is a slender 3D beam built-in at both ends, with one of them subjected to a negative displacement ${\textstyle {\underline {u}}_{d}}$ to produce uniaxial compression, while a perpendicular perturbation ${\textstyle {\underline {F}}_{d}}$ induces buckling out of the plane (Figure 9a). The structure has the following geometry: length ${\textstyle L=10}$ [mm] (${\textstyle x}$-direction), width ${\textstyle a=1}$ [mm] (${\textstyle y}$-direction) and thickness ${\textstyle h=0.1}$ [mm] (${\textstyle z}$-direction). The properties of the material are ${\textstyle E=135}$ [GPa] and ${\textstyle \nu =0.3}$ [-]. The geometry is divided into 100 subdomains and 156 perfect interfaces, therefore, each subdomain has ${\textstyle L_{sst}=0.2}$ [mm], width ${\textstyle a_{sst}=0.5}$ [mm] and thickness ${\textstyle h_{sst}=0.1}$ [mm]. Three meshes are considered, the first two are linear without over discretization (${\textstyle i}$-version), while the mesh (c) is a ${\textstyle p}$-refinement. More details of the meshes and their results are shown in Table 3.

The simulation was performed in 1000 steps. Figure 9a shows the initial configuration and the final deformation at the last time step. The evolution of the compression axial load (${\textstyle P/P_{crit}}$) is obtained as a function of the transverse displacement in ${\textstyle x=L/2}$, as shown in Figure 9b.

It is noticed that simulations (a) and (b) has respectively ${\textstyle 8.68\%}$ and ${\textstyle 2.61\%}$ of error in the critical load when the transverse displacement over ${\textstyle L_{0}}$ is 0.005 [-], while mesh (c) is the closest (only ${\textstyle 1.46\%}$ of error at the same point). In addition, the time spent for mesh (c) is only the ${\textstyle 2.9\%}$ used in (b).

Figure 9. (a) The initial configuration and the final deformation after the last time step. (b) The load-displacement curve

3.3 Delamination

In this section we study the effect of discretization when problems involve CZM. The example to be simulated is a 3D double cantilever beam (DCB), whose length is ${\textstyle L=20}$ [mm] (${\textstyle x}$-direction), width ${\textstyle a=2}$ [mm] (${\textstyle y}$-direction), thickness ${\textstyle h=1}$ [mm] (${\textstyle z}$-direction) and pre-crack ${\textstyle a_{0}=10}$ [mm] located at the end of the beam along the ${\textstyle x}$-direction (Figure 11a). The properties of the material are ${\textstyle E=135}$ [GPa] and ${\textstyle \nu =0.3}$ [-]; the cohesive interface parameters are ${\textstyle k_{n}=100\cdot 10^{3}}$ [N/mm${\textstyle ^{3}}$], ${\textstyle \alpha =1}$ [-], ${\textstyle Y_{c}=0.4}$ [N/mm] and ${\textstyle n=0.5}$ [-].

The geometry is divided into 160 subdomains and 324 interfaces such as each subdomain has ${\textstyle L_{sst}=0.5}$ [mm], width ${\textstyle a_{sst}=1}$ [mm] and thickness ${\textstyle h_{sst}=0.5}$ [mm]. Four meshes were considered (see details in Table 4) and the simulation was performed in 50 time steps.