Modelling of High Velocity Impact on Concrete Structures Using a Rate-Dependent Plastic-Damage Microplane Approach at Finite Strains

Concrete is known as a quasi-brittle material and the microplane model has been proven to be a powerful method to describe its constitutive features. For some dynamic cases, however, numerous microplane models used successfully at small strains are not sufficient to predict the nonlinear behaviour of damaged concrete due to large deformations. In this contribution at hand, a combined plasticity-damage microplane model extended to the finite strain framework is formulated and regularised using implicit gradient enhancement to achieve mesh insensitivity and to obtain more stable finite element solutions. A modified smooth three surface Drucker–Prager yield function with caps is introduced within the compression-tension split. Moreover, a viscoplastic consistency formulation is implemented to deliver rate dependency at dynamic cases. In case of penetration into concrete materials, the proposed model is equipped with an element erosion procedure to yield a better approximation of crack patterns. Numerical examples on impact cases are performed to challenge the capability of the newly proposed model to existing experimental data.


Introduction
Throughout a building's service period, complex behaviour can be observed in concrete structures, such as for instance military base anchorages, nuclear reactors, long-span bridges, and oil platforms, which are subjected to impact phenomena, natural disasters such as earthquake or seismic loading, detonation, or blast loading with high intensity. These circumstances lead to dynamic fractures in concrete material up to the final failure of its structure. Finite element simulations are capable of modelling various phenomena in concrete using material modelling developed extensively up to now. Concrete itself is known as a quasi-brittle material with heterogeneous complex nature as a result of its constituent materials, such as cement paste and different type of aggregates. The behaviour of concrete changes from linear elastic to highly nonlinear inelastic due to several physical phenomena as microcracking, crushing or even plasticity caused by tension, compression, or combination of both tension and compression. These conditions yield difficulties in modelling concrete, especially in terms of constitutive laws and numerical aspects regarding the finite element method.
Several approaches for modelling these phenomena have been employed in order to describe concrete responses realistically in the context of continuum mechanics. The most renowned approaches are using continuum damage mechanics, plasticity, or a combination of damage and plasticity. Continuum damage mechanics is characterised by the degradation of the material stiffness describing the inelastic behaviour at post-peak range. Meanwhile, plasticity deals with elastic-plastic conditions as a result of the stress state, which are evaluated with respect to the applied yield criterion. At monotonic loadings, in fact, either damage or plasticity models are sufficient to represent softening behaviour of concrete. Nonetheless, by the fact that a structure is subjected not only to monotonic loading but also to unloading-reloading states altogether, the necessity of a coupled plastic-damage formulation is inevitable. Various plastic-damage approaches for concrete can be found in [1][2][3][4][5][6][7][8][9]. Moreover, other plastic-damage models with different features are also available, e.g., for impact failure [10], for concrete spalling considering creep-shrinkage [11], for ductile failure [12], and for granular materials at low confining pressure [13].
In the last three decades, a powerful model for concrete materials, namely the microplane approach, has been studied and developed extensively since it was pioneered by [14]. The concept of microplanes is characterised by a projection of the strain tensor to random planes of a microscale. This concept is straightforward to understand the relation between strain and stress vectors in each microplane, which is eventually equivalent to macroscopic strains and stresses. This approach was then improved by [15], namely the model M1 and subsequently, it was developed further up to the recent version M7 employing stress-strain bounds [16]. As discussed in [17], however, the conventional microplane models using either the normal-tangential (N-T) or the volumetric-deviatoric-tangential (V-D-T) split have several mechanical deficits. The full range of Poisson's ratio cannot be represented by the N-T split, while microplane elastic constants with the unique macroscopic-microplane relation is not possible to be considered by the V-D-T split. As observed in [18], the ambiguity of the deviatoric-tangential projection allows a lot of freedom in the V-D-T split and it leads to problems in numerical implementation. Therefore, the deviatoric and tangential components in the V-D-T split were merged and resulted in the volumetric-deviatoric (V-D) split, see [19].
For strain softening at the post-peak regime, however, problems relating to ill-posed differential equations and strain localisation occur. In consequence, these issues lead to pathological mesh dependencies and unstable finite element solutions eventually indicated by slow convergence rates. The so-called implicit gradient enhancement introduced in [20,21] has been proven to tackle these problems, and it has been successfully implemented to regularise the microplane damage model using the V-D split in [22]. Moreover, the gradient-enhanced microplane model was combined with Drucker-Prager plasticity in order to show the softening response of concrete [23], while the coupled plastic-damage model was developed in [24] for modelling concrete in the state of unloading and reloading.
Notwithstanding that the established microplane approaches at small strains have been effectively used to describe concrete behaviour, but finite deformations occur in concrete at dynamic cases. In particular, at high confined pressure on tube-squash tests according to experiments conducted in [25], approximately 30-50% of concrete strains were exhibited at deformed shapes with no visible damage or cracks observed. This phenomenon has been supported by numerical simulations using a finite strain microplane damage model in [26], while simulation results obtained by the small strain version showed that the concrete lost almost its entire stiffness. Meanwhile, using softening microplane plasticity at finite strains proposed in [27,28], large plastic deformations were demonstrated for similar cases. Furthermore, strains up to 100% may happen in the state of impact loading due to bored injection piles or anchors, whereas strains achieve approximately 30% while subjected to blast loading or ground shock waves, as mentioned in [25]. Meanwhile, under seismic or earthquake loading, columns of structures were deformed at around 15% of strains until collapse based on [25]. In the case of projectile or missile penetration into a concrete wall, beyond 100% of effective plastic strain occurred [29,30].
Due to all aforementioned circumstances, the small strain microplane model should be extended to the finite strain framework. Hence, the plastic-damage microplane formulation at finite strains regularised by the implicit gradient-enhanced method is proposed within the V-D split. The described model is also equipped with rate dependency to deliver strain rate effects at dynamic cases. Moreover, a coupled contact-adaptive element erosion is implemented herein for the penetration situation. The objective of the present work is to introduce and formulate the rate-dependent plastis-damage microplane model for concrete at finite strains in order to achieve more accurate and reliable responses, particularly at high velocity impact simulations. This work is organised as follows. First, the research motivation is summarised and followed by the proposed constitutive description as well as algorithmic features implemented. Subsequently, numerical examples of impact load cases are simulated and, finally, the results are compared to existing experimental investigations in order to challenge the capability of the newly proposed formulation.

