Downwind and upwind approximations for mesh and model adaptivity of elasto‐plastic composites

The use of heterogeneous materials, such as composites with Prandtl‐Reuss‐type material laws, has increased in industrial praxis, making finite element modeling with homogenization techniques a well‐accepted tool. These methods are particularly advantageous to account for microstructural mechanisms which can be related to nonlinearities and time‐dependency due to elasto‐plasticity behavior. However, their advantages are diminished by increasing computational demand. The present contribution deals with the balance of accuracy and numerical efficiency of nonlinear homogenization associated with a framework of goal‐oriented adaptivity, which takes into account error accumulation over time. To this end, model adaptivity of homogenization methods is coupled to mesh adaptivity on the macro scale. Our new proposed adaptive procedure is driven by a goal‐oriented a posteriori error estimator based on duality techniques using downwind and upwind approximations. Due to nonlinearities and time‐dependency of the plasticity, the estimation of error transport and error generation is obtained with a backward‐in‐time dual method despite a high demand on memory capacity. In this contribution, the dual problem is solved with a forward‐in‐time dual method that allows estimating the full error during the resolution of the primal problem without the need for extra memory capacity. Finally, a numerical example illustrates the effectiveness of the proposed adaptive approach.


INTRODUCTION
A fully coupled adaptive strategy is established to simultaneously control model and spatial discretization errors of the finite element method (FEM) by simulating nonlinear homogenization processes in accordance with a specific quantity of interest and duality techniques [1][2][3][4][5][6].The goal of this contribution is to develop a framework of goal-oriented adaptivity, taking into account errors due to space discretization and elasto-plastic mean-field and full-field homogenization models over time using downwind and upwind approximations.Exemplarily, downwind and upwind approximations have been used in the literature to solve initial value problems [7].In this contribution, the application of upwind and downwind approximations leads to the resolution of the primal and dual problems with forwards-in-time and backwards-in-time methods.
F I G U R E 1 Illustration of a two-scale problem.

A TWO-SCALE PROBLEM
A two-scale mechanical problem of first order is considered by assuming a scale separation, that is, ∕ ≪ 1, as illustrated in Figure 1, where the macro and micro domains are characterized by the symbols Ω and Ω, respectively.To each material point of the macro domain Ω ⊂ ℝ   , a micro domain Ω ⊂ ℝ   is associated, where the subscripts 0 and  denote the reference and the current configuration, respectively.A framework of small strains is assumed.The micro domain Ω is represented by a representative volume element (RVE) and a two-dimensional problem is considered (  = 2).We let  = ( 0 , ] be the total time interval of interest, where  0 is the initial time and  is the total time.The micro displacement vector  =  −  where  and  are the position vectors with respect to a reference point , respectively in the reference and the current state. The underlying microscopic equilibrium problem  and the macroscopic initial value boundary problem (IVBP) P read The microscopic equilibrium problem: In ,  denotes the microscopic stress tensor; Equation (1.3) describes the constitutive relation in rate form, where the effective material tangent ℙ() relates the stress rate σ() and the strain rate ε().Equation (1.4) prescribes some proper boundary conditions like periodic boundary conditions [2].The scale transition between the micro problem  and the macro problem P is established by where Equation (2c) is the Hill-Mandel condition [8] and the volume averaging operator on a domain Ω is defined as the space integral ⟨•⟩ = where  is referred to as the micro strain localization tensor and  ∈ Ω,  ∈ .
Mean-field formulation for composites: The mean-field method resolves the relevant microscopic fields to their means over  r individual material phases distinguished by a subscript  = 0, 1, ⋯,  r .The volume fraction of each phase is denoted by  r , such that ∑  r r=0 = 1.The subscript 0 is assigned to the matrix material, while the subscript  = 1, ⋯,  r denotes the r-th inclusion, such that the localization rule Equation (3.2) becomes In Equation (5.1), ̇ is the average strain over the local domain Ω r occupied by the r-th material phase.The average strain localization tensor  r in Equation (5.3) is obtained by inserting Equation (3.1) into Equation (5.2).The mean-field version of effective properties Equation ( 6) is obtained from Equation (4) as Here, ℙ r represents the Prandtl-Reuss material tensor of the phase  obtained in combination with Equation (5.4).

