A General, Implicit, Large-Strain FE 2 Framework for the Simulation of Dynamic Problems on Two Scales
LLarge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 1
A General, Implicit, Large-Strain FE Framework for theSimulation of Dynamic Problems on Two Scales
Erik Tamsen , (cid:63) , Daniel Balzani Institute of Mechanics and Shell Structures, Technische Universität Dresden, 01062 Dresden Chair of Continuum Mechanics, Ruhr University Bochum, 44801 Bochum (cid:63)
Email address of corresponding author: [email protected]
Abstract
In this paper we present a fully-coupled, two-scale homogenization method for dynamicloading in the spirit of FE methods. The framework considers the balance of linear mo-mentum including inertia at the microscale to capture possible dynamic effects arisingfrom micro heterogeneities. A finite-strain formulation is adapted to account for geo-metrical nonlinearities enabling the study of e.g. plasticity or fiber pullout, which maybe associated with large deformations. A consistent kinematic scale link is establishedas displacement constraint on the whole representative volume element. The consistentmacroscopic material tangent moduli are derived including micro inertia in closed form.These can easily be calculated with a loop over all microscopic finite elements, only ap-plying existing assembly and solving procedures. Thus, making it suitable for standardfinite element program architectures. Numerical examples of a layered periodic materialare presented and compared to direct numerical simulations to demonstrate the capabilityof the proposed framework. Keywords:
Computational Homogenization, RVE, Microscopic Inertia, Volume Constraint,Consistent Tangent Modulus
1. Introduction
Micro-heterogeneous materials can give rise to wave propagation under dynamic loading, lead-ing to distinct stress distributions at the microscale. Thus, resulting in a complex macroscopicmaterial behavior. There are various examples from different fields of application which cover abroad range of length scales. Currently the interest is high in metamaterials in general, but espe-cially in locally resonant metamaterials exhibiting special properties like band gaps and negativebulk moduli. The applications range from cloaking devices [8, 21] over tunable sound attenua-tion [25, 30] to earthquake protection [5, 35]. More classical materials are being investigated aswell. One example is metaconcrete, which replaces aggregates by rubber-coated lead inclusionsto weaken impact waves [23, 36]. A different approach for an improved impact resistance arestrain-hardening cement-based composites (SHCC), that show a pronounced energy dissipationunder dynamic loading, as well as a change in fiber failure and overall crack pattern [9, 10, 11]. Inaddition to that, porous materials have shown an influence of microscopic inertia on voids underhigh strain rates [37, 44]. This list illustrates the possible influence of the material microstructureon the macroscopic response under high dynamic loading for a wide range of materials and ap-plications. In general, any material under impact loading, which has large variation in stiffness,e.g. rubber-coated particles, or materials with a pronounced variation in density, e.g. if pores orcracks are present at the microscale, can exhibit distinct effective macroscopic properties resultingfrom the dynamics in the microstructure. a r X i v : . [ c s . C E ] O c t . Tamsen, D. Balzani 2 To model the before-mentioned effects, the dynamics at the microscale need to be accounted for.Computational homogenization methods for quasi-static loading have become a common tool innumerical material analysis, see e.g. [14, 15, 33, 38, 46, 49]. A recent overview of computationalhomogenization methods in general is given in [17]. In addition a recapitulation of the FE method in particular is presented in [45]. However, with the rise in interest in metamaterials moreand more dynamic homogenization frameworks have been published in the last years. One is themethod of asymptotic expansion, see e.g., [7, 16, 19, 20], which is mainly based on the original workof Bensoussan et. al. [1]. Then there is the more general theory of elastodynamic homogenizationby Willis [51], which has been applied in e.g. [34, 40, 52] and others. The two approaches arebased on related ideas and their similarities are studied in [39]. Both methods are limited toelastic, periodic media. A more general approach is the micro-macro simulation based on arepresentative volume element (RVE). In such homogenization methods, macroscopic quantitiessuch as the deformation gradient and the displacement vector at a macroscopic integration pointare projected onto a microscopic boundary value problem whose homogenized response replacesthe constitutive material law. In the case where the finite element method (FEM) is used onboth scales, this procedure is called the FE method. A good theoretic introduction to thistheory including dynamics is given by [13], which has been the foundation on which this paper isbuilt. The framework in [16] calculates a quasi-static microstructure but then applies an inertia-induced eigenstrain based on the microstructure as an extra body force at the macroscale toaccount for micro-inertia effects. This was extended by [22] to account for matrix cracking atthe microscale under impact loading. Other, rather FE -type schemes as [28, 29, 42, 43, 47, 50]calculate the full balance of linear momentum at the microscale. In [28] an explicit, periodic,small-strain framework is presented, which was extended to an implicit time integration methodfor modeling resonant elastic metamaterials in [29]. [42, 43] use the assumption of linear elasticityto improve the computational performance, by splitting the problem into a purely static and aspecial dynamic boundary value problem (BVP). To better capture a wider range of appliedfrequencies, [47] use a Floquet-Bloch transformation to build a base of eigenmodes to analyzeelastic, periodic metamaterials. The mentioned frameworks all use at least one of the followingidealizations: small strains, linear elasticity, periodic or symmetric microstructures. In addition,many require quite elaborate implementations.The aim of this paper is to build a multiscale framework for dynamic loading as general aspossible, while still being compatible with standard FE architecture. To enable the analysisof micro-mechanical processes as plasticity or fiber pullout, as well as to incorporate effectsresulting from geometric nonlinearities, the proposed framework uses a finite-strain formulation.Its importance is supported by the analytical example in [24] of a nonlinear, elastic metamaterial,where finite strains were shown to be relevant for large wave amplitudes. To permit a flexibledamage evolution in the RVE, which is not dominated by periodic boundary conditions [6], wepropose to apply kinematic scale links as constraints on the whole RVE. This allows us to modelany type of RVE morphology.The paper starts by discussing the general ideas of the framework in Section 2. The used no-tations of the large-strain framework are presented and the averaging relations derived in [12]are briefly recapitulated to include the full balance of momentum at the microscale. Then inSection 3 the FE formulations of the microscale are presented and the kinematic constraints arederived. In Section 4 the respective macroscale formulations are displayed. There, a focus islaid on the derivation of consistent macroscale tangent moduli in closed form. These enable aquadratically converging macroscopic Newton-iteration, resulting in a robust and efficient algo-rithm compared to numerical differentiation through perturbation of the macroscopic quantities[32]. Once the whole theoretical background of the framework is explained, Section 5 providesnumerical examples, demonstrating the convergence behavior and analyzing some properties ofRVEs under dynamic loading. The paper is completed in Section 6 with concluding remarks. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 3
2. Homogenization Framework Including Microscopic Inertia
The general idea of the homogenization framework for dynamics is to consider the full balance oflinear momentum including the inertia terms at the microscale. The enables not only the anal-ysis of full dynamic fields at the microscale but a direct study of microscopic inertia effects onthe macroscale. By using appropriate averaging relations and kinematic links, a consistent scalebridging for dynamic loading is established. In the FE method, each macroscopic Gauss pointis associated with a separate microscopic RVE simulation which uses the macroscopic mechan-ical quantities to define the microscopic BVP. In order to differentiate the two scales, variablesassociated with the macroscale are denoted with a bar • . Here, the finite-strain framework istaken into account in order to enable the analysis of a wide range of material behavior and micro-mechanical effects under dynamic loading. In the following sections the fundamental ingredientsof the scale-coupling framework are explained, see also the schematic illustration in Figure 1. The connection of the (undeformed) reference and the (deformed) current configuration is de-scribed by the displacement u as u = x − X , where X ∈ B denotes the coordinates in thereference configuration and x ∈ S the deformation. The link between the two configurations interms of transformations of vector elements is described by the deformation and displacementgradients, respectively F = ∂ X x = + H and H = ∂ X u , such that x = F X . In this work, theorigin of the microscopic coordinates is chosen as the geometrical center of the RVE, with (cid:90) B X d V = , (1)where (cid:82) B d V is the volume integral over the microscopic reference body B . Note that this choiceof origin has no influence on the results but it simplifies the notation. Now the microscopicdeformation x can be split into the sum of terms, x = u + F X + (cid:101) u . (2)Herein, two terms result directly from the macroscale: a constant part u , which describes themacroscopic rigid body translations, and a homogeneous part F X , defined in terms of the macro-scopic deformation gradient. F , u , ¨ F , ¨ uP , f ρ , A P,F , A P,u , A f,F , A f,u s o l v ee q u ili b r i u m , a pp l y c o n s t r a i n t s BS Macroscale Microscale l macro l micro Figure 1.:
Overview of the FE framework in-cluding microscale dynamics. Here, an example ofmacroscopic impact on SHCC is illustrated. X xF = + Hx = X + u B S (a)
Macroscale
X xF = F + (cid:102) Hx = u + F X + (cid:101) u B S (b)
Microscale
Figure 2.:
Large-strain continuum mechanicson both scales.. Tamsen, D. Balzani 4
The difference of these homogeneous deformations u + F X to the actual deformations x is themicroscopic displacement fluctuation field (cid:101) u . This is the field the microscopic BVP will be solvedfor. Now the microscopic displacements u can be written as u = u + HX + (cid:101) u . (3)Analogously, the microscopic deformation gradient can be split as F = F + (cid:102) H with (cid:102) H = ∂ X (cid:101) u . (4)At the macroscale the kinematics are standard, see Figure 2, where the kinematic relations forboth the micro- and macroscale are illustrated. To expand quasi-static homogenization frameworks to the case of dynamics, an extended versionof the Hill-Mandel condition of macro homogeneity [18, 31], which takes into account inertia andbody forces at the microscale, is adopted P : δ F − f · δ u = (cid:104) P : δ F − f · δ u (cid:105) , (5)see [3] and [13]. Herein, (cid:104)•(cid:105) = V (cid:82) B • d V is an expression used to abbreviate the volume averageof a microscopic quantity. P is the first Piola-Kirchhoff stress tensor, f is the microscopic bodyforce vector, and f is its macroscopic counterpart. The variational equation (5) is also calledthe Principle of Multiscale Virtual Power. It ensures that the virtual work of the macroscalecoincides with its respective microscopic volume average, thus ensuring energetic consistencyacross the scales. By decomposing δ F and δ u , following (3) and (4), three important equationscan be derived. First, one obtains (cid:68) P : δ (cid:101) F − f · δ (cid:101) u (cid:69) = 0 , (6)which is automatically fulfilled for mechanical equilibrium. Then, more importantly, the averagingequation for the effective macroscopic stress P and the effective macroscopic body force vector f are derived as P = (cid:104) P − f ⊗ X (cid:105) and (7) f = (cid:104) f (cid:105) . (8)These averaging relations can be found in multiple frameworks dealing with homogenization ofdynamics, see e.g., [13, 27, 28, 37, 42, 43]. It was shown in [26] that the extended Hill-Mandelaveraging relation can be applied in a discretized setting without introducing additional error byscale transition. A principal ingredient of the Hill-Mandel condition, as well as the presented extension (5), is theassumption of scale separation. This means that the equation only holds if the length scale of themicroscopic mechanical fields is significantly smaller than the one of the macroscopic mechanicalfields. For dynamic homogenization this means in practice, that the macroscopic wavelength issufficiently larger than the size of the RVE. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 5
In this paper, an implicit numerical time integration method of first order is applied. Consider-ing ¨ • = ( α • − α ˆ • )∆ t , (9)we can express the second time derivative of any quantity • as the difference of the current timestep and the last ˆ • , divided by the square of the time step ∆ t . The parameters α and α definethe specific time integration scheme. For an explicit time integration, α can be set to zero.Applying this to derivatives of acceleration terms with respect to a quantity (cid:3) in the currentconfiguration leads to ∂ ¨ • ∂ (cid:3) = 1∆ t ∂ ( α • − α ˆ • ) ∂ (cid:3) = α ∆ t ∂ • ∂ (cid:3) . (10)Analogously, the derivative of a value • with respect to the second time derivative ¨ (cid:3) reads ∂ • ∂ ¨ (cid:3) = (cid:32) ∂ ¨ (cid:3) ∂ • (cid:33) − = (cid:32) t α ∂ ( (cid:3) − α ˆ (cid:3) ) ∂ • (cid:33) − = ∆ t α ∂ • ∂ (cid:3) . (11)For the numerical simulations in the numerical analysis section the widely used Newmark method[41] is applied, with α = β . Herein, β is one of the two parameters of the Newmark methodinfluencing the type and stability of the time integration.
3. The Microscopic Problem
We start with the formulation of the microscopic problem, which includes the necessary kinematiclinks to the macroscale. Furthermore, an algorithmic treatment for the associated constraintconditions is given.
For a dynamic analysis, the microscopic balance of linear momentum is given by
Div P + f = . (12)The body force vector f can be decomposed into an inertia part f ρ and a body force vectorrepresenting e.g. the gravitational pull f b . As this framework is intended to model impactloading, gravitational forces are assumed to be negligible compared to the inertia forces. Thus,the relevant body force vector is defined as f := f ρ = − ρ ¨ u with ρ referring to the densityof the microscale constituents in the reference configuration. If gravitational forces have to beconsidered, they can be included in the standard way by the additional force vector f b = ρ g ,where g is the gravitation field. Since this force however, does not depend on the displacements,it does not represent any specialty with view to the proposed homogenization framework and isthus omitted from the presentation to avoid unnecessary complications.Using standard FE procedures for the discretization and linearization of the weak form of thebalance of linear momentum, the well-known equation δ (cid:101) D T (cid:99) K ∆ (cid:101) D = δ (cid:101) D T (cid:98) R (13)is obtained, where (cid:99) K and (cid:98) R are the global tangent stiffness matrix and the residuum matrix,respectively. In the following, the hat (cid:98) • is used to highlight quantities that include dynamic . Tamsen, D. Balzani 6 terms. After incorporating the Dirichlet boundary conditions, the global matrix of incrementalnodal displacements ∆ (cid:101) D are computed from (cid:99) K ∆ (cid:101) D = (cid:98) R in each Newton iteration step in orderto obtain the updated displacements until convergence of the Newton scheme is achieved, i.e.until | ∆ (cid:101) D | < tol . For the classical scheme, the global tangent stiffness matrix (cid:99) K is assembledfrom the element matrix defined as (cid:98) k e = k e + α ∆ t m e , with (14) k e = (cid:90) B e B e T A B e d V and m e = (cid:90) B e N e ρ N e T d V. (15)Herein, N e is the classical element matrix of shape functions, B e denotes the classical B-matrixcontaining the derivatives of the shape functions, and A is the matrix representation of thematerial tangent modulus, defined as the sensitivity of the microscopic stress with respect to themicroscopic deformation gradient as A = ∂ F P . Analogously, the global residuum matrix (cid:98) R isobtained by the assembly of the element-wise counterparts, given as (cid:98) r e = (cid:90) B e (cid:16) B e T P + N e ρ ¨ u (cid:17) d V. (16)Herein, P denotes the matrix representation of the first Piola-Kirchhoff stresses. It can be notedthat both, (cid:98) k e and (cid:98) r e have dynamic terms related to the density ρ , which directly enables theevaluation of inertia at the microscale. As depicted in Figure 1, the macroscopic displacements and deformation gradient and their timederivatives are used to define boundary conditions on the RVE. Inserting them into (2) is only thefirst step. If no additional constraint is considered, the BVP will find an equilibrium where thefluctuations (cid:101) u oppose the applied displacements which results in zero effective displacement of themicrostructure. This might not seem obvious at first, however without any imposed constraint,the energetically most favorable position of each node will be its initial configuration, as anydeviation from it requires energy. Thus, resulting in a microscopic displacement fluctuation fieldof (cid:101) u = u + HX , c.f. (3). Based on the principal of kinematic admissibility, described in detailin [3] and applied in dynamic settings e.g. in [13, 43], two kinematic links are chosen here, i.e. F = (cid:104) F (cid:105) and u = (cid:104) u (cid:105) . (17)The first constraint (17) is well known from quasi-static RVE homogenization frameworks. Itpostulates, that the volume average of the microscopic deformation gradients must equal themacroscopic deformation gradient. This is usually enforced by choosing appropriate boundaryconditions, e.g., linear displacement or periodic boundary conditions. The second constraint(17) is a necessary expansion for the dynamic microscopic problem. This link to the macro-scopic displacements is essential in order to prevent the RVE from moving arbitrarily in space. Inquasi-static calculations, fluctuations e.g. of a corner node in the RVE are restricted, which doesnot influence the results. This is, however, not directly possible for the dynamic case withoutartificially restricting the fluctuations. This is due to the fact that the microscopic deformationsare already influenced by just moving the RVE. Thus, restricting the movement at selected loca-tions will yield different deformation fields. Some dynamic homogenization frameworks, e.g. theone in [42], apply a displacement constraint only on the boundary, which can be a reasonableassumption for dynamic metamaterials. There, the boundary lies in the matrix phase, whichbehaves quasi statically compared to the dynamically active inclusions. The framework proposedhere does not make any a priori assumptions on which part of the microstructure will be dynam-ically significant, as this cannot always be determined in advance for arbitrary problems. Using arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 7 the displacement split (3) and the definition of the origin of the local coordinate system as thecenter of the volume (1), the displacement constraint (17) can be reduced to (cid:104) (cid:101) u (cid:105) = , (18)which states that the constraint is fulfilled, once the volume average of the fluctuations equalszero. To enforce the displacement constraint in (18) on the whole RVE domain, we propose to usethe method of Lagrange multipliers [2]. Similar applications can be found in [4, 43]. For thispurpose, the mechanical boundary value problem is recast in terms of the principle of minimumpotential energy. By adding the potential Π λ associated with the Lagrange multipliers λ and theconstraint (18) to the function of potential energy Π , one obtains Π = Π int + Π ext + Π λ , with Π λ = λ · (cid:90) B (cid:101) u d V. (19)In the following, only the terms concerning the Lagrange multiplier will be regarded, as theother terms capturing the internal potential energy Π int and the external potential energy Π ext will result in the standard FE formulation (13) described above. However, note that due to theLagrange term, the additional degrees of freedom λ appear. The potential energy is varied once with respect to the displacement fluc-tuations (cid:101) u and once with respect to the Lagrange multipliers λ , i.e. δ (cid:101) u Π λ = λ · (cid:90) B δ (cid:101) u d V and (20) δ λ Π λ = δ λ · (cid:90) B (cid:101) u d V. (21) Using (cid:101) u ≈ N e (cid:101) d e and δ (cid:101) u ≈ N e δ (cid:101) d e as FE approximations, the dis-cretized expressions can be written as δ (cid:101) u Π λ = λ T A e (cid:20)(cid:90) B e N e d V δ (cid:101) d e (cid:21) and (22) δ λ Π λ = δ λ T A e (cid:20)(cid:90) B e N e d V (cid:101) d e (cid:21) . (23)To obtain the equivalent of the global volume integral in terms of the elements, the assemblyoperator A is applied for the respective matrices. For better readability, a new element matrixis defined as g e T = (cid:90) B e N e d V. (24)This simplifies the formulations to δ (cid:101) u Π λ = λ T A e (cid:104) g e T δ (cid:101) d e (cid:105) and (25) δ λ Π λ = δ λ T A e (cid:104) g e T (cid:101) d e (cid:105) . (26) . Tamsen, D. Balzani 8 To write the whole system of equations as a global prob-lem, the global matrices are defined in terms of the element matrices, i.e. G = A e g e , (cid:101) D = U e (cid:101) d e and δ (cid:101) D = U e δ (cid:101) d e , (27)where A is the afore-mentioned assembly operator and U a unification operator, as the nodedisplacement fluctuations shared by different elements are not added up, but belong to the samedegree of freedom. Now the expressions (25) and (26) can be reformulated in global fields as δ (cid:101) u Π λ = λ T G T δ (cid:101) D = δ (cid:101) D T Gλ and (28) δ λ Π λ = δ λ T G T (cid:101) D = . (29)Since the Lagrange multiplier only appears in Π λ , no terms result from the variation of Π int and Π ext with respect to λ . It follows, that the second expression has to vanish, see (29). In order to solve the nonlinear global system of equations, the Newton-Raphson method is utilized. For that purpose, we not only need the equations in weak formas in (28) and (29) but also their linearizations. They are used to iteratively compute thenodal displacement fluctuations as well as the Lagrange multipliers. Here the definition of the ∆ operator is used. When linearizing a function f ( x ) = 0 at ˆ x as Lin f = f | ˆ x + ∆ f | ˆ x , then ∆ f = ∂f∂x (cid:12)(cid:12) ˆ x ∆ x . Applying this to the weak forms results inLin δ (cid:101) u Π λ = δ (cid:101) D T Gλ + δ (cid:101) D T G ∆ λ and (30)Lin δ λ Π λ = δ λ T G T (cid:101) D + δ λ T G T ∆ (cid:101) D = . (31) From the linearized variations of Π λ we define the global residua R (cid:101) u = − Gλ = e A r (cid:101) u e with r (cid:101) u e = − g e λ and (32) R λ = − G T (cid:101) D = e (cid:88) r λ e with r λ e = − g e T (cid:101) d e . (33)Including all linearized variations of Π int + Π ext + Π λ yields the discrete equation (cid:104) δ (cid:101) D T δ λ T (cid:105) (cid:34) (cid:99) K GG T (cid:35) (cid:34) ∆ (cid:101) D ∆ λ (cid:35) = (cid:104) δ (cid:101) D T δ λ T (cid:105) (cid:34) (cid:98) R + R (cid:101) u R λ (cid:35) (34)as expansion of (13). After including Dirichlet boundary conditions and applying standard argu-ments of variational calculus, the resulting discrete system of equations reads (cid:34) (cid:99) K GG T (cid:35) (cid:34) ∆ (cid:101) D ∆ λ (cid:35) = (cid:34) (cid:98) R + R (cid:101) u R λ (cid:35) . (35)Note that in contrast to K , the new tangent stiffness matrix is not necessarily positive definite,which needs to be taken into account when choosing and setting up a solver. In general, Lagrangemultipliers have the disadvantage of adding new degrees of freedom to the system of equations. Forthe presented displacement constraint, only one extra degree of freedom is added for each spacialdirection This is due to the fact that the constraint is applied on the whole RVE which avoidsthe approximation of the Lagrange multipliers as field variables. Thus, for three-dimensionalproblems, λ will only add three additional degrees of freedom. Compared to the displacementfluctuations which are linked to the nodes and which may thus easily reach thousands of degreesof freedom, the number of three additional degrees of freedom over the whole RVE is negligible,making it computationally cheap. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 9
The constraint related to the deformation gradient (17) , can be derived and applied in thesame manner as just presented for (17) in the previous section. The only change in the finalformulation is that the related g e T (cid:104) F (cid:105) matrix needs to be computed as the volume average of theelement B-Matrix instead of the shape functions, g e T (cid:104) F (cid:105) = (cid:90) B e B e d V. (36)Applying the constraint regarding the deformation gradient on the volume instead of enforcing itusing periodic boundary conditions will lead to minimally invasive boundary conditions enablinge.g. arbitrary damage propagation without artificial restrictions imposed by periodic boundaryconditions. As shown in [12], such minimally invasive boundary conditions result in a softerconstraint compared to periodic boundary conditions. To simplify the numerical examples in thispaper, only the displacement constraint is applied and the constraint related to the deformationgradient is enforced by using periodic boundary conditions.
4. The Macroscopic Problem
For the solution of the dynamic macroscopic boundary value problem, the associated linearized,discretized balance equations are derived. Herein, specific macroscopic tangent moduli appearwhich are consistently derived for the case where the displacement constraint proposed in theprevious section is taken into account.
The complete macroscopic balance of linearmomentum including inertia is given by
Div P + f = . (37)Applying the Galerkin method with a test function δ u on the entire domain B leads to the weakform of linear momentum (cid:82) B δ u T (cid:0) Div P + f (cid:1) d V = 0 . By applying Div( P ) δ u = Div( P T δ u ) − P : Grad δ u and the Gauss theorem (cid:82) B Div( P T δ u ) d V = (cid:82) ∂ B δ u · t d A , the weak form is writtenas G := (cid:90) B δ F : P d V + (cid:90) B δ u T f ρ d V = 0 . (38)Herein, zero traction forces are taken into account at the boundary and δ F = Grad δ u . Analogousto the microscale, only body forces related to inertia, not gravitation, are considered at themacroscale. Thus, the body force vector is set to f := f ρ = (cid:104) f ρ (cid:105) . To solve the weak form of equilibrium by using the standard Newton-Raphson scheme, the linearized balance of linear momentum is obtained asLin G = G + ∆ G = 0 with ∆ G = (cid:90) B δ F : ∆ P d V + (cid:90) B δ u T ∆ f ρ d V. (39)Now the ∆ operator is applied again to ∆ P and ∆ f ρ : ∆ P = ∂ P ∂ F : ∆ F + ∂ P ∂ ¨ u · ∆¨ u and (40) ∆ f ρ = ∂ f ρ ∂ F : ∆ F + ∂ f ρ ∂ ¨ u · ∆¨ u . (41) . Tamsen, D. Balzani 10 Here an interesting property of the two-scale homogenization framework for dynamics is observed.The macroscopic stress not only depends on the deformation gradient but on the acceleration aswell. In turn, the inertia forces can also depend on the deformation gradient in addition to theacceleration. We define the resulting four sensitivities as A P,F = ∂ F P , A P,u = ∂ ¨ u P , A f,F = ∂ F f ρ and A f,u = ∂ ¨ u f ρ . (42)These moduli (42) are inserted into the linearized weak form which results inLin G = (cid:90) B δ F : P d V + (cid:90) B δ u T f ρ d V + (cid:90) B δ F : A P,F : ∆ F d V + (cid:90) B δ F : A P,u · ∆¨ u d V + (cid:90) B δ u T A f,F : ∆ F d V + (cid:90) B δ u T A f,u · ∆¨ u d V. (43) The linearization is of the weak form of the balance of linearmomentum is now discretized in terms of finite elements. First, the linear increment ∆ G = (cid:90) B δ F : A P,F : ∆ F d V + (cid:90) B δ F : A P,u · ∆¨ u d V + (cid:90) B δ u T A f,F : ∆ F d V + (cid:90) B δ u T A f,u · ∆¨ u d V, (44)is discretized using standard FE formulations. Then, to get rid of the dependence on the timederivatives, the numerical time integration in terms of (10) is used, which results in ∆ G = n ele (cid:88) e =1 δd eP (cid:18)(cid:90) B e B eijP A P,Fijmn B emnQ + α ∆ t B eijP A P,uijk N eQk + N eP i A f,Fimn B emnQ + α ∆ t N eP i A f,uik N eQk d V (cid:19) ∆ d eQ . (45)Herein, the matrix representation of the moduli in index notation has been used. Lowercase in-dices refer the spacial dimension n dm , whereas uppercase indices to the total degrees of freedomof a element n edf . Again, standard element B-matrix B e and shape function matrix N e are con-sidered. By extracting the nodal virtual and incremental displacements, this yields the definitionof the full macroscopic element tangent stiffness matrix (cid:98) k eP Q = (cid:90) B e (cid:18) B eijP A P,Fijmn B emnQ + α ∆ t B eijP A P,uijk N eQk + N eP i A f,Fimn B emnQ + α ∆ t N eP i A f,uik N eQk (cid:19) d V. (46)Now the remaining part of the linearization (43), the residuum R , is discretized as (cid:98) R = n ele (cid:88) e =1 (cid:18) δd P (cid:90) B e B eijP P ij d V + δd P (cid:90) B e N eiP f ρi d V (cid:19) . (47)By extracting the nodal virtual displacements, the element residuum is identified as (cid:98) r eP = (cid:90) B e (cid:16) B eijP P ij + N eP i f ρi (cid:17) d V, (48)where again matrix representation and index notation is used. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 11
For the dynamic homogenization framework, four macroscopic tangent moduli (42) need to bedetermined. To obtain the sought-after moduli in closed form, we start with taking the derivativeof the incremental linearized weak form of linear momentum at the microscale with respect to thetwo relevant macroscopic quantities, the deformation gradient F and the acceleration ¨ u . Thenwe derive the moduli by considering the microscopic problem in its equilibrium state. As it willbe shown later, the derivatives of the microscopic fluctuations with respect to the macroscopicdeformation gradient and the acceleration will be required. For their calculation, the incrementallinearized weak forms have to be derived with respect to these two quantities. In order to accountfor the proposed displacement constraint, the associated increments of the weak form of linearmomentum and of the Lagrange multiplier potential with respect to the microscopic displacementfluctuations as well as the Lagrange multipliers, i.e. (20) and (21), are identified as ∆ G (cid:101) u = (cid:90) B δF ij A ijmn ∆ F mn d V + (cid:90) B δ (cid:101) u i ρ ∆¨ u i d V + ∆ λ i (cid:90) B δ (cid:101) u i d V and (49) ∆ G λ = δλ i (cid:90) B ∆ (cid:101) u i d V. (50)Using the decompositions (3) and (4), equation (49) can be reformulated as ∆ G (cid:101) u = (cid:90) B δF ij A ijmn ∆ F mn d V + (cid:90) B δF ij A ijmn ∆ (cid:101) H mn d V + (cid:90) B δu i ρ ∆¨ u i d V + (cid:90) B δu i ρ ∆ ¨ F ij X j d V + (cid:90) B δu i ρ ∆¨ (cid:101) u i d V + ∆ λ i (cid:90) B δ (cid:101) u i d V. (51) By taking the derivatives of the incre-ments (50) and (51) in the equilibrium state ∆ G = 0 , a closed form formulation of the tangentmoduli will be obtained later. Thus, the associated derivatives are computed in the following. Derivative with Respect to F : Taking the derivatives of (50) and (51) with respect to themacroscopic deformation gradient, while considering (10) results in kl = (cid:90) B δF ij A ijkl d V + (cid:90) B δF ij A ijmn ∂ (cid:101) H mn ∂F kl d V + α ∆ t (cid:90) B δu k ρ X l d V + α ∆ t (cid:90) B δu i ρ ∂ (cid:101) u i ∂F kl d V + ∂λ i ∂F kl (cid:90) B δ (cid:101) u i d V and (52) kl = δλ i (cid:90) B ∂ (cid:101) u i ∂F kl d V. (53)Using standard FE discretization yields kl = n ele (cid:88) e =1 δ (cid:101) d eP (cid:32)(cid:90) B e B eijP A ijkl d V + (cid:90) B e B eijP A ijmn B emnQ d V ∂ (cid:101) d eQ ∂F kl + α ∆ t (cid:90) B e N eP k ρ X l d V + α ∆ t (cid:90) B e N eP i ρ N eQi d V ∂ (cid:101) d eQ ∂F kl + (cid:90) B N eP i d V ∂λ i ∂F kl (cid:33) and (54) kl = n ele (cid:88) e =1 δλ i (cid:32)(cid:90) B N eP i d V ∂ (cid:101) d eP ∂F kl (cid:33) . (55) . Tamsen, D. Balzani 12 Rewriting this in global notation using the abbreviations defined in Appendix A Table 2 leads tothe expressions = L + α ∆ t Z + (cid:16) K + α ∆ t M (cid:17) ∂ (cid:101) D ∂ F + G ∂ λ ∂ F and (56) = G T ∂ (cid:101) D ∂ F . (57)By combining the nodal fluctuations and the Lagrange multipliers into one column matrix D ∗ ,the two equations can be written as one system of equations K ∗ ∂ F D ∗ = − L ∗ with the matrices D ∗ T = (cid:104) (cid:101) D T λ T (cid:105) , (58) L ∗ T = (cid:104) L T + α ∆ t Z T (cid:105) and (59) K ∗ = (cid:34) K + α ∆ t M GG T (cid:35) . (60)Then, the required derivative can be computed from ∂ D ∗ ∂ F = − K ∗ − L ∗ . (61)Note that K ∗ is the microscopic tangent stiffness matrix in (35), which is already available fromsolving the microscopic boundary value problem. Derivative with Respect to ¨ u : Analogously, the derivative of (51) with respect to ¨ u can beobtained by applying (11), i.e. one obtains k = ∆ t α (cid:90) B δF ij A ijmn ∂ ¨ (cid:101) H mn ∂ ¨ u k d V + (cid:90) B δu k ρ d V + (cid:90) B δu i ρ ∂ ¨ (cid:101) u i ∂ ¨ u k d V + ∂λ i ∂ ¨ u k (cid:90) B δ (cid:101) u i d V (62) k = δλ i (cid:90) B ∂ (cid:101) u i ∂ ¨ u k d V. (63)Standard FE discretization using matrix representation and index notation yields k = n ele (cid:88) e =1 δ (cid:101) d eP (cid:32)(cid:90) B e B eijP A ijmn B emnQ d V ∂ (cid:101) d eQ ∂ ¨ u k + (cid:90) B e N P k ρ d V + α ∆ t (cid:90) B e N eP i ρ N eQi d V ∂ (cid:101) d eQ ∂ ¨ u k + (cid:90) B N eP i d V ∂λ i ∂ ¨ u k (cid:33) and (64) k = n ele (cid:88) e =1 δλ i (cid:32)(cid:90) B N eP i d V ∂ (cid:101) d eP ∂ ¨ u k (cid:33) . (65)Again, the two equations are combined by joining the displacements and the Lagrange multipliersinto a single matrix, c.f. (58). By solving the resulting system of equations with respect to therequired derivatives ∂ ¨ u D ∗ , we obtain ∂ D ∗ ∂ ¨ u = − K ∗ − W ∗ . (66)Herein, the definitions (58)-(60), the abbreviations defined in Appendix A Table 2, as well as W ∗ T = (cid:104) W T (cid:105) are used. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 13
In this subsection the four moduli will be derivedby inserting the derivatives computed in the last subsection. Note that all moduli are onlyconsistent for a microscopic equilibrium state. Thus, quadratic convergence of the macroscopicNewton-Raphson iteration is only ensured if the microscopic boundary value problem is solvedfor each macroscopic iteration step. After the last microscopic iteration, the consistent tangentmoduli can be computed.
Derivation of A P,F : To derive the sensitivity of the macroscopic stresses with respect to themacroscopic deformation gradient, the derivative is rewritten using the definition of the macro-scopic stresses in terms of the microscopic fields, i.e. A P,Fijmn = ∂P ij ∂F mn = ∂ (cid:104) P ij + ρ ¨ u i X j (cid:105) ∂F mn . (67)Using the chain rule ∂ P ( F ) ∂ F = ∂ P ∂ F : ∂ F ∂ F , applying (3), (4), (10) and inserting FE discretization,the equation can be written as A P,Fijmn = n ele (cid:88) e =1 (cid:18) V (cid:90) B e A ijmn d V + α ∆ t V (cid:90) B e ρ δ im X j X n d V + 1 V (cid:90) B e A ijkl B eklP d V ∂ (cid:101) d eP ∂F mn + α ∆ t V (cid:90) B e ρ N eP i X j d V ∂ (cid:101) d eP ∂F mn (cid:33) . (68)By using the global abbreviations defined in Appendix A Table 2, Table 3 and inserting (61), theclosed form result is obtained as A P,F = (cid:28) A + 1 β ∆ t Y (cid:29) − V L ∗ T K ∗ − L ∗ . (69)Note that this result has already been presented in [48] for a special scenario of dynamic homog-enization, which did not include macroscopic acceleration and the displacement constraint. Derivation of A P,u : The derivation of the sensitivity of the macroscopic stresses with respect tothe macroscopic accelerations is analogous to that of A P,F . First the derivative is rewritten usingthe definition of the macroscopic stresses in terms of the microscopic fields as A P,uijk = ∂P ij ∂ ¨ u k = ∂ (cid:104) P ij + ρ ¨ u i X j (cid:105) ∂ ¨ u k . (70)Then using the chain rule, (3), (4), (10) and FE discretization, the equation reads A P,uijk = n ele (cid:88) e =1 (cid:32) V (cid:90) B e ρ δ ik X j d V + 1 V (cid:90) B e A ijmn B emnP d V ∂ (cid:101) d eP ∂ ¨ u k + α ∆ t V (cid:90) B e ρ X j N eP i d V ∂ (cid:101) d eP ∂ ¨ u k (cid:33) . (71)Finally, using the global abbreviations in Appendix A Table 2, Table 3 and inserting (66), themodulus is obtained as A P,u = (cid:104) V (cid:105) − V L ∗ T K ∗ − W ∗ . (72) . Tamsen, D. Balzani 14 Macroscopic problemloop over all macroscopic elements loop over element Gauss points
Microscopic probleminput: F , u , ¨ F , ¨ u boundary conditions: x = u + F X + (cid:101) u with (cid:104) (cid:101) u (cid:105) = on B Eq. (2),(18) loop micro-iteration until | ∆ D ∗ | < tol compute global microscopic fields: (cid:99) K , (cid:98) R , R (cid:101) u , R λ , G Eq. (14),(16), (32),(33),(24) solve K ∗ ∆ D ∗ = R ∗ Eq. (35) update D ∗ ⇐ D ∗ + ∆ D ∗ compute homogenized fields and moduli: P , f ρ , A P,F , A P,u , A f,F , A f,u Eq. (7),(8),(69),(72),(75),(78) compute element matrices (cid:98) k e and (cid:98) r e , using Gauss integration Eq. (46),(48) solve (cid:99) K ∆ D = (cid:98) R Figure 3.:
Algorithm for single macroscopic iteration of the dynamic FE framework withrespective equation references. It should be noted that the overall structure of the standardFE procedure does not change, only some additional fields need to be computed. Furthermore,for the implementation of the microscopic problem, the macroscopic displacments u may beomitted from the code. It is the second derivative ¨ u , computed in the macroscopic problem,which influences the microscopic results. Derivation of A f,F : The derivation of the sensitivity of the macroscopic inertia with respectto the macroscopic deformation gradient is again similar to that of A P,F . First, the derivative isrewritten as A f,Fimn = ∂f ρi ∂F mn = ∂ (cid:104) ρ ¨ u i (cid:105) ∂F mn (73)and by using (3), (4), (10) and FE discretization, the equation reads A f,Fimn = n ele (cid:88) e =1 (cid:32) α ∆ t V (cid:90) B e ρ δ im X n d V + α ∆ t V (cid:90) B e ρ N eP i d V ∂ (cid:101) d eP ∂F mn (cid:33) . (74)Using the global abbreviations in Appendix A Table 2, Table 3 and plugging in (61), the modulusis identified as A f,F = α ∆ t (cid:10) V T (cid:11) − V α ∆ t W ∗ T K ∗ − L ∗ . (75) Derivation of A f,u : Analogously, the derivative is rewritten as A f,uik = ∂f ρi ∂ ¨ u k = ∂ (cid:104) ρ ¨ u i (cid:105) ∂ ¨ u k . (76)Then, using (3) and FE discretization, the expression becomes A f,uik = n ele (cid:88) e =1 (cid:32) V (cid:90) B e ρ δ ik d V + 1 V α ∆ t (cid:90) B e ρ N eP i d V ∂ (cid:101) d eP ∂ ¨ u k (cid:33) . (77) arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 15
Taking into account the global abbreviations in Appendix A Table 2, Table 3 and inserting (66),the modulus is derived as A f,u = (cid:104) ρ (cid:105) − V α ∆ t W ∗ T K ∗ − W ∗ . (78)Note that if α and α are equal to zero, which would be equivalent to a quasi-static calculation,the first tangent moduli in (69) take the same form as e.g. found in [33]. Here, the closed formmoduli (69),(72), (75) and (78) extend this consistently to the dynamic regime. An overview overthe algorithm of the proposed framework is presented in Figure 3.
5. Numerical Examples
This section presents numerical examples as a proof of concept, as well as an initial analysis ofdifferent RVE choices. As it turns out, for dynamic homogenization the definition of RVE is evenmore complex than for quasistatic cases. Single-scale comparisons are calculated to asses thereliability of the homogenization framework. First, a rather arbitrary example is shown to ana-lyze the macroscopic Newton iteration and demonstrate the quadratically converging algorithm,which is based on the tangent moduli derived in Section 4. Then the concept of a unit cell asRVE is analyzed under dynamic conditions. Finally, a comparison of two different displacementconstraints, including the proposed one, is presented. All numerical examples make use of theNewmark scheme with the parameters γ = 0 . and β = 0 . , resulting in an unconditionallystable algorithm. For more details on time integration methods in the context of nonlinear FEmethods see e.g. [53].A one-dimensional model of a layered structure with the total length of L is investigated. Thestudied heterogeneous material consists of two alternating phases, a soft and light phase, anda stiff and heavy phase. Each phase has a length of l M , a Young’s modulus E and E , and adensity ρ and ρ , respectively. All calculations are run using E = 2 · Nmm , E = 2 · Nmm , ρ = 1 · kgm and ρ = 1 · kgm . The Poisson’s ratio is chosen to be negligible, i.e. ν = 10 − ,to enable a quasi-1D investigation. The left boundary is fixed, on the right end an impactload is applied in terms of a displacement boundary condition using the polynomial function u ( t ) = u max T ( t ) ( t − T ) , where u max is the amplitude of the impact wave and T the durationin which the load is applied. Initially, the bar is at rest. The problem will be solved using both,a standard single-scale finite element problem referred to as direct numerical simulation (DNS),as well as the proposed dynamic FE framework, c.f. Figure 4 respectively. The DNS discretizesthe microscopic phases at the macroscale into a large number of finite elements with a length Ll M l E u ( t ) E , ρ E , ρ Direct Numerical Simulation (DNS)(a) l E example RVEIntegration Point l E u ( t ) Multiscale Simulation(b)
Figure 4.:
Illustration of the numerical calculations including (a) a 1D single-scale FE Modeland (b) the macroscopic model and the RVE of the FE approach.. Tamsen, D. Balzani 16 of l E . It thereby serves as overkill reference for the multiscale approach. The FE simulationshave a macroscopic element length of l E and make use of the same element length l E at themicroscale for better direct comparability of the microscopic fields to the DNS. To approximatethe displacement fields of the elements, linear shape functions and two Gauss points are used forall scales. As shown in Figure 4(b), each microscopic RVE calculation is associated to a singlemacroscopic integration point. The corresponding parameters for each simulation, regardinggeometry, material parameters and loading will be listed in the caption of each figure. This example analyzes the convergence behavior of the macroscopic Newton iteration. In Fig-ure 5(a), the distribution of macroscopic displacement fields is shown at three different timeinstances for both, the DNS and the FE calculation. As RVE, the basic unit cell of the type A(c.f. Figure 6) is used. It can be seen, that the dynamic multiscale framework approximates theoverall behavior well and even captures some of the smaller waves arising from the microstructure.A better representation of the wave propagation might be achieved by using finer time steps, butthis would generally make the calculation converge faster as the initial values are already closer tothe solution, defying the objective to properly test the tangent moduli. The convergence behaviorfor the three arbitrarily chosen time frames is depicted in Figure 5(b). Quadratic convergenceof the norm of the updates of the nodal displacements | ∆ D ∗ | is observed. This demonstratesthat the macroscopic tangent moduli, incorporating both the microscale inertia forces as well aspossible constraints, have been derived in a consistent manner. − − Coordinate X in mm D i s p l a c e m e n t u ( X ) i n mm t = 0 . s DNS t = 0 . s FE t = 0 . s DNS t = 0 . s FE t = 0 . s DNS t = 0 . s FE (a) − − − − Number of Iterations | ∆ D ∗ | t = 0 . s t = 0 . s t = 0 . s (b) Figure 5.:
Analysis of algorithmic consistency: (a) Comparison of displacement fields and (b)Convergence of the macroscopic Newton iteration with a tolerance of − . The simulationparameters are L = 10000 mm, l M = 10 mm, l E = 33 . mm, u max = 100 mm, T = 0 . s, ∆ t = 5 · s and basic unit cell type A as RVE. Unit Cell Type A Unit Cell Type B
Figure 6.:
Selection of RVE choices with different numbers of basic unit cells.arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 17
For quasi-static homogenization simulations of periodic microstructures, it is known that theresulting macroscopic answer, as well as the corresponding microscopic fields are invariant withrespect to the specific choice of unit cell, as long as an admissible periodic unit cell is chosen.In contrast to quasi-static cases, the distribution of the mass relative to the geometrical cen-ter matters in a dynamic setting. An extreme example is shown in Figure 7, which comparesthe macroscopic displacement field at t = 0 . s presented in the first example in Figure 5awith a simulation using the basic unit cell type B as an RVE (c.f. Figure 6). To properlymeasure the influence of different RVE choices on the FE simulation, an objective error mea-sure (cid:15) is considered. It is defined as average difference of the macroscopic displacement fields (cid:15) = (cid:80) n nodes i (cid:12)(cid:12) u I i ( t j ) − u II i ( t j ) (cid:12)(cid:12) /n nodes . This measure can be evaluated for each time step and thus,the average is once more computed over the number of time steps (cid:15) time = (cid:80) n timesteps j (cid:15) j /n timesteps .Figure 8 shows the calculated errors of different choices of the RVE. For the comparison, RVEswith multiple periods of the same unit cell type, as depicted in Figure 6, were considered. Twoeffects can be observed: The first, presented in Figure 8a, is that the difference in macroscopicdisplacements between different choices of unit cell type decreases, when the number of unit cellsper RVE is increased. This means that the choice of particular basic unit cell type does notmatter as long as the RVE is chosen large enough. The second effect, shown in Figure 8b, isthat the error, computed as difference to the DNS reference, increases when the size of the RVE,relative to the macroscopic element length, gets too large. Then, errors resulting from a violationof the scale separation assumption are obtained. Generally, the second effect can be neglected,as calculations with RVE sizes larger than the macroscopic element length have little practicalapplication when using FE . At this point it is favorable to use domain decomposition approachesinstead of a homogenization method in order to avoid the scale separation assumption. Finally the proposed displacement link u = (cid:104) u (cid:105) is analyzed for the examples in Figure 8b, incomparison to the standard displacement link for quasi-static periodic homogenization, where thefluctuations at the RVE corner nodes is set to zero, i.e. (cid:101) u corner = . For the quasi-1D example Coordinate X in mm D i s p l a c e m e n t u ( X ) i n mm FE - Basic Unit Cell AFE - Basic Unit Cell B Figure 7.:
Comparison of macroscopic displacement fields for two different RVEs with dif-ferent basic unit cells types. The simulation parameters are L = 10000 mm, l M = 10 mm, l E = 33 . mm, u max = 100 mm, T = 0 . s, t = 0 . s and ∆ t = 5 · s.. Tamsen, D. Balzani 18 . . . . Basic Unit Cells per RVE E rr o r (cid:15) t i m e i n mm (a) . . . . Basic Unit Cells per RVE E rr o r (cid:15) t i m e i n mm Unit Cells Type AUnit Cells Type B (b)
Figure 8.:
Analysis of different RVE choices: (a) Direct comparison of RVEs with unit cell typeA and B shown by with increasing number of basic unit cells per RVE (the error is computedas difference between the response of the two unit cell types, not with respect to the DNS). (b)Error of unit cell type A and B, compared to DNS as reference. The simulation parametersare L = 10000 mm, l M = 2 . mm, l E = 20 mm, u max = 100 mm, T = 0 . s, ∆ t = 5 · s, n timesteps = 400 . analyzed here, this is equivalent to setting the integral over the surface equal to the correspondingmacroscopic displacements, which has been taken into account in other dynamic homogenizationschemes.The first observation is, that using the proposed volume constraint u = (cid:104) u (cid:105) results in a morerobust framework in terms of stability of the Newton-Raphson iterations. Furthermore, slightlysmaller error values are obtained compared to the DNS reference. Table 1 shows the number oftime steps reached before either the calculations crashed (due to diverging Newton iterations atthe microscale) or they were finished successfully after the intended complete set of 1000 timesteps. Especially the calculations using the unit cell type B in combination with the zero fluc-tuations of the corner nodes, underperformed the other scenarios. To understand the differencebetween the performance of the displacement links, it is necessary to examine the behavior atthe RVE level.Here the examples with the RVEs consisting of three periods of the basic unit cell B are fur-ther analyzed. Figure 9 compares the microscopic displacements for four relevant time instancesNumber of Unit Cells per RVE 1 3 5 7Unit Cell Type u -Link Number of Time StepsA u = (cid:104) u (cid:105) (cid:101) u corner = B u = (cid:104) u (cid:105) (cid:101) u corner =
944 420 192 166
Table 1.:
Number of time steps before either the simulation crashed (divergence of Newtoniteration at microscale) or the intended complete set of 1000 time steps was successfully reached.Different choices of RVEs and constraints were analyzed.arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 19 right before the peak of the input wave passes through the RVEs. More specifically, the differ-ences between the microscopic displacement fields u of an RVE and the respective macroscopicdisplacements u . To compare the DNS, an effective u has been computed as the average displace-ment over the associated length. Thereby, the quality of the microscopic fields can be analyzedindependently from the macroscopic displacements. With this , the two different displacementconstraint options can be effectively compared with the reference solution obtained from DNS.The graphs show, that the fixed corner constraint leads to artificially increased displacementintensities at the microscale due to the constricted boundary. These increased displacementseventually lead to extreme deformations in single elements at the microscale, crashing the simu-lation. The proposed displacement volume constraint however, leads to a softer constraint whichresults in a more robust computation while still enabling dynamic effects which agree well withthe ones from the reference DNS. In the presented examples, the only rate-dependent influenceare inertia forces. In cases where also rate-dependent material properties are included we expectthe influence of different displacement constraints on the overall simulation to increase, in favorfor the proposed volume constraint. − . − . − .
25 6 . . . − . − . − . . . . Position X of RVE in mm u − u i n mm DNS FE VC FE FC u = − . mm u = − . mm u = − . mm (a) t = 0 . s − . − . − .
25 6 . . . − . − . − . − . . . . Position X of RVE in mm u − u i n mm DNS FE VC FE FC u = − . mm u = − . mm u = − . mm (b) t = 0 . s − . − . − .
25 6 . . . − . − . . . Position X of RVE in mm u − u i n mm DNS FE VC FE FC u = − . mm u = − . mm u = − . mm (c) t = 0 . s − . − . − .
25 6 . . . − . − . . . . Position X of RVE in mm u − u i n mm DNS FE VC FE FC u = − . mm u = − . mm u = − . mm (d) t = 0 . s Figure 9.:
Comparison of microscopic displacement fields obtained from the FE simulationswith the DNS. The displacements have been normalized by u (for the DNS, average displace-ment) to analyze the quality of microscopic displacements more or less independent of themacroscopic displacements. The two different displacement links, the volume constraint (VC)and the fixed corners (FC) are analyzed. The simulation parameters are L = 10000 mm, l M = 2 . mm, l E = 20 mm, u max = 100 mm, T = 0 . s, ∆ t = 5 · s, location of themacroscale integration point X = 7504 , mm, section of the DNS displacement field from X = 7496 . mm to X = 7511 . mm, RVE with three periods of basic unit cell types B.. Tamsen, D. Balzani 20
6. Conclusion
In this paper, a general purpose, consistent, two-scale homogenization framework for dynamicsat the macro- and microscale was proposed in the sense of the FE method. The frameworkdoes not include any simplifications such as linearized strains, explicit time integration or partlyquasi-static scenarios. The only assumption taken into account is a sufficiently pronounced scaleseparation, which is anyway essential requirement of any FE approach. Therefore, it enablesthe simulation of various structural problems of complex micro-heterogeneous materials underdynamic loading such as impact. The main aspects are: (i) the incorporation of the completebalance of momentum at the micro- and macroscale including inertia forces, (ii) a large-strainformulation enabling the simulation of a wide range of macro- and micro-mechanical phenomena,(iii) extended microscopic boundary conditions resulting in a more robust scheme and giving thepossibility for non-periodic RVE boundaries, and (iv) the derivation of consistent macroscopictangent moduli ensuring quadratically converging macroscopic iterations. The presented resultson different choices of RVEs show that different admissible unit cells of a periodic microstructurelead also to a different mechanical response - even if periodic boundary conditions are employed.This is in contrast to quasi-static scenarios. However, these results should not entail that themost basic unit cell is necessarily a bad choice for a simulation, but the specific choice of thisbasic unit cell is not unique. Future work will focus on the analysis of more complex RVEs fornon-periodic microstructures under dynamic loading. Acknowledgment
The authors gratefully acknowledge financial support from the German Research Foundation(DFG) within project B1 of the Training Research Group (GRK) 2250 “Mineral-Bonded Com-posites for Enhanced Structural Impact Safety”. Furthermore, fruitful discussions with CeliaReina (University of Pennsylvania) are greatly appreciated. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 21
AppendixA. Matrix Abbreviations
Global Element
Definition Size K = n ele A e =1 k e n edf × n edf k eP Q = (cid:90) B e B eijP A ijkl B eklQ d V L = n ele A e =1 l e n edf × n dm l eP ij = (cid:90) B e B eklP A klij d V M = n ele A e =1 m e n edf × n edf m eP Q = (cid:90) B e N eP i ρ N eQi d V W = n ele A e =1 w e n edf × n dm w eP i = (cid:90) B e ρ N eP i d V Z = n ele A e =1 z e n edf × n dm z eP ij = (cid:90) B e ρ N eP i X j d V G = n ele A e =1 g e n edf × n dm g eP i = (cid:90) B e N eP i d V V n dm × n dm V ijk = ρ δ ik X j Y n dm × n dm Y ijkl = ρ δ ik X l X j Table 2.:
Overview of the used fields at element and global level, with n edf : number of DOFat element level and n dm : spacial dimension. Matrix Size D ∗ = (cid:104) (cid:101) D T λ T (cid:105) T ( n edf + n lgr ) × L ∗ = (cid:104) L T + α ∆ t Z T (cid:105) T ( n edf + n lgr ) × n dm L ∗ = (cid:104) L T + α ∆ t Z T (cid:105) T ( n edf + n lgr ) × n dm W ∗ = (cid:104) W T (cid:105) T ( n edf + n lgr ) × n lgr K ∗ = (cid:34) K + α ∆ t M GG T (cid:35) ( n edf + n lgr ) × ( n edf + n lgr ) Table 3.:
Overview of the extended fields, with n edf : number of DOF at element level, n lgr : number of DOF of Lagrange constraint and n dm : spacial dimension.. Tamsen, D. Balzani 22 References [1] A. Bensoussan, J.-L. Lions, and G. Papanicolaou.
Asymptotic analysis for periodic structures .North-Holland Publishing Company, Amsterdam, 1978.[2] D. P. Bertsekas.
Constrained optimization and lagrange multiplier methods . Athena Scientific,1996.[3] P. J. Blanco, P. J. Sánchez, E. A. de Souza Neto, and R. A. Feijóo. Variational foundationsand generalized unified theory of RVE-based multiscale models.
Archives of ComputationalMethods in Engineering , 23:191–253, 2016.[4] P. J. Blanco, A. Clausse, and R. A. Feijóo. Homogenization of the navier-stokes equations bymeans of the multi-scale virtual power principle.
Computer Methods in Applied Mechanicsand Engineering , 315:760–779, 2017.[5] S. Brûlé, E. H. Javelaud, S. Enoch, and S. Guenneau. Experiments on seismic metamaterials:Molding surface waves.
Physical Review Letters , 112:133901, 2014.[6] E. W. C. Coenen, V. G. Kouznetsova, and M. G. D. Geers. Novel boundary conditions forstrain localization analyses in microstructural volume elements.
International Journal forNumerical Methods in Engineering , 90:1–21, 2012.[7] R. V. Craster, J. Kaplunov, and A. V. Pichugin. High-frequency homogenization for periodicmedia.
Proceedings of the Royal Society A , 466:2341–2362, 2010.[8] S. A. Cummer and D. Schurig. One path to acoustic cloaking.
New Journal of Physics , 9:45, 2007.[9] I. Curosu, V. Mechtcherine, and O. Millon. Effect of fiber properties and matrix compositionon the tensile behavior of strain-hardening cement-based composites (SHCCs) subject toimpact loading.
Cement and Concrete Research , 82:23–35, 2016.[10] I. Curosu, V. Mechtcherine, D. Forni, and E. Cadoni. Performance of various strain-hardeningcement-based composites (SHCC) subject to uniaxial impact tensile loading.
Cement andConcrete Research , 102:16–28, 2017.[11] I. Curosu, V. Mechtcherine, M. Hering, and M. Curbach. Mineral-bonded composites forenhanced strucutural impact safety - Overview of the format, goals and achievements of theresearch training group GRK 2250. , 2019.[12] E. A. de Souza Neto and R. A. Feijóo.
Variational foundations of large strain multiscalesolid constitutive models: kinematical formulation , chapter 9, pages 341–378. John Wileyand Sons, Ltd, 2010.[13] E. A. de Souza Neto, P. J. Blanco, P. J. Sánchez, and R. A. Feijóo. An RVE-based multiscaletheory of solids with micro-scale inertia and body force effects.
Mechanics of Materials , 80:136–144, 2015.[14] F. Feyel. Multiscale FE elastoviscoplastic analysis of composite structures. ComputationalMaterials Science , 16:344–354, 1999.[15] J. Fish and K. Shek. Finite deformation plasticity for composite structures: computationalmodels and adaptive strategies.
Computer Methods in Applied Mechanics and Engineering ,172:145–174, 1999. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 23 [16] J. Fish, W. Chen, and G. Nagai. Non-local dispersive model for wave propagation in het-erogeneous media: multi-dimensional case.
International Journal for Numerical Methods inEngineering , 54:347–63, 2002.[17] M. G. D. Geers, V. G. Kouznetsova, K. Matouš, and J. Yvonnet.
Homogenization methodsand multiscale modeling: Nonlinear problems , pages 1–34. John Wiley and Sons, Ltd., 2017.[18] R. Hill. On constitutive macro-variables for heterogeneous solids at finite strain.
Proceedingsof the Royal Society A , 326:131–147, 1972.[19] R. Hu and C. Oskay. Multiscale nonlocal effective medium model for in-plane elastic wavedispersion and attenuation in periodic composites.
Journal of the Mechanics and Physics ofSolids , 124:220–243, 2019.[20] T. Hui and C. Oskay. A high order homogenization model for transient dynamics of hetero-geneous media including micro-inertia effects.
Computer Methods in Applied Mechanics andEngineering , 273:181–203, 2014.[21] M. Kadic, T. Bückmann, R. Schittny, and M. Wegener. Experiments on cloaking in optics,thermodynamics and mechanics.
Philosophical Transactions of the Royal Society A , 373:20140357, 2015.[22] A. Karamnejad and L. J. Sluys. A dispersive multi-scale crack model for quasi-brittle het-erogeneous materials under impact loading.
Computer methods in applied mechanics andengineering , 278:423–444, 2014.[23] C. Kettenbeil and G. Ravichandra. Experimental investigation of the dynamic behavior ofmetaconcrete.
International Journal of Impact Engineering , 111:199–207, 2018.[24] R. Khajehtourian and M. I. Hussein. Dispersion characteristics of a nonlinear elastic meta-material.
AIP Advances , 4:124308, 2014.[25] J. Li and C. T. Chan. Double-negative acoustic metamaterial.
Physical Review E , 70(5):055602, 2004.[26] C. Liu and C. Reina. Discrete averaging relations for micro to macro transition.
Journal ofApllied Mechanics , 83(8):081006, 2016.[27] C. Liu and C. Reina. Discrete averaging relations for micro to macro transition.
Journal ofApplied Mechanics , 83(8):081006, 2016.[28] C. Liu and C. Reina. Variational coars-graining procedure for dynamic homogenization.
Journal of the Mechanics and Physics of Solids , 104:187–206, 2017.[29] C. Liu and C. Reina. Dynamic homogenization of resonant elastic metamaterials with space/-time modulation.
Computational Mechanics , 64(1):147–161, 2018.[30] Z. Liu, X. Zhang, Y. Mao, Y. Y. Zhu, Z. Yang, C. T. Chan, and P. Sheng. Locally resonantsonic materials.
Science , 289(5485):1734–1736, 2000.[31] J. Mandel.
Plasticité classique et viscoplasticité . Springer, Udine, 1971.[32] C. Miehe. Numerical computation of algorithmic (consistent) tangent moduli in large-straincomputational inelasticity.
Computer Methods in Applied Mechanics and Engineering , 134(3):223–240, 1996.[33] C. Miehe, J. Schotte, and J. Schröder. Computational micro-macro transitions and overallmoduli in the analysis of polycrystals at large strains.
Computational Materials Science , 16:372–382, 1999. . Tamsen, D. Balzani 24 [34] G. W. Milton and J. R. Willis. Minimum variational principles for time-harmonic wavesin a dissipative medium and associated variational principles of Hashin-Shtrikman type.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 466(2122):3013–3032, 2010.[35] M. Miniaci, A. Krushynska, F. Bosia, and N. M. Pugno. Large scale mechanical metamate-rials as seismic shields.
New Journal of Physics , 18:083041, 2016.[36] S. J. Mitchell, A. Pandolfi, and M. Ortiz. Effect of brittle fracture in a metaconcrete slabunder shock loading.
Journal of Engineering Mechanics , 142(4), 2016.[37] A. Molinari and S. Mercier. Micromechanical modelling of porous materials under dynamicloading.
Journal of the Mechanics and Physics of Solids , 49:1497–1516, 2001.[38] H. Moulinec and P. Suquet. A numerical method for computing the overall response ofnonlinear composites with complex microstructure.
Computer Methods in Applied Mechanicsand Engineering , 157:69–94, 1998.[39] H. Nassar, Q.-C. He, and N. Auffray. On asymptotic elastodynamic homogenization ap-proaches for periodic media.
Journal of the Mechanics and Physics of Solids , 88:274–290,2016.[40] S. Nemat-Nasser and A. Srivastava. Overall dynamic constitutive relations of layered elasticcomposites.
Journal of the Mechanics and Physics of Solids , 59(10):1953–1965, 2011.[41] N. M. Newmark. A method of computation for structural dynamics.
Journal of EngineeringMechanics , 85:67–94, 1959.[42] K. Pham, V. G. Kouznetsova, and M. G. D. Geers. Transient computational homogenizationfor heterogeneous materials under dynamic excitation.
Journal of the Mechanics and Physicsof Solids , 61:2125–2146, 2013.[43] D. Roca, O. Lloberas-Valls, J. Cante, and J. Oliver. A computational multiscale homog-enization framework accounting for inertia effects: Application to acoustic metamaterialsmodelling.
Computer methods in applied mechanics and engineering , 330:415–446, 2018.[44] C. Sartori, S. Mercier, N. Jacques, and A. Molinari. Constitutive behavior of porous ductilematerials accounting for micro-inertia and void shape.
Mechanics of Materials , 80:324–339,2015.[45] J. Schröder.
Plasticity and Beyond - Microstructures, Crystal-Plasticity and Phase Transi-tions (CISM Lecture Notes 550, Eds. J. Schröder, K. Hackl) , chapter A numerical two-scalehomogenization scheme: the FE -method. Springer, 2013.[46] R. Smit, W. Brekelmans, and H. Meijer. Prediction of the mechanical behavior of nonlinearheterogeneous systems by multi-level finite element modeling. Computer Methods in AppliedMechanics and Engineering , 155:181–192, 1998.[47] A. Sridhar, V. G. Kouznetsova, and M. G. D. Geers. A general multiscale framework for theemergent effective elastodynamics of metamaterials.
Journal of the Mechanics and Physicsof Solids , 111:414–433, 2018.[48] E. Tamsen, W. Weber, and D. Balzani. First steps towards the direct micro-macro simula-tion of reinforced concrete under impact loading.
Proceedings in Applied Mathematics andMechanics , 18(1):e201800181, 2018. arge-Strain FE Framework for the Simulation of Dynamic Problems on Two Scales 25 [49] K. Terada, M. Hori, T. Kyoyac, and N. Kikuchi. Simulation of the multi-scale convergencein computational homogenization approaches.
International Journal of solids and structures ,37:2285–2311, 2000.[50] Z.-P. Wang and C. T. Sun. Modeling micro-inertia in heterogeneous materials under dynamicloading.
Wave Motion , 36:473–485, 2002.[51] J. R. Willis.
Dynamics of composites . Springer Vienna, 1997.[52] J. R. Willis. The construction of effective relations for waves in a composite.
Comptes RendusMécanique , 340(4):181–192, 2012.[53] P. Wriggers.