Microplane Approach with V-D Split
The present work continues to use the V-D split for extending the small strain microplane model to the finite strain regime employing the conjugate strain-stress pair, namely the Green-Lagrange strain tensor and the second Piola-Kirchhoff stress tensor, which is transformed to the Cauchy stress tensor at post-processing [26,28]. Both predecessor models showed that large deformations occur in concrete either at the damaged region or at the plasticity domain. Due to these phenomena, this subsequent work will combine both damage and plasticity within the microplane approach at the finite strain framework using the proposed model in [26,28]. As discussed in [29,31], the Green-Lagrange strain tensor E is chosen to extend the small strain microplane description for concrete to the finite strain formulation computed as with the right Cauchy-Green tensor C obtained by the deformation gradient F as follows where X and x are the initial and spatial Cartesian coordinates, respectively. Rather than the Biot or Hencky strain tensor, the Green-Lagrange strain tensor is selected since it is the only admissible strain tensor for extending the small strain microplane model to the finite strain formulation due to several reasons. A direct physical interpretation, such as the V-D or N-T decompositions, is clearly explained by the Green-Lagrange strain tensor. Stretch characterisations in the normal strain component should be independent of the stretch in other directions. Similarly, the representation of shear angles in the shear strain component should also be independent of the shear angle in other microplanes as well as the stretch in other directions. Solely the Green-Lagrange strain tensor satisfies these circumstances and, hence, it is chosen to be the strain measure in this extended model.
In general, the V-D decomposition at finite strain theories is multiplicative. Nevertheless, by the fact that the volumetric strain of concrete is always small at all conditions, the multiplicative formulation could be simplified additively for the extension using the Green-Lagrange strain tensor [29,31]. Please note that for other materials with large volume changes, the multiplicative decomposition should be employed instead of the additive split. According to the V-D-T split, the microplane normal strain E N is related to its volumetric and deviatoric parts, E V and E ND , as follows where E N and E V are computed as While E is obtained from Equation (1) and n denotes the normal vector, E 0 is described by the Biot strain As the kinematic constraint of the V-D split is used here, the deviatoric component E ND is modified into vector form of the deviatoric strain E D where the tangential strain vector E T is calculated from its third order projection tensor T as with the symmetrical fourth order identity tensor defined as I sym = 1 2 [I + I T ]. In terms of elastoplasticity, the multiplicative decomposition into elastic and plastic parts is performed using the deformation gradient F [32]. Since concrete undergoes very small elastic strains, the elastic part of the finite strain tensor could be replaced by small strains and it leads non-negative dissipation by plastic strains [29] and, thus, the additive elastic-plastic decomposition is admissible. This approach has been introduced in [28], and it yields accurate plastic behaviour of concrete at finite strains. For hyperelastic materials such as rubbers, the multiplicative formulation should be implemented instead of using the proposed additive approach.
The formulation of microplane models, in general, starts with the macroscopic free-energy function as the integration of its microplane quantities As seen in Equation (10), the microplane free-energy is the function of the Green-Lagrange strain tensor, its conjugate should be the second Piola-Kirchhoff S stress tensor. Due to material orientations, the stress quantity obtained by applying material laws to the Green-Lagrange strain tensor is the co-rotated Cauchy stress tensor as suggested in [29]. However, the second law of thermodynamics is not satisfied by the non-conjugate pair of the Green-Lagrange strain tensor and the co-rotated Cauchy stress tensor. Besides thermodynamical restrictions, the use of the second Piola-Kirchhoff stress tensor as the stress measure is appropriate for concrete, since this tensor refers to the initial configuration. Concerning the physical meaning of strength and frictional limit within the normal and shear components, the second Piola-Kirchhoff stress tensor is allowed to be the stress measure since no large elastic strain is observed in concrete, yet the finite strain may occur at the inelastic regime. For those reasons, as discussed in [26,28], the approximation of the second Piola-Kirchhoff stress tensor to be equal to the co-rotated Cauchy stress tensor is considered.
By coupling damage and plasticity, the second Piola-Kirchhoff stress tensor S is then computed as where d mic denotes the damage variable separated into compression and tension, while the superscript p indicates the plastic quantities. The constants K mic and G mic define the elastic microplane material parameters related to the bulk and shear moduli as K mic = 3K and G mic = G, respectively. Moreover, the volumetric and deviatoric projection tensors, V and Dev, are obtained as Next, the plastic strains evolve following the elastoplastic flow ruleṡ whereλ denotes the plastic multiplier, while m V and m D define the flow directions expressed as where the yield function F mic originates from the Drucker-Prager yield criterion as used in [28], which is now equipped with compression and tension caps. Meanwhile, the superscript e indicates the effective stress measure in which S e V and S e D are Furthermore, in order to visualise damage and plastic contributions, the scalar values corresponding to the damage variable d mic and the hardening variable κ mic are calculated Numerical integrations over the surface of the sphere as found in Equations (10), (11), and (17) are performed using only 21 microplanes reduced from 42 microplanes due to symmetry.

