A TWO-SCALE RATE DEPENDENT CRACK MODEL FOR QUASI-BRITTLE MATERIALS UNDER DYNAMIC LOADING

A multi-scale numerical approach for modeling cracking in heterogeneous quasibrittle materials under dynamic loading is presented. In the model, a discontinuous crack model is used at macro-scale to simulate fracture and a gradient enhanced damage model has been used at meso-scale to simulate diffuse damage. The traction-separation law for the cohesive zone model at macro-scale is obtained from the meso-scale through the discontinuous computational homogenization method [1] which is developed based on the so-called failure zone averaging scheme [2] in which the averaging theorem is used over the active damage zone of the meso-scale. Unlike standard averaging, this method is objective with respect to the local-scale sample size in the softening regime. In order to evaluate the macroscopic traction at each integration point on the crack, at each time step of the macro model solution, a static boundary value problem is solved for the representative volume element (RVE) whose size is significantly smaller than the macro length-scale and the macroscopic wave-length. The effect of the crack opening rate on the macro cohesive law is taken into account by relating the material properties of the meso-scale model to the macro crack opening rate. The objectivity of the model response with respect to the representative volume element size is demonstrated for wave propagation problems. The model is verified by comparison with a direct numerical simulation (DNS).


INTRODUCTION
Optimal design of quasi-brittle materials is of great interest nowadays for civil and defense structures.Many of the materials which are used in such engineering structures are heterogeneous and have more than one length scale.Modeling heterogeneous materials using a direct numerical simulation (DNS) in which detailed heterogeneities are modeled directly in the macro-scale can give accurate results but this method needs enormous computational efforts and is in most cases not practical.
Homogenization-based multi-scale methods are widely used to obtain the macroscopic behavior of the heterogeneous materials by averaging local-scale properties.The multi-scale methods can be analytical or computational.Analytical methods are limited to simple problems and cannot be applied to more complex structures.Therefore, numerical [3] and computational homogenization [4,5] methods have received significant attention.In computational homogenization methods, heterogeneous material is replaced by a homogeneous one with unknown macroscopic constitutive behavior.Then, a representative volume element (RVE) is associated to each material point and the constitutive law is obtained by solving a boundary value problem for the RVE.
Existence of the RVE is a crucial issue in such techniques.A sample volume can be defined as RVE when homogenized properties do not change significantly with varying RVE size.When a standard homogenization method is used, in linear and hardening regimes, the RVE can be defined but in the softening regime a RVE cannot be defined [6].The existence of a RVE for quasi-brittle materials with random complex heterogeneous structures is shown in [2], where a failure zone averaging scheme is introduced.In [1], based on the failure zone averaging scheme, a homogenization model for modeling cohesive crack in quasi-brittle materials has been developed.Using this scheme, macroscopic cohesive laws which are independent of RVE size can be obtained from the local-scale with localized deformation.
Wave propagation problems in heterogeneous materials using the multi-scale method are studied by many researchers.A dispersive wave model for wave propagation in heterogeneous materials is developed in [7] based on higher order mathematical homogenization with multiple spatial and temporal scales.In [8] a multi-scale model for composite materials under impact loading is developed in which damage is introduced in the model using locally homogeneously distributed microcracks through a cohesive zone model.In their multi-scale model, the macroscopic problem is solved dynamically while the local scale is considered a static problem.
In the present work, the proposed homogenization scheme by Nguyen et.al. [1] is extended to wave propagation problems.Furthermore, rate effects are added to the model by relating the material properties of the RVE to the rate of the macroscopic crack opening.Numerical results are given to show the capability of the model.

