A general, implicit, finite-strain FE2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} framework for the simulation of dynamic problems on two scales

In this paper we present a fully-coupled, two-scale homogenization method for dynamic loading in the spirit of FE2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} methods. The framework considers the balance of linear momentum including inertia at the microscale to capture possible dynamic effects arising from micro heterogeneities. A finite-strain formulation is adapted to account for geometrical nonlinearities enabling the study of e.g. plasticity or fiber pullout, which may be associated with large deformations. A consistent kinematic scale link is established as displacement constraint on the whole representative volume element. The consistent macroscopic 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 applying existing assembly and solving procedures. Thus, making it suitable for standard finite element program architectures. Numerical examples of a layered periodic material are presented and compared to direct numerical simulations to demonstrate the capability of the proposed framework. In addition, a simulation of a split Hopkinson tension test showcases the applicability of the framework to engineering problems.


Introduction
Under dynamic loading, micro-heterogeneities can give rise to additional wave effects at the microscale leading to complex microscopic stress distributions. This then contributes to the complex macroscopic dynamic material behavior. There are various examples from different fields of application which cover a broad range of length scales. Currently the interest is high in metamaterials in general, but especially in locally resonant metamaterials exhibiting special properties like band gaps and negative bulk moduli. The applications range from cloaking devices [8,24] over tunable sound attenuation [29,34] to earthquake protection [5,43]. More classical materials are being investigated as well. One example is metaconcrete, which replaces aggregates by rubber-coated lead inclusions to weaken impact waves [26,44]. A different B Daniel Balzani daniel.balzani@rub.de 1 approach for an improved impact resistance are strainhardening cement-based composites (SHCC), that show a pronounced energy dissipation under dynamic loading, as well as a change in fiber failure and overall crack pattern [10][11][12]. In addition to that, porous materials have shown an influence of microscopic inertia on voids under high strain rates [45,53]. This list serves only to illustrate the possible influence of the material microstructure on the macroscopic response under high dynamic loading for a wide range of materials and applications. 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 or cracks are present at the microscale, can exhibit distinct effective macroscopic properties resulting from the dynamics in the microstructure. In principle, locally large deformations may occur at the microscale in such materials. One example is the fiber-pullout process in SHCC under dynamic loading. Due to the progress of cracks through the cementitious matrix, structural problems arise at the microscale with deformations which are large compared to the structural components such as the fibers or the crack. This requires the consideration of finite strains in a geometrically nonlinear setting where the description is not limited to linearized strains.
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 in numerical material analysis, see e.g. [15,16,40,46,55,59]. A recent overview of computational homogenization methods in general is given in [18]. In addition a recapitulation of the FE 2 method in particular is presented in [54]. However, with the rise in interest in metamaterials more and more dynamic homogenization frameworks have been published in the last years. Early works date back to elastodynamic problems [62], with more recent work confirming that the effective macroscopic density under the influence of inertia is not the simple average density as for quasi static applications, c.f. [37,41]. Then there are works dealing with microdynamics analytically [64,67], motivating further research. One widely applied homogenization method is asymptotic expansion, see e.g., [7,17,[20][21][22], which is mainly based on the original work of Bensoussan et al. [1]. In addition to that, there is the more general theory of elastodynamic homogenization by Willis [63], which has been applied in e.g. [42,47,49,65] and others. The two approaches are based on related ideas and their similarities are studied in [48]. Both mentioned methods are limited to elastic, periodic media, a considerable limitation when dealing with various composites. A further approach is based on the work of Irving and Kirkwood [23]. It is called the continuum homogenization theory, where the main extensive quantities as mass, momentum and energy are computed as weighted averages of their microscopic counterparts, c.f. [35,38]. A more general approach is the micro-macro simulation based on a representative volume element (RVE). In such homogenization methods, macroscopic quantities such as the deformation gradient and the displacement vector at a macroscopic integration point are projected onto a microscopic boundary value problem whose homogenized response replaces the constitutive material law. In the case where the finite element method (FEM) is used on both scales, this procedure is called the FE 2 method. A good theoretic introduction to this theory including dynamics is given by [14], which has been the foundation on which this paper is built. The framework in [17] calculates a quasi-static microstructure but then applies an inertiainduced eigenstrain based on the microstructure as an extra body force at the macroscale to account for micro-inertia effects. This was extended by [25] to account for matrix cracking at the microscale under impact loading. Other, rather FE 2 -type schemes as [32,33,51,52,56,60,61] calculate the full balance of linear momentum at the microscale. In [32] an explicit, periodic, small-strain framework is presented, which was extended to an implicit time integration method for modeling resonant elastic metamaterials in [33]. [51,52] use the assumption of linear elasticity to improve the computational performance, by splitting the problem into a purely static and a special dynamic boundary value problem (BVP). To better capture a wider range of applied frequencies, [56] use a Floquet-Bloch transformation to build a base of eigenmodes to analyze elastic, periodic metamaterials. The mentioned frameworks all use at least one of the following idealizations: small strains, linear elasticity, periodic or symmetric microstructures. In addition, many require quite elaborate implementations. A recent publication [60] does fulfill the named requirements, however for the derivation of the macroscopic tangent moduli, specific nodes are assumed to experience no microscale fluctuation. This assumption is quite reasonable in the context of metamaterials, whose analysis was primary goal in [60], because there the regions at the unitcell boundaries are particularly stiff. For microscopic problems not specifically addressing classical metamaterial microstructures, however, more flexibility at the boundaries may be more reasonable which is why we propose an alternative approach here.
The aim of this paper is to build a multiscale framework for dynamic loading as general as possible, while still being compatible with standard FE architecture. To enable the analysis of micro-mechanical processes as plasticity or fiber pullout, as well as to incorporate effects resulting from geometric nonlinearities, the proposed framework uses a finite-strain formulation. Its importance is supported by the analytical example in [27] of a nonlinear, elastic metamaterial, where finite strains were shown to be relevant for large wave amplitudes. The conducted numerical studies are so far not specifically chosen to highlight this aspect of the framework. To permit a flexible damage evolution in the RVE, which is not dominated by periodic boundary conditions [6] or aforementioned assumptions regarding fixated fluctuations at specific boundary nodes [60], we propose to apply kinematic scale links as constraints on the whole RVE. This allows us to model any type of RVE morphology.
The paper starts by discussing the general ideas of the framework in Sect. 2. The used notations of the large-strain framework are presented and the averaging relations derived in [13] are briefly recapitulated to include the full balance of momentum at the microscale. Then in Sect. 3 the FE formulations of the microscale are presented and the kinematic constraints are derived. In Sect. 4 the respective macroscale formulations are displayed. There, a focus is laid on the derivation of consistent macroscale tangent moduli in closed form. These enable a quadratically converging macroscopic Newton-iteration, resulting in a robust and efficient algorithm compared to numerical differentiation through perturbation of the macroscopic quantities [39]. Once the whole theoretical background of the framework is explained, Sect. 5 provides numerical examples, demonstrating the convergence behavior and analyzing some properties of RVEs under dynamic loading. In Sect. 6, the applicability of the framework to construction materials is given. A split Hopkinson tension test is conducted on a sample of strain-hardening cementitious composite (SHCC). The presented example is an extension of the example published in [58], presenting a more complex, 3D microstructure. The publication [58] is focused on the comparison to experimental results and the understanding of the composite behavior, without going into detail on the developed homogenization framework. The paper is completed in Sect. 7 with concluding remarks.

Homogenization framework including microscopic inertia
The general idea of the homogenization framework for dynamics is to consider the full balance of linear momentum including the inertia terms at the microscale. This enables not only the analysis of full dynamic fields at the microscale but a direct study of microscopic inertia effects on the macroscale. By using appropriate averaging relations and kinematic links, a consistent scale bridging for dynamic loading is established. In the FE 2 method, each macroscopic Gauss point is associated with a separate microscopic RVE simulation which uses the macroscopic mechanical quantities to define the microscopic BVP. In order to differentiate the two scales, variables associated with the macroscale are denoted with a bar •. Here, the finite-strain framework is taken 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 ingredients of the scale-coupling framework are explained, see also the schematic illustration in Fig. 1.

Kinematics at different scales
The connection of the (undeformed) reference and the (deformed) current configuration is described by the displacement u as u = x − X, where X ∈ B denotes the coordinates in the reference configuration and x ∈ S the deformation. The link between the two configurations in terms of transformations of vector elements is described by the deformation and displacement gradients, respectively F = ∂ X x = 1 + H and H = ∂ X u, such that x = FX. In this work, the origin of the microscopic coordinates is chosen as the geometrical center of the RVE, with where B dV is the volume integral over the microscopic reference body B. Note that this choice of origin has no influence on the results but it simplifies the notation. Now the microscopic deformation x can be split into the sum of terms, Herein, two terms result directly from the macroscale: a constant part u, which describes the macroscopic rigid body translations, and a homogeneous part FX, defined in terms of the macroscopic deformation gradient. The difference of these homogeneous deformations u + FX to the actual deformations x is the microscopic displacement fluctuation field u. This is the field the microscopic BVP will be solved for. Now the microscopic displacements u can be written as Analogously, the microscopic deformation gradient can be split as At the macroscale the kinematics are standard, see Fig. 2, where the kinematic relations for both the micro-and macroscale are illustrated.

Averaging relations
To expand quasi-static homogenization frameworks to the case of dynamics, an extended version of the Hill-Mandel condition of macro homogeneity [19,36], which takes into account inertia and body forces at the microscale, is adopted see [3] and [14]. Herein, • = 1 V B • dV is an expression used to abbreviate the volume average of a microscopic quantity. P is the first Piola-Kirchhoff stress tensor, f is the microscopic body force vector, and f is its macroscopic counterpart. The variational Eq. (5) is also called the Principle of Multiscale Virtual Power. It ensures that the virtual work of the macroscale coincides with its respective microscopic volume average, thus ensuring energetic consistency across the scales. By decomposing δF and δu, following (3) and (4), three important equations can be derived. First, one obtains which is automatically fulfilled for mechanical equilibrium. Then, more importantly, the averaging equation for the effective macroscopic stress P and the effective macroscopic body force vector f are derived as P = P − f ⊗ X and (7) These averaging relations can be found in multiple frameworks dealing with homogenization of dynamics, see e.g., [14,31,32,45,51,52,60]. It was shown in [30] that the extended Hill-Mandel averaging relation can be applied in a discretized setting without introducing additional error by scale transition.

Scale separation
A principal ingredient of the Hill-Mandel condition, as well as the presented extension (5), is the assumption of scale separation. This means that the equation only holds if the length scale of the microscopic mechanical fields is significantly smaller than the one of the macroscopic mechanical fields. For dynamic homogenization this means in practice, that the macroscopic wavelength is sufficiently larger than the size of the RVE.

Time integration
In this paper, an implicit numerical time integration method of first order is applied. Considering we can express the second time derivative of any quantity • as the difference of the current time step and the last•, divided by the square of the time step t, bearing in mind that• includes the first and second time derivatives of the last time step. The parameters α 1 and α 2 define the specific time integration scheme. For an explicit time integration, α 1 can be set to zero. Applying this to derivatives of acceleration terms with respect to a quantity in the current configuration leads to Analogously, the derivative of a value • with respect to the second time derivative¨ reads For the numerical simulations in the numerical analysis section the widely used Newmark method [50] is applied, with α 1 = 1 β . Herein, β is one of the two parameters of the Newmark method influencing the type and stability of the time integration.

The microscopic problem
We start with the formulation of the microscopic problem, which includes the necessary kinematic links to the macroscale. Furthermore, an algorithmic treatment for the associated constraint conditions is given.

Microscopic balance of linear momentum
For a dynamic analysis, the microscopic balance of linear momentum is given by The body force vector f can be decomposed into an inertia part f ρ and a body force vector representing e.g. the gravitational pull f b . As this framework is intended to model impact loading, gravitational forces are assumed to be negligible compared to the inertia forces. Thus, the relevant body force vector is defined as f := f ρ = − ρ 0ü with ρ 0 referring to the density of the microscale constituents in the reference configuration. If gravitational forces have to be considered, they can be included in the standard way by the additional force vector f b = ρ 0 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 is thus omitted from the presentation to avoid unnecessary complications. Using standard FE procedures for the discretization and linearization of the weak form of the balance of linear momentum, the well-known equation is obtained, where K and R are the global tangent stiffness matrix and the residuum matrix, respectively. In the following, the hat • is used to highlight quantities that include dynamic terms. After incorporating the Dirichlet boundary conditions, the global matrix of incremental nodal displacements D are computed from K D = R in each Newton iteration step in order to obtain the updated displacements until convergence of the Newton scheme is achieved, i.e. until | D| < tol. For the classical scheme, the global tangent stiffness matrix K is assembled from the element matrix defined as (15) Herein, N e is the classical element matrix of shape functions, B e denotes the classical B-matrix containing the derivatives of the shape functions, and A is the matrix representation of the material tangent modulus, defined as the sensitivity of the microscopic stress with respect to the microscopic deformation gradient as A = ∂ F P. Analogously, the global residuum matrix R is obtained by the assembly of the element-wise counterparts, given as Herein, P denotes the matrix representation of the first Piola-Kirchhoff stresses. It can be noted that both, k e and r e have dynamic terms related to the density ρ 0 , which directly enables the evaluation of inertia at the microscale.

Kinematic links to the macroscale
As depicted in Fig. 1, the macroscopic displacements and deformation gradient and their time derivatives are used to define boundary conditions on the RVE. Inserting them into (2) is only the first step. If no additional constraint is considered, the BVP will find an equilibrium where the fluctuations u oppose the applied displacements which results in zero effective displacement of the microstructure. 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 any deviation from it requires energy. Thus, resulting in a microscopic displacement fluctuation field of u = u + HX, c.f. (3). Based on the principal of kinematic admissibility, described in detail in [3] and applied in dynamic settings e.g. in [14,52], two kinematic links are chosen here, i.e.
The first constraint (17) 1 is well-known from quasi-static RVE homogenization frameworks. It postulates, that the volume average of the microscopic deformation gradients must equal the macroscopic deformation gradient. This is usually enforced by choosing appropriate boundary conditions, e.g., linear displacement or periodic boundary conditions. The second constraint (17) 2 is a necessary expansion for the dynamic microscopic problem. This link to the macroscopic displacements is essential in order to prevent the RVE from moving arbitrarily in space. In quasi-static calculations, fluctuations e.g. of a corner node in the RVE are restricted, which does not influence the results. This is, however, not directly possible for the dynamic case without artificially restricting the fluctuations. This is due to the fact that the microscopic deformations are already influenced by just moving the RVE. Thus, restricting the movement at selected locations will yield different deformation fields. Some dynamic homogenization frameworks, e.g. the one in [51], apply a displacement constraint only on the boundary, which can be a reasonable assumption for dynamic metamaterials. There, the boundary lies in the matrix phase, which behaves quasi statically compared to the dynamically active inclusions. The framework proposed here does not make any a priori assumptions on which part of the microstructure will be dynamically significant, as this cannot always be determined in advance for arbitrary problems.
Using the displacement split (3) and the definition of the origin of the local coordinate system as the center of the volume (1), the displacement constraint (17) 2 can be reduced to which states that the constraint is fulfilled, once the volume average of the fluctuations equals zero.

Algorithmic treatment of kinematic constraints
To enforce the displacement constraint in (18) on the whole RVE domain, we propose to use the method of Lagrange multipliers [2]. Similar applications can be found in [4,52]. For this purpose, the mechanical boundary value problem is recast in terms of the principle of minimum potential energy. By adding the potential λ associated with the Lagrange multipliers λ and the constraint (18) to the function of potential energy , one obtains In the following, only the terms concerning the Lagrange multiplier will be regarded, as the other 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 the Lagrange term, the additional degrees of freedom λ appear.

Variation
The potential energy is varied once with respect to the displacement fluctuations u and once with respect to the Lagrange multipliers λ, i.e.

Discretization
Using u ≈ N e d e and δ u ≈ N e δ d e as FE approximations, the discretized expressions can be written as To obtain the equivalent of the global volume integral in terms of the elements, the assembly operator A is applied for the respective matrices. For better readability, a new element matrix is defined as This simplifies the formulations to

Global matrix notation
To write the whole system of equations as a global problem, the global matrices are defined in terms of the element matrices, i.e.
where A is the aforementioned assembly operator and U a unification operator, as the node displacement fluctuations shared by different elements are not added up, but belong to the same degree of freedom. Now the expressions (25) and (26) can be reformulated in global fields as 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).