Modified Smooth Three Surface Drucker-Prager Yield Criterion
The smooth cap model within the Drucker-Prager yield criterion is implemented for capturing complex triaxial behaviour of concrete. As discussed in [33], a non-smooth yield surface leads to the singularity in the tangent terms and causes numerical instabilities with slow convergence rates. Accordingly, the smooth yield surface function depicted in Figure 1, based on the study in [34] then implemented in [24], is now extended to the finite strain framework and written in terms of the reference configuration as Here, the function F 1 contains the Drucker-Prager yield criterion with hardening with the linear hardening function defined as where D is the material constant and the hardening variable κ evolves following the flow rule aṡ κ =λ.
Moreover, the functions F c and F t describe the compression and tension caps, respectively, which are computed as where the intersection points of corresponding caps are The function H c and H t are defined by the Heaviside function for activating the respective caps as the stress state is within their domains hardened yield surface Furthermore, in accordance with [35], the calibration of several parameters in the smooth three-surface Drucker-Prager yield function is required in order to simplify the numerical implementation. By knowing concrete data in terms of uniaxial compressive strength f uc in MPa, other empirical formulas can be obtained such as the uniaxial tensile strength f ut and the biaxial compressive strength f bc Now, the parameters S 0 and α in Equation (19) can be determined as Substituting Equation (27) into Equation (28) yields α = 0.2. Moreover, f bc can be related to the yield function as (22) and (24), the material constant S C V defines the abscissa of the intersection point between the compression cap and the Drucker-Prager yield function, whereas R denotes the ratio between the major to the minor axes of the compression cap. The constant S C V is difficult to be determined. However, according to [24], S C V = − 2 3 f bc can be set as a minimum value mentioned before if no triaxial test data are provided. Meanwhile, for tension caps as in Equations (23) and (25), the material constant S T V defines the abscissa of the intersection point between the tension cap and the Drucker-Prager yield function, whereas T 0 is the initial intersection point of the tension cap with the volumetric axis. These two parameters are easier to find, and they can also be related to the empirical formula as

For compression caps as in Equations
while parameter R t controls the increase of the current intersection point T. As mentioned in [24], parameter R t can be taken as R t = 1 if no data for uniaxial cyclic tensile tests are available.

Damage Evolution Law
As mentioned earlier, the damage evolution law in this work is decomposed into compression and tension parts, d mic where γ c0 and γ t0 denote the damage thresholds for the state of compression and tension, respectively. Meanwhile, the equivalent finite strains η mic are obtained in terms of rates aṡ By activating constraints in Equations (37) and (38), this proposed model covers two conditions such as damage prevention of concrete at high confined pressure loading ifĖ p V > 0, and plastic volumetric compaction of concrete in the compression cap ifĖ p V ≤ 0.

Rate Dependency
To deliver rate dependency of concrete at dynamic cases, the proposed model is also supplemented by a viscoplastic function, a so-called consistency type formulation, adopted from [36] and applied in [37,38], which is expressed as where η vp denotes the viscosity parameter with an arbitrary value. As a consequence, considering rate dependency, Equation (18) yields where the function F 1 is modified as and, now, the intersection point T increases controlled not only by the parameter R t as in Equation (25) but also by the viscoplastic function F vp as All functions F 1 , F c , and F t as seen in Equation (40) are now a function ofκ, which are rate-dependent functions. Furthermore, the effect of rate dependency on the yield function is shown in Figure 2.