COMPUTATIONAL HOMOGENIZATION APPROACH
In the multi-scale method, local scale averaged properties are used to obtain macroscopic behavior of the heterogeneous material.Figure 1 shows a schematic description of the classical computational homogenization method.As shown in this figure, the heterogeneous material is substituted with a homogeneous one.For each point in the macro model with strain M , the strain is imposed as a boundary condition on the external boundary of the RVE associated to this point.After solving the boundary value problem for the RVE, the corresponding macroscopic stress σ M and tangent moduli C M are obtained.
To use computational homogenization theory, the problem must meet the following requirements.Firstly, the RVE should exist for the heterogeneous material.The RVE exists if an increase in size does not change homogenized properties and the sample is large enough so that the meso-/micro-structure randomness does not affect the homogenized properties.The second important issue in computational homogenization is the principle of seperation of scales which indicates that the macroscopic characteristic length scale, l M , is assumed to be much larger than the local-scale length, l m (see figure 1).
The third requirement is the strain averaging theorem which states that the macroscopic strain (stress) at any point is defined as the volume average of strain (stress) over the RVE associated with that point and mathematically can be written as: where M (σ M ), m (σ m ) and Ω m are macroscopic strain (stress), local-scale strain (stress) and RVE volume, respectively.Energy consistency in transition of scales is satisfied by the Hill-Mandel macro-homogeneity principle which states that the macroscopic work rate must be equal to the volume average of local-scale work rate over the RVE, according to: The macro-scale mass density can be defined as: where ρ M and ρ m are macro-scale and meso-scale mass densities, respectively.
In wave propagation problems, if the length of the propagating wave at macro-scale is significantly larger than the meso-scale length, l m , one can neglect wave effects at mesoscale and it is possible to simplify the meso-scale problem to a quasi-static one.By using this assumption, it is not possible to model wave dispersion effects.But, in case of low frequencies in which this simplification will be used, dispersion effects are insignificant.

DISCONTINUOUS COMPUTATIONAL HOMOGENIZATION SCHEME
When strain localization occurs in the material, strain and stress homogeneity over the meso-scale domain is not valid and in such cases the multi-scale solution depends on the size of the fine scale which means that a RVE cannot be defined.Nguyen et.al. [1] developed an objective homogenization scheme for such strain localization problem in which stress/strain averaging and the Hill-Mandel theories are performed over the localization band instead of the whole domain of the fine scale.In the discontinuous homogenization scheme, at macroscale the crack is exhibited using a cohesive zone model in a homogeneous solid.Then, the traction-separation law for the crack is obtained from the meso-scale model.

Macro-scale model
The macro-scale model is illustrated in figure 2. Macro cracking is modeled using the XFEM.In the finite element model, the momentum equation can be written as: where üM represents the macroscopic acceleration vector, M is the mass matrix.f ext M is the external force vector, f Bulk M and f Coh M represent the bulk force vector and the cohesive force vector, respectively and are given as: in which t M is the macro-scale traction and N and B are the matrix of nodal shape functions and the matrix of derivatives of the shape functions, respectively.The bulk macro-stress can be computed as: The fourth-order tensor D 0 is the bulk homogenized tensor which can be computed using a standard homogenization technique.The macro traction, t M , is obtained from the cohesive law via: where [[u]] M is the displacement jump for the macro crack and T M is the macro cohesive tangent.At each time step, the displacement jump is obtained for each integration point on the crack and the corresponding macro traction, t M , and macro cohesive tangent, T M , are computed using the discontinuous homogenization scheme from the meso-scale model.

Meso-scale model
At meso-scale, failure is modeled using the implicit gradient-enhanced damage model [9].The stress-strain relation is given as: where ω is the scalar damage variable (0 ≤ ω ≤ 1) and D is a fourth-order tensor which contains the elastic moduli.The damage evolution law is written as: where γ, β and κ I denote residual stress, softening slope and damage threshold, respectively.κ is a scalar measure of the largest strain ever reached and is defined by loading function f as: f and κ satisfy the KuhnTucker conditions: ¯ eq is the nonlocal equivalent strain which is implicitly related to the local one according to: In this equation, c is defined as c = 1 2 l 2 c and l c represents the length scale.The local equivalent strain is defined as: where i are the principle strains and x means the positive part of x.The discrete system of equations for meso-scale model can be written as: The internal force vector for the meso-scale model can be obtained from f int m = Ωm B T σ m dΩ and f ext m is a function of the macroscopic displacement jump.For time step t i , at each iteration of the macro-scale model solution, equation ( 14) is solved for the meso-scale model.