Linearization
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 form as in (28) and (29) but also their linearizations. They are used to iteratively compute the nodal displacement fluctuations as well as the Lagrange multipliers. Here the definition of the operator is used. When linearizing a function f (x) = 0 atx as Lin Applying this to the weak forms results in

Global discretized microscopic problem including constraint
From the linearized variations of λ we define the global residua Including all linearized variations of int + ext + λ yields the discrete equation as expansion of (13). After including Dirichlet boundary conditions and applying standard arguments of variational calculus, the resulting discrete system of equations reads 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, Lagrange multipliers have the disadvantage of adding new degrees of freedom to the system of equations. For the presented displacement constraint, only one extra degree of freedom is added for each spacial direction This is due to the fact that the constraint is applied on the whole RVE which avoids the approximation of the Lagrange multipliers as field variables. Thus, for three-dimensional problems, λ will only add three additional degrees of freedom. Compared to the displacement fluctuations which are linked to the nodes and which may thus easily reach thousands of degrees of freedom, the number of three additional degrees of freedom over the whole RVE is negligible, making it computationally cheap.

Coupling of the deformation gradient
The constraint related to the deformation gradient (17) 1 , can be derived and applied in the same manner as just presented for (17) 2 in the previous section. The only change in the final formulation is that the related g eT F matrix needs to be computed as the volume average of the element B-Matrix instead of the shape functions, Applying the constraint regarding the deformation gradient on the volume instead of enforcing it using periodic boundary conditions will lead to minimally invasive boundary conditions enabling e.g. arbitrary damage propagation without artificial restrictions imposed by periodic boundary conditions. As shown in [13], such minimally invasive boundary conditions result in a softer constraint compared to periodic boundary conditions. To simplify the numerical examples in this paper, only the displacement constraint is applied and the constraint related to the deformation gradient is enforced by using periodic boundary conditions.

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 appear which are consistently derived for the case where the displacement constraint proposed in the previous section is taken into account.