Implicit Gradient Enhancement
The microplane approach, which incorporates nonlocal damage or plasticity, or the combination of both, requires a regularisation method to eliminate pathological mesh dependencies and to yield a more stable finite element solution. As implemented successfully in previous models, the present work also employs the so-called implicit gradient-enhanced formulation for regularisation. Two governing equations are then used, such as the balance of linear momentum for the dynamic case ∇ · σ + f = ρü, (43) and the modified Helmholtz equation introducing a nonlocal field with its boundary condition ∇η m · n b = 0.
In Equation (43), ∇· is the divergence, while σ and f denote the Cauchy stress tensor and the body force vector, respectively. Mass density is defined by ρ, whereasü is the acceleration vector. Moreover, ∇ 2 and ∇ are the Laplace operator and the gradient, respectively. While the parameter c governs the nonlocal interaction range, n b describes the unit normal of outside boundaries. Furthermore, η m defines the local variable containing the equivalent strains η mic c and η mic t for compression and tension, whereasη m is the nonlocal counterpart. For adding only two extra degrees of freedom, according to [24], the variable η m is taken into consideration as Nonetheless, as found in [39,40], the over-nonlocal method can achieve a full regularisation in plastic-damage descriptions. Therefore, the equivalent strains η mic c and η mic t are enhanced aŝ It is stated that localisation still occurs at the strain softening regime if the value m = 1 is applied, which means the standard regularisation is implemented. Hence, the value of m should be taken as larger than 1 to eliminate localisation and achieve regularisation. Accordingly, η mic c and η mic t in Equations (35) and (36) are now replaced byη mic c andη mic t .

Finite Element Formulation
As standard procedure, the weight functions δu and δη m are used for obtaining the weak forms of Equations (43) and (44) With the help of the Gauss theorem and the Cauchy theorem as well as the boundary condition in Equation (45), Equations (49) and (50) become For carrying out the spatial discretisation using shape functions N, the displacement and the variational field as well as their gradients are employed Analogously, the discretisation for the nonlocal field using the nonlocal shape function N is considered as where d and E denote nodal displacements and nodal nonlocal equivalent finite strains, respectively. The discretisation using the above equations and the linearisation for Equations (51) and (52) should be used in order to be solved using an iterative Newton-Raphson method. Now, the coupled linearised system is obtained as where the residual vectors for the mechanical and the nonlocal parts are The submatrix S uu,i in Equation (58) is evaluated by the Newmark transient solver to analyse dynamic cases as where β n and ∆t denote the Newmark parameter and the incremental time step, respectively. Moreover, the mass matrix M and the submatrix K uu,i are The other submatrices in Equation (58) are then All required derivations for the tangent terms can be found in Appendix A.

Stress Return Algorithm
Two conditions may occur, either the stress state is still in the elastic region, which lies inside the yield surface, or the stress state enters the plastic regime, which is outside the yield surface. If the stress fulfils the yield criterion F mic ≤ 0, the condition is elastic. Therefore, all plastic strains at the current time step are equal to the condition at the previous time step as where the subscript n + 1 indicates the current time step, while the subscript n denotes the previous one. Then, the evaluation of the trial stresses and the yield function denoted by the subscript tr are The return algorithm needs to be enforced as F mic = 0 if the yield function is not fulfilled, which means that the step enters the plastic condition For rate dependency, the evolution ofκ is simply aṡ Using Equation (15), the flow directions m V and m D are obtained as follows Substituting Equation (77) into Equation (73) yields As can be seen from Equation (78), S e D,n+1 and S e D,tr are collinear and one can rewrite Equation (78) as Then, substituting Equation (79) into Equation (71) results in For obtaining the plastic multiplier ∆λ n+1 and the volumetric stress S e V,n+1 , Equation (80) can be solved simultaneously with the help of Equation (72), so that By deriving both Equations (81) and (82) with respect to ∆λ n+1 and S e V,n+1 , one obtains Subsequently, the updated deviatoric stress and plastic strains are

Contact Mechanism
To limit the discussion here, a general node-to-node contact mechanism is applied and it is briefly explained. This contact model is simply used since forces are computed directly at the nodes, and it provides a straightforward procedure to find the contact nodes. Moreover, a frictionless model is implemented so that only normal components of contact forces are considered, which implies that the tangential components are negligible. Since the contact model is not the main focus in the present work, no special properties are provided so that general codes for contact models can be used herein. Slight modifications are conducted in order to exclude contact nodes in damaged elements when the contact model is coupled to the adaptive element erosion procedure.