Macro-meso transition
Figure 3 depicts the discontinuous homogenization scheme.The connection between macro-scale and meso-scale is given by: where u R is the total displacement at the right edge of the RVE.The first term in the RHS represents the linear displacement and u 0 dam is the compatibility displacement.w and l denote the width of the RVE and the averaged width of the localization band, respectively (figure 3).Matrix C 0 is the projection of compliance tensor D −1 0 on the crack plane and is obtained as: The failure zone averaging scheme is used to compute averaged quantities for the meso-scale model.It should be noted that in this scheme, the averaged quantities are calculated over the active damaged zone which contains integration points which are damaged and are loading.The active damage zone, Ω d , can be expressed mathematically as The meso-scale quantities can be defined through: where | • | represents the area of the domain and h is the height of the RVE.u 0 dam is calculated at the moment of crack initiation using above equations.By solving the system of equations ( 14) and (15), one can find the macroscopic traction, t M , and cohesive tangent, T M .More details on theoretical and computational aspects can be found in [1,10].

NUMERICAL RESULTS
To study wave propagation in strain localization problems, a heterogeneous beam is subjected to a constant velocity at both ends (figure 4).Tensile waves propagate through the beam and after superposition of the waves at the center of the beam, the stress at this point exceeds the tensile strength and a crack initiates.Figure 5 shows the multi-scale model of the problem.Voided structures with different sizes are chosen as RVE for this problem.It should be mentioned that the multi-scale scheme is applied only on the crack and the bulk part is solved using the standard finite element method.The material properties for the RVE and the bulk material are given in table 1.A constant velocity equal to 0.3 (m/s) is applied at both ends of the beam.Cohesive laws computed from different RVE sizes, according to the failure zone averaging scheme, are illustrated in figure 6.It can be observed that the results are objective with respect to RVE size.In order to verify the multi-scale model, the results are compared with a DNS model.Figure 7 depicts the DNS model in which the material properties of the voided part and bulk part are similar to those of the RVE and the bulk part of the multi-scale model.Averaged stress over active damage zone (similar to the averaged strain in equation ( 17)) versus damage opening, u dam , for the DNS model and the multi-scale model are shown in figure 8 which shows good agreement.The difference between the results in elastic branch is due to the fact that the mesostructure is not present in the multi-scale model before crack initiation and we do not use averaged properties for the bulk part before crack presence.
The RVE failure mode using the multi-scale model is also compared to that of the DNS model at time step t=0.0128e (ms) in figure 9.This comparison also demonstrates that the development of the damage zone for both models is similar.