Macroscopic equilibrium equation
The complete macroscopic balance of linear momentum including inertia is given by Applying the Galerkin method with a test function δu on the entire domain B leads to the weak form of linear momentum B δu T Div P + f dV = 0. By applying Div(P)δu = Div(P T δu) − P : Grad δu and the Gauss theorem B Div(P T δu) dV = ∂B δu · t d A, the weak form is written as Herein, zero traction forces are taken into account at the boundary and δF = Grad δu. Analogous to the microscale, only body forces related to inertia, not gravitation, are considered at the macroscale. Thus, the body force vector is set to f :=f ρ = f ρ .

Linearization
To solve the weak form of equilibrium by using the standard Newton-Raphson scheme, the linearized balance of linear momentum is obtained as Now the operator is applied again to P and f ρ : 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 as well. In turn, the inertia forces can also depend on the deformation gradient in addition to the acceleration. We define the resulting four sensitivities as These moduli (42) are inserted into the linearized weak form which results in

FE discretization
The linearization is of the weak form of the balance of linear momentum is now discretized in terms of finite elements. First, the linear increment is discretized using standard FE formulations. Then, to get rid of the dependence on the time derivatives, the numerical time integration in terms of (10) Herein, the matrix representation of the moduli in index notation has been used. Lowercase indices refer the spacial dimension n dm , whereas uppercase indices to the total degrees of freedom of a element n edf .
By extracting the nodal virtual displacements, the element residuum is identified as where again matrix representation and index notation is used.