Adaptive Element Erosion
As mentioned earlier, damaged elements should be removed if they are destroyed and the impactor may penetrate into the concrete structure. To accommodate this feature, so-called adaptive element erosion is introduced in the present work supplementing the proposed material models. Here, the elimination method is implemented in a straightforward scheme as shown in Figure 3. A certain damage value as failure threshold is firstly determined prior to start the removal procedure as where d hom is a scalar damage measure which is homogenised from all microplanes as computed using Equation (17), while d val is determined as the critical value.  From the material level, the value of d hom is obtained, then it is looped over Gauss points on the element level. Here, an 8-node brick element in the three-dimensional model is used, so that eight Gauss points are taken into account. If d val is achieved in all eight Gauss points, one element is completely damaged and it can be considered to be a failed element indicated by the red-coloured square in Figure 3. Since the stiffness of the damaged element has been degraded up to almost zero, the element is now inactive, and it has to be excluded from the finite element computation without re-meshing. Hence, only active elements are involved in the computation. Furthermore, a special treatment is performed only if all neighbouring elements of a node are completely damaged. This condition causes the node connecting to its surrounding elements being unsupported and it becomes a hanging node, denoted by the yellow-coloured circle as shown in Figure 3. Now, this node is also skipped and the profile is reconstructed without re-meshing as previously.
To examine the adaptive element erosion, a four-point bending test at static loading is simulated. The proposed gradient-enhanced plastic-damage microplane model is used for modelling concrete herein. To evaluate and analyse the effect of the damage value tolerance d val to concrete responses, the values of 0.9, 0.97, and 0.99 are chosen. As can be seen in Figure 4, a good approximation is achieved by d val = 0.99 which coincides to the curve without d val , while the simulation with d val = 0.90 gives a different response among others. Using d val = 0.90, the stiffness is degraded too early and, thus, it leads to the spurious damage since the undamaged elements may categorise as damaged elements. Another example is a compact tension test in order to observe that the proposed adaptive element erosion also works well for crack branching. The simulations with three different velocities of 0.035 m/s, 1.4 m/s, and 4,3 m/s are performed according to [26] based on experiments in [41] using d val = 0.99, as obtained before. Figure 6a displays the crack paths of the compact tension specimens, while Figure 6b depicts the removed elements corresponding to its crack paths.

Coupled Contact and Adaptive Element Erosion
Since the final goal of this work is to simulate high velocity impact on concrete plates, the contact model needs to be combined with adaptive element erosion introduced previously. Therefore, the contact mechanism should be applied when the impactor hits the concrete plate surface, whereas the adaptive element erosion is used to eliminate damaged elements as the impactor may penetrate into the structure.
An initial condition with a gap between an impactor and concrete plate is shown schematically in Figure 7a, while the contact process starts when the impactor hits the plate surface as in Figure 7b and contact nodes are indicated by blue-coloured circles. The regular contact mechanism works during this process. If some elements are damaged as explained before, these elements are then eliminated and a hanging node can probably exist, see Figure 7c. Once the unsupported node is also removed, the impactor continues to move down following the applied velocity and the new contact nodes are found automatically based on the gap tolerance value as in the beginning process, see also Figure 7d. Hence, the penetration process continues and the next step proceeds repeatedly up to failure.

Model Parameters Adjustment
Several types of model parameter are investigated by numerical examples performed in here, i.e., parameters for elasticity, plasticity, damage, and nonlocal features. Young's modulus E and Poisson's ratio ν as the elasticity parameters are simply taken from the experimental data. Next, the plasticity parameters in this proposed model consist of f uc , S C V , R, R t , and D. The uniaxial compressive strength f uc is mostly given in the experimental data, whereas the parameter R can be approximated by R = X 0 / f 1 (S C V ) as mentioned in [24]. Meanwhile, the hardening D is related to the damage parameters, hence, it needs to be identified corresponding to β c and γ c0 by performing a uniaxial cyclic compression test, whereas R t is related to the damage parameter β t and identified by carrying out a uniaxial cyclic tension test. As a starting point, R t = 1 and β t 1.5β c can be considered, see [24]. Moreover, due to the fact that softening in the tensile state occurs right away subsequent to the elastic limit, the tension damage threshold is considered to be γ t0 = 0 in all simulations.
Moreover, the nonlocal parameter c is quite challenging to be identified as c = l 2 . Nevertheless, several attempts were conducted in [42][43][44] to find the parameter c. For the parameter m, as explained earlier, it should be taken larger than 1 to achieve a full regularisation. According to [45], the value of m is considered to be 1 < m ≤ 1.1 , while [40] used the value slightly beyond 1 as m = 1.005. However, m = 2.5 is used in all simulations herein. For rate dependency, an arbitrary value larger than zero of the viscosity parameter η vp gives a rate effect in the modelling. A benefit of the consistency type formulation in the viscoplastic model as in Equation (39) is to consider the constitutive law with or without strain rate effects in a straightforward manner. To produce a rate-independent model, η vp = 0 can be simply taken, whereas in order to activate the rate-sensitive formulation, η vp = 1 s/MPa is implemented in all subsequent simulations.