The macroscopic IVBP:
In the equilibrium condition (1.6), the body force is symbolized by b.
We also assume a continuous functional dependency such that In this way, the macroscopic stress tensor σ is dependent on the displacement vector ū( X) at the position X with initial value ū0 .Consequently, Equation (7.2) becomes the constitutive relation in rate form Equation (1.8), where the effective material tangent P( ū) relates the stress rate σ( ū) and the strain rate ̇( ū) of a Prandtl-Reuss type material law.Equation (1.7) prescribes the small strain tensor ε as a symmetric gradient of the displacement vector ū.As indicated in Equation (1.5), the boundary conditions of the microstructure have in common a prescription in terms of the macroscopic strain ε( ū).In order to solve the macro problems P with the FEM, we introduce a weak-form equivalent to the formulations for the macro problem P in Equation ( 1): Additionally, use has been made of the general formulation for the inner product in < •, • >∶= ∫ Ω < •, • > dV, now the comma (,) replaced by (⋅) for F( ū) and replaced by (∶) for B( ū;  ū).Generally, F(⋅) is a linear form, while B(⋅; ⋅) is a bilinear form.Observe from Equation (10) that the stress rate bilinear form B ( ū;  ū) in Equation ( 11) is defined.Furthermore, the strain rate bilinear form B ( ū;  ū) can also be defined from Equation (10) as illustrated in ref. [9].

Downwind and upwind temporal approximations:
The variational local form in Equation ( 8) can be used for discretization in space and in time.A general scheme is obtained by so-called cG()dG() methods representing a space-time discretization with continuous piecewise polynomials of degree  in space and discontinuous piecewise polynomials of degree  in time.We begin with a partition of the time interval  = ( 0 , ] as  0 <  1 < ⋯ <   < ⋯ <   =  and use the notations Consequently, the bilinear form in Equation (10) becomes B ( ū,  ū) = ∑  =1 B  ( ū,  ū)+ < ū( 0 ) ∶  ū( 0 ) >.We select the semi-discrete spaces for test and trial functions and we apply the discontinuous approximation ū ≈ ūr of the exact solution ū in Equation ( 8).For the semi-discrete spaces r  associated to each time interval   ,  r (  , V) denotes the space of polynomials of degree less than or equal to r ∈ ℕ 0 taking values in V.Note that  ūr ( 0 ) has to be specified separately in Equation ( 13) since  0 ∉  1 .In the sequel, we give a brief description of the procedure to drive a discretized time stepping scheme.Exemplarily, we consider the stress rate bilinear form in Equation (11).Its integration by parts renders the integral over the subinterval Note that Equation ( 16) gives an exact solution of the stress rate bilinear form in Equation (11).Next, Equation ( 16) is approximated using the jump of the displacement ū at each nodal point   defined as [ ū]  = ū+  − ū−  , where ū+(−)  = lim ↓(↑)0 ū(  + ) [7].Therefore, a downwind approximation and an upwind approximation are applied on the stress tensor σ( ūr ) and the test strain tensor  ε( ūr ) in Equation ( 16).Then, the right-hand side of Equation ( 16) is integrated by parts again and we end up with four stress rate bilinear form approximates.