Consistent macroscopic tangent moduli
For the dynamic homogenization framework, four macroscopic tangent moduli (42) need to be determined. To obtain the sought-after moduli in closed form, we start with taking the derivative of the incremental linearized weak form of linear momentum at the microscale with respect to the two relevant macroscopic quantities, the deformation gradient F and the accelerationü. Then we derive the moduli by considering the microscopic problem in its equilibrium state.

Incremental weak forms including displacement constraint
As it will be shown later, the derivatives of the microscopic fluctuations with respect to the macroscopic deformation gradient and the acceleration will be required. For their calculation, the incremental linearized weak forms have to be derived with respect to these two quantities. In order to account for the proposed displacement constraint, the associated increments of the weak form of linear momentum and of the Lagrange multiplier potential with respect to the microscopic displacement fluctuations as well as the Lagrange multipliers, i.e. (20) and (21), are identified as Using the decompositions (3) and (4), equation (49) can be reformulated as

Derivatives of incremental weak forms
By taking the derivatives of the increments (50) and (51) in the equilibrium state G = 0, a closed form formulation of the tangent moduli 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 the macroscopic deformation gradient, while considering (10) results in Using standard FE discretization yields Rewriting this in global notation using the abbreviations defined in Appendix 1 Table 2 leads to the expressions 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 Then, the required derivative can be computed from Note that K * is the microscopic tangent stiffness matrix in (35), which is already available from solving the microscopic boundary value problem.
Derivative with respect toü: Analogously, the derivative of (51) with respect toü can be obtained by applying (11), i.e. one obtains Standard FE discretization using matrix representation and index notation yields Again, the two equations are combined by joining the displacements and the Lagrange multipliers into a single matrix, c.f. (58). By solving the resulting system of equations with respect to the required derivatives ∂üD * , we obtain Herein, the definitions (58)-(59), the abbreviations defined in Appendix 1 Table 2, as well as W * T = W T 0 are used.