Four-Point Bending Test at Cyclic Loading
A four-point bending test at cyclic loading is simulated. This analysis is performed to identify the material parameters which need to be fitted with respect to experimental data in [46]. The meshed beam with a notch including boundary conditions and applied cyclic loading for the bending simulation is provided in Figure 8a,b, respectively. Table 1 provides all model parameters used in this bending test simulation. A load-displacement relation obtained by the proposed formulation using the fitted parameters is plotted in Figure 9a, while damage and plastic contributions computed by Equation (17) are depicted in Figure 9b. Moreover, mesh insensitivity of the proposed model is achieved as plotted in Figure 10a with accompanying convergence rates for the mesh with 2464 elements at chosen certain time steps shown in Figure 10b.

Parameter
Concrete Two model parameters, for instance the uniaxial compressive strength f uc and the hardening stiffness D, are observed to figure out their effects with respect to the load-displacement response. Three different values are arbitrarily chosen to perform the simulation. As can be seen in Figure 11a, the obvious influence of f uc is shown that increasing f uc affects the higher peak load. This parameter is simply taken from the experiment or, if no data provided, the value can be empirically determined corresponding to the Young's modulus E as considered in many building codes. Meanwhile, the parameter D can be identified relating to other damage parameters, β c and γ c0 , as mentioned previously. However, its effect to the hardening-softening branch of the load-displacement curve can be seen in Figure 11b. The curve slope downgrades significantly by using the lower value of D, otherwise the higher D leads to the higher hardening and the delayed softening so that the curve is going down slowly.

Homogeneous Test
In addition to identify the damage parameters D, β c , and γ c0 by simulating the uniaxial compression test as well as the parameters β t and R t by performing the uniaxial tension test, a homogeneous test can also be used to find out about the effect of strain ratesė as shown in Figure 12. Only one element is simulated subjected to uniaxial tension as illustrated in Figure 12a. The activation of strain rate effects is simply performed by taking η vp = 1 s/MPa as mentioned previously. Here, six strain rates of 0.05 s −1 , 1 s −1 , 10 s −1 , 25 s −1 , 50 s −1 , and 100 s −1 are chosen to demonstrate the influence of the proposed rate-dependent model. Due to the applied strain rates, one can see the growth of peak stresses in terms of stress-strain relations plotted in Figure 12b. In other words, by increasing the strain rate, the stress-strain curve is getting higher accordingly. Furthermore, the effect of different strain rates at cyclic uniaxial tension-compression test with respect to the stress-strain response is depicted in Figure 13.

Numerical Examples
A couple of concrete plate simulations with different thickness subjected to high velocity impacts are carried out using the newly proposed framework. Comparisons between the small and finite strain model simulated by the same material parameters are provided and compared also to the experimental data.

Thick Plate Simulation
First, an example of a thick plate is simulated to show the penetration of an impactor into the concrete plate. As reported in [47], a machine plate of 540 kg and a steel impactor of 50 kg at a height of 2 m hit a concrete plate with the dimensions of 1500 × 1500 × 300 mm. The impactor velocity is considered to be 25,000 mm/s according to [26]. Furthermore, as boundary conditions, the concrete plate is restrained at all four corners denoted by red dots as can be seen in Figure 14. This impact test is modelled using one-fourth of the whole specimen by 5241 nodes and 3272 elements due to the symmetrical condition. The whole geometry of the thick plate impact test in finite element mesh discretisation is depicted in Figure 14a, while the dimension of its impactor is shown in Figure 14b. All material properties of this impact test are provided in Table 2.
The proposed plastic-damage microplane model at finite strains including the contact mechanism coupled to the adaptive element erosion is able to represent the damage evolution including the eroded elements as shown in Figure 15a,b seen from top view of the plate surface at the chosen time step for small and finite strain simulation, respectively. Meanwhile, the visualisation of damage evaluation from the bottom side for both models is depicted in Figure 16a,b. As explained before, the damage value tolerance d val is considered to be 0.99, which means that the element is removed once the value of d hom in all eight Gauss points of one element is equal or larger than 0.99. From both Figures 15 and 16, one can see the difference between the small and finite formulation that the damage of concrete plate at small strains evolves faster than the damage in finite strain simulations. It means that the faster degradation of material stiffness occurs in numerical simulations using the small strain formulation rather than using the finite strain model. Moreover, in order to display a more realistic crack growth, Figure 17 represents only the material view of steel and concrete, which is simulated by the finite strain model. As shown in Figure 17, crack initiation starts immediately on the mid-bottom side of the concrete plate once the steel impactor hits the plate surface. Further cracks propagate diagonally towards the four fixed supports implemented in each plate corner. Crack openings get wider following the impactor penetration into the concrete plate. By these results, the use of damage values obtained from the proposed penetration model has an advantage since it is able to be coupled with the element erosion in a straightforward manner. A good visualisation of the crack growth is valuable to supplement the proposed material model, which has indeed yielded meaningful results without coupling the element erosion method altogether. Besides the damage evolution, the velocity-time and displacement-time relations are also plotted to show the differences between the small and finite strain microplane model. Analogous to the previous explanation regarding the material stiffness degradation, the velocity and displacement responses differ accordingly. The impact velocity is observed at a node located on the contact surface, while the displacement is investigated on the mid-bottom side of the plate. As one can see in Figure 18a, the impact velocity simulated by the small strain model drops immediately, so that the finite strain result has a higher residual velocity, which differs approximately 2.3 m/s. Likewise, since the material stiffness of the small strain model degrades faster as explained earlier, its displacement also deviates accordingly as shown in Figure 18b. Here, more or less 10 mm difference of the displacement is discovered between the small and finite strain simulation result. In this simulation, only a picture of damaged plate is compared to the simulation results since no other experimental data are reported. The comparison for both small and finite strain models to the experimental investigation is depicted in Figure 19. Some deviations of the crack propagation direction obtained by the experiment may occur due to the heterogeneity of material constituents in concrete, for instance diverse aggregate sizes as well as porosity between aggregates and cement paste. However, the numerical result is able to yield quite similar crack patterns compared to the experimental investigation. From Figure 20, one can see that more damage occurs surrounding the impactor in the small strain simulation, whereas the finite strain model gives a closer result to the experimental data. This observed inaccuracy may indicate that a finite strain case happens in the present impact test, so that the small strain formulation is not sufficient for modelling this phenomenon.