GOAL-ORIENTED ADAPTIVITY ON THE MACRO-SCALE
A quantity of interest, which could be any differentiable function of a user's interest, is written as the following general form containing both a time integral part with regards to  1 and further part  2 at the final time .Due to discretization and model errors, a goal-oriented error measure for the quantity of interest in Equation ( 17) is defined as and its secant forms is written as Note that the secant forms in Equation ( 19) and Equation (20) are linear in the error ēr  ∶= ū − ūr which is the error of the primal solution [3].These secant forms allow us to define the dual problem D ∶ Find z ∈ , such that   ( ūr ; z,  ū) ∶=   ( ū, ūr ;  ū) + ρ ( ū, ūr ; z,  ū) = 0 ∀  ū ∈  .
( , zr  =∶   = const, zr −1 =∶  −1 , ūr  =∶   = const,  ūr  =∶   = const for  = 1,  − 1 with regards to the fact that the residual   ( ūr , zr ;  ūr ) can be localized, and can be solved by a time-stepping scheme when a dependency of jump terms [( z)] is taken into account.However, the residual   ( ūr , zr ;  ūr ) cannot be localized for each time increment   when a dependency of jump terms [( ūr  )] is considered.Therefore, Equation ( 22) can lead to the backwards-in-time method for dual solutions due to the terms related to the  + 1st time step as follows Last, we consider an error representation in terms of the projection  z ∈ r ⊂  and which satisfies the global discontinuous-in-time Petrov-Galerkin applied to the primal problem ρ( ūr ;  z) = 0 in Equation ( 8).Consequently, Equation (28.2) may be replaced by ( ū; ūr ) = ρ( ūr ; z) − ρ( ūr ;  z) = ρ( ūr ; z −  z).Note that an accurate estimate of the exact global error  j is introduced as viewed in ref. [4].

NUMERICAL EXAMPLES
We investigated a compact tension (CT) specimen and assumed a plane strain state.Due to symmetry properties, a half model is considered, as illustrated in Figure 2A.In Figure 2B, the initial mesh discretization is shown.The specimen is stretched by a displacement ū * = 0.08 mm in the vertical direction, which is uniformly distributed on the entire boundary of the hole.The quantity of interest is defined as a local type quantity represents a component of the macro stress tensor σ in the main direction 2. Ω ′ is the local green domain in Figure 2A, with  = 4 mm.The CT-specimen is assumed to be inhomogeneous on the micro-scale as its microstructure is composed of composite materials which consist of a matrix and an inclusion.Note that the inclusion is randomly deviated from the center position of the matrix, as illustrated in Figure 3 both, mean-field methods are less accurate due to the fact that they do not consider the location of inclusion inside the matrix, as depicted in Figure 3F.The latter is more accurate as the exact location of inclusion inside the matrix is taken into account (Figures 3A-3E).For the inelastic micro constituent, the classical von Mises plasticity is assumed.A total of 100 time steps are used during each mesh and model adaptivity step.The reference solution of quantity of interest  ref is obtained by using a uniformly refined mesh from the last adaptive mesh with a uniform full-field model distribution n = 2.
Figures 4A-4B and 4C-4D illustrate that the local mesh refinements are initially concentrated within the local domain Ω ′ before spreading to nearby areas for the forwards-in-time method and the backwards-in-time method, respectively.
However, Figures 4E-4F and in 4G-4H show that full-field method are mostly used within the local domain Ω ′ before extending to neighboring elements for both methods.Figure 5 displays mean-field and full-field model distributions by increasing the number of adaptive step  for the forwards-in-time method and the backwards-in-time method.The actual error Êrel decreases from 27.11% to 1.79% and from 31.11% to 2.51% for the forwards-in-time method and the backwardsin-time method, respectively.Note that an actual error Êrel smaller than 7.5% is already reached with 13 adaptive steps for both methods.Figures 6A and 6B illustrates the quantity of interest  j , the relative actual error Êj rel , and the relative error estimate Ẽj rel over each adaptive step . Figure 6A presents a convergence of the quantity of interest  j to the reference value  ref while Figure 6B shows that the relative actual error Êj rel is already reduced significantly and effectively after 13 adaptive steps.Additionally, Figure 6C illustrates that an actual error lower then 10% can be obtained by the use of 30% of elements with the FEM.Furthermore, with regards to Figure 6, the results obtained with the forwards-in-time method are quite similar to the result obtained with the backwards-in-time method.

CONCLUSION
In summary, several downwind and upwind approximations for primal and dual problems of elasto-plasticity have been presented and applied to a goal-oriented framework considering mesh and model error, where a forwards-in-time method has been developed to solve dual problem without additional storage requirements compared to backwards-in-time method.We then compared the results of both methods: The results show that the simulation of homogenization techniques combined with an adaptive approach is quite effective compared to an uniform full-field homogenization method due to the computational cost.

A C K N O W L E D G M E N T S
The investigations are financially supported by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) under grant number MA 1979/30-2.
Open access funding enabled and organized by Projekt DEAL.
. The matrix is an aluminum alloy with elasto-plastic properties characterized by nonlinear isotropic hardening and linear kinematic hardening ( = 76500[MPa],  = 0.3,  0 = 200[MPa],  1 = 300[MPa],  = 111,  = 608.5, 2 = 300[MPa],  = 1), while all inclusion fibers are ceramic with linear elastic properties ( = 400000[MPa] and  = 0.2).The volume fraction of inclusion  is equal to 0.1.The model hierarchy is based on the basic mean-field model self-consistent (SC), the mean-field interaction direct derivative (IDD) and the full-field FEM under periodic boundary condition with hierarchical order (n = 0), (n = 1), and (n > 1), respectively.The former, F I G U R E 3 RVE for Full-Field method (A-E), RVE for Mean-Field method (F).RVE, representative volume element.F I G U R E 4 CT specimen under tension with mesh and model adaptivity: adaptive mesh refinement distribution with the forwards-in-time (A-B) and backwards-in-time (C-D) methods; adaptive refined model distribution with the forwards-in-time (E-F) and backwards-in-time (G-H) methods.CT, compact tension.

F
I G U R E 5 CT specimen under tension with mesh and model adaptivity.Model distributions for different error levels (mean-field in green and full-field in red): forwards-in-time (A-D), backwards-in-time (E-H).CT, compact tension.F I G U R E 6 CT specimen under tension with mesh and model adaptivity: (A) quantities of interest versus adaptive step; (B) relative actual error Êrel and relative error estimate Ẽrel versus adaptive step; (C) relative actual error Êrel and relative error estimate Ẽrel versus the fraction  fem of elements using the FEM.FEM, finite element method.
Equation (3.2) is only valid if  is independent on time .The identity in Equation (3.3) is obtained by inserting Equation (3.1) into Equation (2a) where   denotes the fourth-order symmetric identity tensor.