Derivation of tangent moduli
In this subsection the four moduli will be derived by inserting the derivatives computed in the last subsection. Note that all moduli are only consistent for a microscopic equilibrium state. Thus, quadratic convergence of the macroscopic Newton-Raphson iteration is only ensured if the microscopic boundary value problem is solved for each macroscopic iteration step. After the last microscopic iteration, the consistent tangent moduli can be computed.
Derivation of A P,F : To derive the sensitivity of the macroscopic stresses with respect to the macroscopic deformation gradient, the derivative is rewritten using the definition of the macroscopic stresses in terms of the microscopic fields, i.e.
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 By using the global abbreviations defined in Appendix 1 Tables 2, 3 and inserting (61), the closed form result is obtained as Note that this result has already been presented in [57] for a special scenario of dynamic homogenization, which did not include macroscopic acceleration and the displacement constraint. . First the derivative is rewritten using the definition of the macroscopic stresses in terms of the microscopic fields as

Derivation of
Then using the chain rule, (3), (4), (10) and FE discretization, the equation reads Finally, using the global abbreviations in Appendix 1 Tables 2, 3 and inserting (66), the modulus is obtained as Derivation of A f,F : The derivation of the sensitivity of the macroscopic inertia with respect to the macroscopic deformation gradient is again similar to that of A P,F . First, the derivative is rewritten as and by using (3), (4), (10) and FE discretization, the equation reads Using the global abbreviations in Appendix 1 Tables 2, 3 and plugging in (61), the modulus is identified as Derivation of A f,u : Analogously, the derivative is rewritten Then, using (3) and FE discretization, the expression becomes It should be noted that the overall structure of the standard FE 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 be omitted from the code. It is the second derivativeü, computed in the macroscopic problem, which influences the microscopic results Taking into account the global abbreviations in Appendix 1 Tables 2, 3 and inserting (66), the modulus is derived as Note that if α 1 and α 1 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 [40]. Here, the closed form moduli (69), (72), (75) and (78) extend this consistently to the dynamic regime. Note that for a comparable approach at small strains, most derivative terms are conceptually quite similar. Also the incorporation of the new kinematic constraint in terms of Lagrange multipliers can be considered analogous for small strains. An overview over the algorithm of the proposed framework is presented in Fig. 3.

Numerical analysis: layered structure
This section presents numerical studies as a proof of concept, as well as an initial analysis of different RVE choices. As it turns out, for dynamic homogenization the definition of RVE is even more complex than for quasi-static cases. Single-scale comparisons are calculated to assess the reliability of the homogenization framework. First, a rather arbitrary example is shown to analyze the macroscopic Newton iteration and demonstrate the quadratically converging algorithm, which is based on the tangent moduli derived in Sect. 4. Then the concept of a unit cell as RVE is analyzed under dynamic conditions. Finally, a comparison of two different displacement constraints, including the proposed one, is presented.
All numerical examples make use of the Newmark scheme with the parameters γ = 0.5 and β = 0.25, resulting in an unconditionally stable algorithm. For more details on time integration methods in the context of nonlinear FE methods see e.g. [66].
A one-dimensional model of a layered structure with the total length of L is investigated. The studied heterogeneous material consists of two alternating phases, a soft and light phase, and a stiff and heavy phase. Each phase has a length of l M , a Young's modulus E 1 and E 2 , and a density ρ 1 and ρ 2 , respectively. All calculations are run using E 1 = 2·10 3 N mm 2 , E 2 = 2 · 10 5 N mm 2 , ρ 1 = 1 · 10 3 kg m 3 and ρ 2 = 1 · 10 5 kg m 3 . The Poisson's ratio is chosen to be negligible, i.e. ν = 10 −6 , to enable a quasi-1D investigation. The left boundary is fixed, on the right end an impact load is applied in terms of a displacement boundary condition using the polynomial function u(t) = 2 8 4 , where u max is the amplitude of the impact wave and T the duration in 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 2 framework, c.f. Fig. 4 respectively. The DNS discretizes the microscopic phases at the macroscale into a large number of finite elements with a length of l E . It thereby serves as overkill reference for the multiscale approach. The FE 2 simulations have a macroscopic element length of l E and make use of the same element length l E at the microscale for better direct comparability of the microscopic fields to the DNS. To approximate the displacement fields of the elements, linear shape functions and two Gauss points are used for all scales. As shown in Fig. 4b, each microscopic RVE calculation is associated to a single macroscopic integration point. The corresponding parameters for each simulation, regarding geometry, material parameters and loading will be listed in the caption of each figure.