Small Thin Plate Simulation
Second example deals with a square small thin plate subjected to high velocity impact loading based on experiments by [48] with dimension of 610 × 610 × 30 mm. A steel impactor with a diameter of 100 mm and a height of 100 m installed in a steel frame hits the concrete plate by an initial velocity of 12,300 mm/s. As boundary conditions, fixed supports of 30 mm wide denoted by red dots are applied on each side of the plate specimen. Figure 21a displays the whole geometry of the impact test in finite element mesh discretisation, whereas Figure 21b provides the dimension of the steel impactor used in the present work. Due to the symmetry, only a quarter of the geometry is simulated with 5463 nodes and 3440 elements. Moreover, all material properties of this impact test are shown in Table 3.
1 -Similar to the previous example, the damage value tolerance d val is considered to be d hom = 0.99 obtained from the material level. As discussed previously, the closer d val is approaching 1.0, spurious damage in numerical simulations can be avoided. Furthermore, the damage evolution obtained by the proposed finite strain model is displayed in Figures 22 and 23 for the top and bottom side of the concrete plate, respectively. Damage initiation starts at the middle of the plate after the impactor hits the plate surface. Hereafter, damage evolves continuously, in particular surrounding the steel impactor. Then, it gets wider following the impactor velocity. This condition causes the destruction of elements at the mid-surface as shown in Figure 22, while the impactor starts to penetrate into the plate. Subsequently, as one can see in Figure 23, cracks propagate diagonally towards the four corners of the plate emerged at the bottom side. Moreover, vertical and horizontal cracks are then observed. In a similar manner to the preceding example, the crack growth visualisation for this impact simulation is then provided in Figure 24a,b shown from the top and bottom side of the concrete plate, respectively. For the sake of a clearer and more realistic view, Figure 24 displays only the material colour of the steel and concrete. From the figure, additional crack propagations as seen in Figure 23 are not represented in Figure 24b since the damage value d hom of the diagonal cracks as well as the vertical and horizontal cracks have not achieved 0.99 yet.
To investigate the difference between small and finite strain simulations, the impact test is performed using three different velocities, i.e., 12.3 m/s, 16.5 m/s, and 20.3 m/s. First investigation is the velocity-time relation for both models compared to the experimental data in [48], see Figure 25. At the lower velocity, v = 12.3 m/s, coinciding curves of small and finite strain results are observed. Meanwhile, both curves differs slightly at higher velocity than before, with v = 16.5 m/s. Furthermore, at v = 20.3 m/s, the deviation of both models gets larger with the difference in velocity of approximately 1 m/s. By performing all simulations with the three different velocities, the finite strain microplane model gives better agreement compared to the experimental investigations than the small strain formulation. Here, one can deduce that the higher the impact velocity, the larger the deviation of residual velocity between the small and finite strain simulation. Again, finite strain condition may occur in the impact test simulation at higher velocities, so that the deviation is observed using the small strain simulation.  Figure 27b is almost identical for both models. Moreover, at the higher impact velocity of 20.3 m/s as shown in Figures 26c and 27c, damage patterns obtained by the small strain simulation are focused only in the middle part of the plate. No other crack propagation at the top side is observed unlike the finite strain result. Furthermore, by comparing the simulation results numerically to the existing experimental investigations, one can observe that the finite strain microplane model yields a better accuracy in predicting the damage pattern of the present high velocity impact test than the small strain formulation. Apparently, as can be seen in Figure 28, the difference of crack prediction for both models compared to the experiments is caused by the heterogeneous nature of the concrete material due to its composition as well as aggregate size. In addition, imperfection of the experimental set-up can possibly happen during the impact test in the reality, whereas the numerical simulations are always carried out by considering the model testing in a perfect condition. Throughout the simulations and result analyses of the small thin concrete plate at the three different velocities, both small and finite strain microplane models are used to perform the simulation. As one can see from the aforementioned investigations, the small strain approach provides good results and gives almost identical responses at the impact velocity of 12.3 m/s compared to the proposed finites strain formulation. Hence, both models are able to simulate this condition. However, each formulation yields different results during the impact test at higher velocities, such as 16.5 m/s and 20.3 m/s. The observation in terms of velocity and crack pattern simulated by the finite strain model gives a better approximation compared to the experimental data than the small strain formulation. In this regard, a finite strain case probably occurs so that the small strain microplane model cannot suffice to demonstrate the accurate response.