RATE-DEPENDENT COHESIVE LAW
There are two sources for rate dependency in concrete materials [11]: (i) the viscoelasticity of the material behavior in the bulk of the structure, (ii) the rate process of the bonds breakage in the fracture process zone.Both mechanisms are important for concrete but in high strain rate dynamic loading, the latter is the dominant mechanism which causes the cohesive law to be rate dependent.Bažant [11,12], by considering fracture as a thermally activated   phenomenon, derived a rate-dependent softening law.Here, we consider mode I and for the traction in normal direction x to the crack surface, the rate dependent softening law can be written as: where x M denotes the macro crack opening rate and t 0x M is the traction under static loading condition.c 0 and c 1 are material parameters.Here, we assume that, when a crack initiates, the strain threshold, κ I , in the gradient damage model which is used for meso-scale model, is dependent on the crack opening rate through: in which κ 0 I is the static strain threshold.In order to investigate this assumption, cohesive laws are computed for various values of κ I which are obtained from equation ( 19

) for [[u]]
x M = 0.0, 0.25, 0.5, 1.0 (m/s).Here, c 0 and c 1 are taken equal to 0.8 and 0.5, respectively.In figure 10 .From figure 10, it can be concluded that: The above relation shows that equations ( 18) and ( 19) are almost equivalent.So, in order to capture rate dependency effects in the macro-scale cohesive law, one can insert rate effects in the meso-scale model using equation ( 19).In the solution procedure, at time step t i , for a certain crack in the macro-scale model, the crack opening rate is calculated and then the strain threshold for the RVE corresponding to the integration points on this crack is updated using equation (19).To obtain a more accurate result, the problem is solved again for time step t i with updated values for strain threshold, κ I .
The problem described in figure 4 is now considered for a crack with a rate-dependent cohesive law.The multi-scale problem is solved for different loading rates.Figure 11 illustrates the computed cohesive laws for various RVE sizes at different loading rates.As it can be observed in this figure, for a given crack opening, the traction increases with loading rate.It is also obvious that the obtained softening laws are objective with respect to the RVE size.
In order to verify the model, a DNS model is presented as before.In the DNS model, the relative velocity values between RHS and LHS of the voided part (parts shown with red lines in figure 7), after damage initiation, is taken as the crack opening rate.A comparison of crack opening rate in multi-scale model and DNS model for V 0 =0.3 (m/s) is shown in figure 12.The averaged stress over the active damage zone versus damage opening is given in figure 13 for the multi-scale and the DNS model at various loading rates.It can be observed that for lower velocities the results are in good agreement.But, at higher loading rates, the curve obtained for the DNS model is above the multi-scale curve and the difference between these two curves increases with loading rate.This difference stems from the inertia forces around the damaged zone in the DNS model.In the multi-scale model, as discussed before, the inertia forces at meso-scale are neglected in the present work and as a result the effects of inertia forces cannot be captured.Nevertheless, even at high rates, the multi-scale model is capable of properly calculating the material response.To illustrate this fact, the density of the voided part in the DNS model is assumed to be very small so that the inertia forces around the damaged zone are negligible.Averaged stress-damage opening curves are shown for V 0 =1.0 (m/s) in figure 14.It can be observed that the curves for the DNS model and the multi scale model lie on top of each other when the inertia forces are neglected in the voided part.Figure 12.Crack opening rate at V 0 =0.3 (m/s) for DNS and multi-scale models.

CONCLUSION
For the modeling of wave propagation in heterogeneous materials using multi-scale methods, one can perform a quasi-static analysis at lower-scale as long as the length of the macroscopic propagating wave is significantly larger than the lower-scale length.In case of strain localization problems where standard homogenization methods cannot be used due to the dependence of the results on RVE size, the discontinious homogenization scheme, based on a failure zone averaging technique, is applied to obtain objective results with respect to the RVE size.Good agreement between results obtained from the multi-scale model and the DNS model, not only certifies the capability of the discontinuous homogenization method but also supports the idea of neglecting wave propagation at the lower-scale problem.
Rate effects due to bond breakage in the fracture process zone can also be modeled in the multi-scale scheme by using a rate dependent material model at meso-scale.

Figure 5 .
Figure 5. Multi-scale model and different RVE sizes.

Figure 6 .
Figure 6.Computed cohesive law for different RVE sizes.

Figure 8 .
Figure 8.Comparison of averaged stress over failure zone vs. damage opening for multi-scale model and DNS model.

Figure 9 .
Figure 9.Comparison of failure modes for multi-scale model and DNS model.

Figure 11 .
Figure 11.Computed cohesive laws for different RVE size at various loading rates.

Figure 13 .
Figure 13.Averaged stress over the active damage zone versus damage opening for DNS model and multi-scale model at various loading rates.

Table 1 .
Material properties for bulk material and RVE.