Consistency of numerical framework
This example analyzes the convergence behavior of the macroscopic Newton iteration. In Fig. 5a, the distribution of macroscopic displacement fields is shown at three different time instances for both, the DNS and the FE 2 calculation. As RVE, the basic unit cell of the type A (c.f. Fig. 6) is used. It can be seen, that the dynamic multiscale framework approximates the overall 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, but this would generally make the calculation converge faster as the initial values are already closer to the solution, defying the objective to properly test the tangent moduli. The convergence behavior for the three arbitrarily chosen time frames is depicted in Fig. 5b. Quadratic convergence of the norm of the updates of the nodal displacements D * is observed. This demonstrates that the macroscopic tangent moduli, incorporating both the microscale inertia forces as well as possible constraints, have been derived in a consistent manner. Note that in this example, locally strains appear of over 10% which exceed the range of small strains. However, due to the onedimensionality of the problem, no significant influence from the finite strain setting can be expected.

Analysis of the unit cell concept under dynamic loading
For quasi-static homogenization simulations of periodic microstructures, it is known that the resulting macroscopic answer, as well as the corresponding microscopic fields are invariant with respect 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 center matters in a dynamic setting. An extreme example is shown in Fig. 7, which compares the macroscopic displacement field at t = 0.045 s presented in the first example in Fig. 5a with a simulation using the basic unit cell type B as an RVE (c.f. Fig. 6). To properly measure the influence of different RVE choices on the FE 2 simulation, an objective error measure is considered. It is defined as average difference of the macroscopic displacement fields = n nodes i u I i (t j ) − u II i (t j ) /n nodes , where u I and u II stand for any two displacement signals that are being compared. This measure can be evaluated for each time step and thus, the average is once more computed over the number of time steps time = n timesteps j j /n timesteps . Note, that due to the comparison of nodal displacements, a unitless error value obtained by dividing the difference by the reference value is problematic as the regarded value might be zero. Furthermore, a delay in the response may lead to relatively large errors in the displacements summarizing over time even if the shape of the displacement wave is perfectly fine. Therefore, in principle it is difficult for dynamic problems to obtain a quantitative measure to evaluate the accuracy. The interpretation of the error value should therefore be done with caution, as the single value does not convey any information about the actual shape of the wave. However for the present analysis the considered error measure is sufficient to analyze the general trend. Figure 8 shows the calculated errors of different choices of the RVE. For the comparison, RVEs with multiple periods of the same unit cell type, as depicted in Fig. 6, were considered. Two effects can be observed: The first, presented in Fig. 8a, is that the difference in macroscopic displacements between different choices of unit cell type decreases, when the number of unit cells per RVE is increased. This means that the choice of particular basic unit cell type does not matter as long as the RVE is chosen large enough. The second effect, shown in Fig. 8b, is that 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 violation of 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 practical application when using FE 2 . At this point it is favorable to use domain decomposition approaches instead of a homogenization method in order to avoid the scale separation assumption.

Influence of displacement constraints
Finally the proposed displacement link u = u is analyzed for the examples in Fig. 8b, in comparison to the standard Different choices of RVEs and constraints were analyzed displacement link for quasi-static periodic homogenization, where the fluctuations at the RVE corner nodes is set to zero, i.e. u corner = 0. For the quasi-1D example analyzed here, this is equivalent to setting the integral over the surface equal to the corresponding macroscopic displacements, which has been taken into account in other dynamic homogenization schemes.
The first observation is, that using the proposed volume constraint u = u results in a more robust framework in terms of stability of the Newton-Raphson iterations. Furthermore, slightly smaller error values are obtained compared to the DNS reference. Table 1 shows the number of time steps reached before either the calculations crashed (due to diverging Newton iterations at the microscale) or they were finished successfully after the intended complete set of 1000 time steps. Especially the calculations using the unit cell type B in combination with the zero fluctuations of the corner nodes, underperformed the other scenarios. To understand the difference between the performance of the displacement links, it is necessary to examine the behavior at the RVE level.
Here the examples with the RVEs consisting of three periods of the basic unit cell B are further analyzed. Figure 9 compares the microscopic displacements for four relevant time instances right before the peak of the input wave passes through the RVEs. More specifically, the differences between the microscopic displacement fields u of an RVE and the respective macroscopic displacements u. To compare the DNS, an effective u has been computed as the average displacement over the associated length. Thereby, the quality of the microscopic fields can be analyzed independently of the macroscopic displacements. With this, the two different displacement constraint options can be effectively compared with the reference solution obtained from DNS. The graphs show, that the fixed corner constraint leads to artificially increased displacement intensities at the microscale due to the constricted boundary. These increased displacements eventually lead to extreme deformations in single elements at the microscale, crashing the simulation. The proposed displacement volume constraint however, leads to a softer constraint which results in a more robust computation while still enabling dynamic effects which agree well with the ones from the reference DNS. In the presented examples, the only rate-dependent influence are inertia forces. In cases where also rate-dependent material properties are included we expect the influence of different displacement constraints on the overall simulation to increase, in favor for the proposed volume constraint.