Conclusions and Outlook
Formulating a robust material model for concrete using the well-known microplane approach extended to the finite strain framework as a main objective of the contribution at hand is achieved. The proposed model can be used to analyse concrete structures numerically at various loading states particularly for dynamic cases within implicit finite element codes. Constitutive models are developed in the present work involving a combination of both damage and plasticity. The newly proposed model combines the plastic-damage microplane approach at finite strains due to the fact that large deformations occur in concrete at high velocity impacts. The conjugate pair of the Green-Lagrange strain tensor and the second Piola-Kirchhoff stress tensor is used to extend the small strain microplane model to the finite strain regime. A smooth three-surface Drucker-Prager yield function is introduced and supplemented by rate dependency to deliver strain rate effects at dynamic loading. The proposed model is then used to simulate impact tests with high velocities on concrete plates after being implemented together with a contact mechanism and an adaptive element erosion during penetration.
A penetration model is required in such impact tests due to the fact that an impactor may penetrate into concrete structures at high velocity impact. Damage values obtained by the proposed material model are used as failure criterion to eliminate destroyed elements of concrete during the impactor penetration. To simulate this phenomenon, the adaptive element erosion is also equipped with a general contact mechanism. The accuracy and capability of the new framework are then evaluated by comparing with the experimental data. As a result, the proposed microplane model at finite strains gives better agreement to the experimental investigations than the small strain version since finite strain situation may occur in the present impact test simulations.
The research work at hand is able to represent concrete responses in a good approximation and reliable results compared to existing experimental investigations. However, some improvements can be pursued to obtain a more sophisticated tool for modelling concrete structures in larger examples and different types of loading. For instance, the adjustment of model parameters in an advanced manner is useful to facilitate numerical simulations of other cases. Herein, the plastic parameters can be simplified by identifying the uniaxial compressive strength of concrete. Nonetheless, determination of damage and nonlocal parameters in each numerical example are not simple. It probably needs an optimisation algorithm in order to find accurate values, yet indeed, computational costs may increase accordingly.

Acknowledgments:
The authors gratefully acknowledge the financial support by the Institute for Structural Analysis, TU Dresden for finishing this research topic.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Algorithmic Tangent
Algorithmic tangent moduli for the proposed model contain both elastoplastic and damage tangent terms. Once the elastoplastic tangent is obtained, the total tangent is then taken into account by considering damage.

Appendix A.1. Effective Elastoplastic Tangent
The differentiation of the yield function is required to compute the effective tangent terms δF mic = n V δS e V + n D · δS e D + n κ δκ = 0, with Next, the volumetric and deviatoric stresses are differentiated as Knowing Equation (21), Equation (A3) becomes Then, one can rearrange Equations (A4) and (A5) as Furthermore, Equations (A6) and (A7) can be written as where h 1 , h 2 , and H 3 are, respectively, (A12) Now, Equations (A8) and (A9) are substituted into Equation (A1) as follows or one can write in a convenient way where h 4 is Then, Equation (A14) can be substituted into Equations (A8) and (A9) Subsequently, multiplying Equations (A16) and (A17) by V and Dev T , respectively, one obtains The differentiation of the effective second Piola-Kirchhoff stress tensor for one microplane is obtained by summing Equations (A18) and (A19) as Now, δS mic e is used to obtain the total tangent incorporated with damage.

Appendix A.2. Total Tangent
From Equation (11), one can write where S mic e originates from S mic which is differentiated with respect to both local and nonlocal quantities in terms of compression and tensions part, so that