Numerical analysis: split Hopkinson bar
After presenting an analysis focused on the study of general properties of dynamic homogenization, this section shows the applicability of the framework to an engineering problem. A standard experiment to study material behavior under impact loading is therefore replicated, the split Hopkinson tension bar test, see e.g. [11,58]. By releasing a pre-strained steel bar, a loading pulse is transmitted to an aluminum input bar. The wave travels along the bar and through a specimen, which is sandwiched between the input and another aluminum bar called the output bar. Strain gauges in the two bars record the wave signal. By applying the elastic, uniaxial stress wave propagation theory to the split Hopkinson bar experiment, the time history of the forces and the displacements of the faces of the test specimen are calculated, in the following denoted as σ 1 and σ 2 , cf. [28]. These are then used to approximate the stress and strain within the specimen in loading direction. A schematic visualization is presented in Fig. 10. As target material, a strain-hardening cementitious composite (SHCC) is chosen. This fiber-reinforced concrete exhibits an outstanding ductility and a pronounced energy dissipation under high strain rates. Therefore it is ideal as reinforcing layer to improve the impact resistance of structures [12]. The interpretation of experimentally measured data under dynamic conditions is difficult. Therefore there is a need for accurate simulations including inertia at the macroand the microscale. In [58] a quite simplified microstructure was considered as RVE. Although a three-dimensional discretization of the RVE was taken into account, the multidimensional character with respect to mechanics was strongly limited by the fact that only one fiber in the direction of the main macroscopic traveling wave direction was taken into account. Therefore, in contrast to [58], here we consider an RVE with multiple fibers to more realistically reflect the three-dimensional nature of microscopic problems in real SHCC. In total, 13 randomly oriented fibers are included representing a more or less isotropic distribution of fiber orientations. First the applied micromechanical models and the chosen microstructure are presented. The macroscopic boundary value problem is shown and used to analyze the dynamic effects.

Microscale problem
The microscale problem discretizes the fiber-reinforced concrete. Therefore the relevant micromechanical features need to be replicated. Before the first crack through the cementitious matrix, the concrete matrix itself dominates the material behavior of the composite. Once a crack has formed the fibers are engaged and bridge the crack. Finally, a crack will open when the fibers are pulled out. The fiber properties as well as the pullout behavior are shown to be rate dependent [9,11]. To capture these micromechanical effects two main models are required: (i) a model for the concrete matrix, including cracking and (ii) a representation for the embedded fibers. The concrete is represented as a first approximation by a standard Neo-Hooke material law. The complex compressive behavior of concrete is neglected as only tensile loading is considered within this example. A realistic representation of the crack development is not within the scope of this work. Therefore, to represent the matrix cracking, a simple erosion method is implemented. It sets the material stiffness of an integration point close to zero, once the stress in loading direction surpasses a critical value σ cr . The fibers are represented by a linear truss element, sharing the nodes with the matrix mesh. This allows for a 1D effective material law to be applied, capturing the complex material behavior of the fiber as well as that of the fiber-matrix bond. However, the forces are thus not transferred to the matrix continuously, leading to stress concentrations at the nodes resulting in spurious crack patterns. Therefore, the crack path in the matrix is here defined in advance. The chosen material model uses a general 1D Neo-Hookean material law with an additional strain rate sensitivity of the stress as well as a damage formu- lation, cf. [58]. The stress in terms of the first Piola-Kichhoff stress tensor is given by where E is the Young's modulus, denotes the dynamic increase which takes on only positive values, and D is the damage value which takes on values from 0 to 1. The dynamic increase function is defined by a logarithmic function, which depends on the rate of the deformation gradientḞ. For values ofḞ ≥ α II , it is defined as = α I ln[Ḟ α II ] otherwise it is zero, i.e. no dynamic increase for rates lower than α II . The exponential damage formulation is given by The damage value D is determined by the effective energy considered for damage ψ D . It is defined as the current maximum value of the strain energy function ψ 0 . The damage algorithm is therefore classified as a discontinuous damage approach. The damage evolution is controlled by three parameters. D ∞ defines the maximum damage reached, D rate > 0 sets the velocity of the damage evolution, and D shape modifies the shape of the function.
As RVE, a cubic body with an edge length of 1 mm is chosen, see Fig. 10a. A crack face perpendicular to the loading direction is located in the center, splitting the matrix in two halves, which are connected by a number of 13 randomly oriented fibers. The concrete matrix is modeled by brick elements, whereas in the center, the feature for cracking is included to represent a possible crack through the matrix. This microstructure fulfills minimal geometric requirements to reproduce the relevant micromechanical processes of SHCC. More complex structures are only of value once advanced micromechanical models are utilized. In the subsequent simulations the following material parameters are used: for the concrete matrix E = 29 kN/mm 2 , ν = 0.3, and ρ 0 = 2100 kg/m 3 , additionally for the crack E cr = 10 −3 kN/mm 2 and σ cr = 5 kN/mm 2 , for the fiber model E = 40 kN/mm 2 , A = 0.00385 mm 2 , ρ 0 = 980 kg/m 3 , D ∞ = 0.9982, D shape = 0.36, D rate = 0.2, α I = 0.08, and α II = 0.51.

Macroscopic boundary value problem
The macroscale represents the experimental setup of the split Hopkinson tension bar. The cylindrical equipment and SHCC sample are discretized by truss elements. Standard elements are used to simulate the aluminum input and output bars, the two-scale homogenization framework is applied to model only the SHCC specimen. A sketch of the setup is shown in Fig. 10b. All elements have a length of 10 mm. To simulate the experimental loading conditions a piece-wise polynomial function u BC is chosen. The load is applied as a boundary displacement. The three parts are defined as u . The transitions between the respective functions are at u I (0.592 t vc ) = u II (0.592 t vc ) and u II (t vc ) = u III (t vc ).
The function describes the pulse, which is characterized by two phases. First the acceleration phase and second a phase of constant velocity. The two parameters of the function t vc and v c , define respectively the time when the transition from the first to the second phase is reached and the constant velocity.

Results
By using the presented microstructure in combination with the split Hopkinson tension bar simulation, the capabilities of the developed framework to analyze material behavior of composites under dynamic loading are illustrated. Following the experimental procedure, the stress signals σ 1 and σ 2 of the nodes at the specimen interfaces is averaged as σ . The stress is then plotted against an approximated specimen strain, computed by the difference of displacement of the interfaces divided by the specimen length. First a quasi-static simulation is compared to the dynamic response, shown in Fig. 11. A clear increase in the dynamic simulation is observed. In addition, a shift is visible from the sudden drop of stress under quasi static loading, due to the homogeneous stress state, to a more gradual stress increase under dynamic conditions, mimicking the multiple cracking behavior of SHCC. Then, to understand the origin of the dynamic increase, two parameter studies are conducted. In Fig. 12, the strain rate sensitivity of the fibers is varied. Clearly, the overall stress response is increased by an increase in α I . In Fig. 11 Quasi-static and dynamic results of the split Hopkinson tension bar simulation. The zoom shows the initial cracking of the matrix Fig. 12 Variation of the strain rate sensitivity α I Fig. 13 Influence of microinertia on the macroscopic response addition, the unsteady phase of multiple cracking is completed at smaller overall strains, due to the effective increase in fiber stiffness. A similar effect is observed when comparing microstructures with random fiber direction to perpendicular fibers in loading direction. The fibers at an angle to the applied load result in the softer macroscopic response, c.f. [58]. In Fig. 13, the focus is on the influence of the microinertia. It is evident that the main dynamic influences of this example are the structural inertia and the strain rate sensitivity of the fiber pullout. The microscopic inertia, in contrast to the earlier example of the layered structure in Sect. 5, does not appear to have a significant influence on the macroscopic stress. As the chosen microstructure only allows for moderate dynamic activity, this is not surprising. However with more advanced micromechanical models this might change.

Conclusion
In this paper, a general purpose, consistent, two-scale homogenization framework for dynamics at the macro-and microscale was proposed in the sense of the FE 2 method. The novelty of this framework lays in its generality. The framework does not include any simplifications such as linearized strains, explicit time integration or partly quasistatic scenarios. The only assumption taken into account is a sufficiently pronounced scale separation, which is anyway essential requirement of any FE 2 approach. Therefore, it enables the simulation of various structural problems of complex micro-heterogeneous materials under dynamic loading such as impact. Furthermore, the derived formulations are compatible with standard FE 2 architecture. The main aspects are: (i) the incorporation of the complete balance of momentum at the micro-and macroscale including inertia forces, (ii) extended microscopic boundary conditions resulting in a more robust scheme and giving the possibility for nonperiodic RVE boundaries, (iii) a finite-strain formulation enabling the simulation of a wide range of macro-and micromechanical phenomena, and (iv) the derivation of consistent macroscopic tangent moduli ensuring quadratically converging macroscopic iterations. The presented results on different choices of RVEs show that different admissible unit cells of a periodic microstructure lead 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 the most basic unit cell is necessarily a bad choice for a simulation, but the specific choice of this basic unit cell is not unique. In addition, a more complex simulation of a split Hopkinson tension test on an SHCC specimen was presented. The example at hand, compared to the simplified microstructure in [58], exhibits a softer macroscopic response, as not all fibers are oriented in loading direction. This is apparent due to the more pronounced phase of multiple cracking up to 1.5% strain. The example shows the capability of the multiscale framework to further the understanding of experimental material behavior under impact.
Acknowledgements The authors gratefully acknowledge financial support from the German Research Foundation (DFG) within project B1 of the Training Research Group (GRK) 2250 "Mineral-Bonded Composites for Enhanced Structural Impact Safety". Furthermore, fruitful discussions with Celia Reina (University of Pennsylvania) are greatly appreciated.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix: Matrix abbreviations
See Tables 2 and 3.