In order to simulate large elastoplastic deformation of metal, the rate of deformation is decomposed into elastic and plastic parts as D = D e + D P (1) where D e is the elastic deformation rate, and D P is the plastic deformation rate. The elastic constitutive relationship is defined as D e = 1 + ν E σ ∇ J − ν E tr ( σ ∇ J ) I (2) where E   is Young's modulus, ν is Poisson's ratio, σ ∇ J is the Jaumann rate of Cauchy stress, and I is the second order identity tensor.For isotropic materials, the yield function is defined as f ( σ , ε p , T ) = σ e q − σ y ( ε p , T ) (3) where σ e q   is the equivalent stress, ε p is the equivalent plastic strain, T is the temperature, and σ e q = 3 2 σ ′ : σ ′ (4) σ ′ = σ − 1 3 t r ( σ ) I (5) where σ '   is the stress deviation, and the plastic strain increment is related to the stress deviator.The yield stress function is defined as (Aune et al., 2017). σ y = [ σ 0 + ∑ i = 1 2 Q i ( 1 − exp ( − C i ε p ) ) ] [ 1 + ε ˙ p ε ˙ 0 p ] c [ 1 − T ∗ m ] (6) where σ 0   is the initial yield stress, Q 1 , Q 2 , C 1 , C 2 are the material hardening parameters, ε ˙ 0 p is the reference plastic strain rate, and the expression of T ∗ is given as T ∗ = T − T r T m − T r (7) where T r is the ambient temperature, and T m is the melting temperature of the material. Assuming that the plastic potential function coincides with the yield function, the plastic deformation rate tensor can be expressed as D P = ε ˙ p ∂ f ∂ σ = 3 2 ε ˙ p σ ′ σ e q (8) Since materials undergo rapid plastic deformation subjected to impulsive loads, the adiabatic process is assumed. The rate of temperature rise is defined as T ˙ = χ ρ C p σ : D p = χ ρ C p σ e q ε ˙ p (9) where χ   is the Taylor-Quinney coefficient, ρ is the density, and C p is the specific heat. The material parameters used in the metal constitutive relation are shown in Table 2.

Because of geometrical feature of the FML plates having a thickness significantly smaller than the other dimensions, the aluminum and composites layers are modeled using shell elements. Composite layups can be composed of plies made of different materials in different orientations. To simplify the model, the shell composite layup model in ABAQUS that discretizes only the reference surfaces for each ply is used, and thus fibers and matrix are not explicitly modeled.

The failure criteria of fiber-reinforced composite materials adopt the Hashin’s theory (Hashin and Rotem, 1973; Hashin, 1980) that has been implemented in the commercial finite element package ABAQUS. Four different modes of failure are considered including fiber rupture in tension, fiber buckling and kinking in compression, matrix cracking under transverse tension and shearing, and matrix crushing under transverse compression and shearing. The initiation criterion is expressed as d f = { ( σ ^ 11 X T ) 2 + α ( τ ^ 12 S L ) 2 , σ ^ ≥ 0 ( σ ^ 11 X C ) 2 , σ ^ 11 < 0 ( σ ^ 22 Y T ) 2 + ( τ ^ 12 S L ) 2 , σ ^ 22 ≥ 0 ( σ ^ 22 2 S ) 2 + [ ( Y C 2 S T ) 2 − 1 ] σ ^ 22 Y C + ( τ ^ 12 S L ) 2 , σ ^ 22 < 0 where X T   denotes the longitudinal tensile strength, X C   denotes the longitudinal compressive strength, Y T   denotes the transverse tensile strength, Y C   denotes the transverse compressive strength, S L   denotes the longitudinal shear strength, S T   denotes the transverse shear strength, α   is a coefficient that determines the contribution of the shear stress to the fiber tensile initiation criterion, and σ ^ 11 , σ ^ 22 , τ ^ 12 are components of the effective stress tensor. A value of d f ≥ 1.0 indicates that the initiation criterion has been met.

When damage begins to accumulate, each failure mode has its own damage evolution and a damage index in the range of 0–1 as the measure of damage. Using the energy-based linear damage evolution law, the fracture energies   G f t c ,   G f c c , G m t c   and G m c c associated with the energy consumed in the damage process of each mode, are set in the simulations. In addition, material models with stiffness degradation usually lead to difficulties in the convergence of analyses, so the viscosity regularization scheme is used to facilitate convergence. Viscosity coefficients,   η f t , η f c , η m t and η m c   for each damage mode, are summarized in Table 1.

The degradation of bonding adhesive layer between metal and composite layers in an FML significantly reduce the rigidity and eventually lead to interlaminar delamination. Thus, the adhesive behavior of bonded neighboring layers must be taken into consideration. In this study, the surface-based cohesive behavior model available in ABAQUS is used.

The linear elastic traction-separation behavior is assumed in the surface-based cohesive model. Material properties used to define the elastic behavior of adhesive are summarized in Table 3. The debonding process consists of three phases, i.e., damage initiation, damage evolution and complete failure. This study uses the quadratic separation criterion to characterize damage initiation represented as { 〈 δ n 〉 δ n 0 } 2 + { δ s δ s 0 } 2 + { δ t δ t 0 } 2 = 1 , 〈 δ n 〉 = { δ n       δ n > 0 0           δ n < 0 (11) where δ n 0 ,   δ s 0 , and δ t 0   represent the values of separation, when the separation is either purely along the contact normal or purely in the first or the second shear direction, respectively. δ n , δ s , and δ t   represent distances of separation in different directions.

The law of damage evolution describes the degradation rate of stiffness once the initiation criterion is reached. The energy-based exponential damage evolution law is adopted, and the power law criterion is used to describe the interaction of mixed modes as (Camanho and Davila, 2002) { G n G n C } 2 + { G s G s C } 2 + { G t G t C } 2 = 1 (12) where G n C , G s C and G t C refer to the critical fracture energies required to cause failure in the normal, the first, and the second shear directions, respectively. G n , G s and G t represent the work done by the tractions and their conjugate separations. The mixed-mode fracture energy is G C = 1 ( { m 1 G n C } α + { m 2 G s C } α + { m 3 G t C } α ) 1 α (13) where m 1 = G n / G T , m 2 = G s / G T , m 3 = G t / G T , and G T = G n + G s + G t . Exponential softening is expressed by the evolution of damage variable D as D = ∫ δ m 0 δ m f T eff G C − G 0 d δ (14) where δ m f denotes the displacement at complete failure, δ m 0   denotes the displacement at the beginning of the damage, T eff   and δ   are the traction and separation, respectively, and G 0   is the elastic energy at damage initiation.

A primary goal of this study is to investigate the anomalous dynamic response that is highly sensitive to the applied loads, hence the load function is of critical importance. While many studies use simplified pulse shapes such as triangles or rectangles in their numerical simulations, an in-house code specifically developed for modeling the impulsive pressure is ported to commercial finite element software via user defined subroutines in the present study.

The pressure-time history is comprised of positive and negative phases as shown in Figure 2. The pressure profile is characterized by the ambient pressure p 0 , peak reflected overpressure p rα with the subscript α denoting the angle of incidence, negative overpressure p r , min , arrival time t a , positive duration t d + , negative duration t d − , and decay parameter b .

The required input parameters for the impulsive load function are the location of source and equivalent charge weight. The scaled distance using the Hopkinson-Cranz scaling law can be used to describe the intensity of shock loading as Z = R W 1 3 (15) where R is the stand-off distance and W is the equivalent weight. The positive phase of the pressure-time history can be described by the modified Friedlander equation expressed as p ( t ) = p 0 + p rα ( 1 − t − t a t d + ) exp ( − b ( t − t a ) t d + )                 t a < t < t a + t d + (16) The reflected pressure considering the angle of incidence α can be calculated by p rα = p inc ( 1 + cos ⁡ α − 2 ⁡ cos 2 ⁡ α ) + p ref ⁡ cos 2 ⁡ α where p inc and p ref are incidence and reflected pressures, respectively. By integrating the positive phase of pressure over time, the impulse i r + is obtained i rα + = ∫ t a t a + t d + p rα ( 1 − t − t a t d + ) exp ( − b ( t − t a ) t d + ) = p rα t d + b 2 ( b − 1 + exp ( − b ) ) (17) Since p inc , p ref , i rα + are given as functions of the stand-off distance R and weight W , the decay parameter b can be implicitly solved. The negative phase of the pressure is critical to the occurrence of counter-intuitive behavior; thus, it must be considered in the loading function. A bilinear approximation of the negative phase is adopted (Aune et al., 2017), and the full expression of pressure is defined as p ( t ) = { p 0 , t < t a p 0 + p rα ( 1 − t − t a t d + ) exp ( − b ( t − t a ) t d + ) , t a < t < t a + t d + p 0 − p r , min ( t − ( t a + t d + ) c t d,lin − ) , t a + t d + < t < t a + t d + + a t d,lin − p 0 − p r , min ( t − ( t a + t d + + c t d − ) ( 1 − c ) t d − ) , t a + t d + + a t d,lin − < t < t a + t d + + t d,lin − p 0 , t > t < t a + t d + + t d,lin − (18) Where t d − is the time duration of negative phase, which is typically longer than that of the positive phase, and the coefficient c is introduced to control the proportion of the negative pressure rise time in the entire negative pressure duration.

In this subsection, the dynamic response of an aluminum plate under impulsive loading is studied and compared against experimental results available in literature. The dimensions of the plate are 300 × 300 × 0.8   mm 3 with edges fully clamped. Two impulsive loads of varying intensity defined in terms of the Hopkinson-Cranz scaled distance, Z = 1.46   m / kg 1 / 3 and Z = 1.82   m / kg 1 / 3 , are considered.

The deflection at the plate center is shown in Figure 3, indicating radically different dynamic responses of the plate subjected to impulsive loads of varying intensity. The final deflection of the plate under the impulsive load of Z = 1.46   m / kg 1 / 3 is around 25 mm. On the other hand, the impulsive load of Z = 1.82   m / kg 1 / 3 yields a negative final deflection around - 20 mm with the plate resting on the opposite direction to the incoming impulsive load. This phenomenon is referred to as anomalous or counter-intuitive behavior. Figure 4 shows the final configurations of the plate. Different signs in the out-of-plane (i.e., the z -direction) displacement can easily be noticed. The predicted permanent deflections are validated against the experimental results reported in (Aune et al., 2016). Good agreement between numerical predictions and experiments is reached. It must be emphasized that since the counter-intuitive behavior only occurs in a tight range of impulsive loads, the accuracy of loading function is crucial to the prediction of counter-intuitive behavior. A slight discrepancy might lead to wrong results in predictions of anomalous responses. Consistency with the experimental results in both intuitive and counter-intuitive behaviors underpins the following studies of FMLs subjected to impulsive loading.

We then study dynamic responses of FML plates and examine the difference in comparison to the responses of the monolithic metal plate. The geometries and material properties have been summarized in Materials and Methods. The thickness of the FML plates under consideration is identical to that of the aluminum plate studied in Dynamic Response of a Monolithic Metal Plate, so a quantitative comparison in dynamic responses can be made between monolithic and FML plates.

Figure 5A plots the deflection time histories at the plate center with the impulsive load of Z = 1.82   m / kg 1 / 3 . Once reaching the maximum positive deflection, all three FML plates rebound toward the opposite direction to the incoming impulsive load. After about two cycles of oscillation from the lowest negative deflection, FML plates progressively arrive at their final configuration. Like the monolithic aluminum plate, the FML plate with a 0.18  mm composite layer reaches a final configuration with a negative deflection at the plate center. But the deflection is reduced by 42 % from 19.16  mm to 11.20  mm. Two FML plates with a 0.22  mm composite layer approach their initial configuration of zero deflection at the plate center. What stands out in Figure 5B is that the thickness of composite layer in an FML only mildly affects the maximum positive deflection but imposes a strong effect on the final configuration.

Figure 6 shows the final deformation profiles of FML plates. For the 0.31 / 0.18 [ 0 ° ] / 0.31 FML plate, the deformation profiles along the x - and y -directions demonstrate a significant difference due to the directionality of composite layer. On the other hand, the FML with two orthogonally laid plies of fibers, i.e., the 0.22 / 0.36 [ 0 ° / 90 ° ] / 0.22 plate, demonstrate a much smaller difference in the deformation profiles along the x - and y -directions. More importantly, the slope of the y -direction deformation profile near the edge is considerably higher in the FML plate with only unidirectional lamina. As large deflection slope might eventually lead to tearing of the plate at the boundary, special attention in structural design is thus desired.

Given the considerable difference in the slope of deflection near the plate edges, we further inspect the state of stress triaxiality defined as the ratio of the hydrostatic mean stress to the von Mizes equivalent stress, i.e., σ ∗ = σ H / σ e q in which the hydrostatic mean stress defined in terms the principal stresses as σ H = ( σ 1 + σ 2 + σ 3 ) / 3 . It is generally accepted that shear fracture is the failure mechanism in the negative stress triaxiality regime and fracture due to void formation in the high stress triaxiality regime (Bao and Wierzbicki, 2004). Figure 7A plots the evolution of stress triaxiality at the mid-point of edge perpendicular to the fiber direction. Upon the arrival of pressure pulse, the stress triaxiality abruptly increases to a state of 0.6 in the vicinity of the top cohesive interface, and a negative stress triaxiality of −0.4 is reached near the bottom interface. During the dynamic responses to the impulsive loading, the stress triaxiality intensely oscillates between positive and negative regimes, indicating different failure mechanisms can be activated once a threshold value, such as critical strain, is reached. When the plate eventually rests on its final configuration, the stress triaxialities are positive on the top and negative on the bottom interfaces, respectively. Figure 7B shows the evolution of stress triaxiality at the mid-point of edge parallel to the fiber direction. Positive stress triaxialities develop near both top and bottom cohesive interfaces upon the arrival of pressure pulse. Two phases of oscillation in triaxiality are noticed, first in the range of 0.6 to −0.6 and −0.2 to −0.6 subsequently. When the plate rests on its final configuration, the stress triaxialities are negative on the top and positive on the bottom interfaces, respectively. The clear differences in the state of stress triaxiality observed in Figures 7A,B is closely related with the fiber orientation.

Next, we examine the dynamic responses of FML plates subjected to higher intensity of impulsive loading with the scaled distance Z = 1.46   m / kg 1 / 3 . Figure 8 shows the deflection versus time at the plate center. Recall that the permanent displacement of the aluminum plate is in the intuitive direction under the impulsive load of Z = 1.46   m / kg 1 / 3 (see Figure 3). The 0.31/0.18 [0°]/0.31 plate, however, demonstrate a counter-intuitive behavior with the permanent deflection in the opposite direction to the incident impulsive loading. Hence, it is inferred that the dynamic responses of monolithic and FML plates can be entirely different. Furthermore, unlike the 0.31/0.18 [0°]/0.31 plate, both 0.22/0.36 [0°/0°]/0.22 and 0.22/0.36 [0°/90°]/0.22 plates demonstrate an intuitive permanent deformation with positive deflection at the plate center, suggesting significant influences of FML configurations on dynamic responses. The presence or absence of counter-intuitive behavior is affected by many factors, such as material properties (Li et al., 2008), load magnitudes and profiles (Symonds and Yu, 1985; Aune et al., 2017), boundary conditions (Bassi et al., 2003), etc. The ultimate state is extremely sensitive to the initial conditions (Ng and Daolin, 2002). The analyses in the present study highlight additional, but potentially enormous, intricacy induced by FMLs to the nonlinear responses that deserve attention in structural design and applications.

For all three FML plates subjected to the scaled pulse Z = 1.82   m / kg 1 / 3 , composite layers suffer matrix tension damage as shown in Figure 9. Comparing Figures 9A,B reveals the influence of the ratio of the metal to composite layer thickness while the total plate thickness is kept the same for both cases. It is found that the 0.31 / 0.18 [ 0 ° ] / 0.31 plate suffers less damage than the 0.22 / 0.36 [ 0 ° / 0 ° ] / 0.22 plate. Therefore, a thicker aluminum layer attenuates the matrix damage resulted from impulsive loading. Additionally, reducing aluminum layer thickness tends to shift damage from the plate center to the edges along the direction perpendicular to the fiber direction as shown in Figure 9B. Two orthogonally laid plies have a beneficial effect in reducing damage as shown in Figure 9C. In general, the plate center and clamping edges are of high propensity to accumulate damage. To further reduce the extent of damage, the thickness of the FML plates needs to increase.

The degradation and eventual failure of bonding between two cohesive surfaces are major concerns for applications using FMLs. Figure 10 shows the damage evolution at various locations of the 0.22/0.36 [0°/0°]/0.22 plate. Overall, the top cohesive interface, which is closer to the source of impulsive loading, suffers severer damage than the bottom cohesive interface; plate edges sustain more damage than the center, which is closer to the source though. Although damage is produced, no debonding takes place for both Z = 1.82   m / kg 1 / 3 and Z = 1.46   m / kg 1 / 3 .

It is worth pointing out that although the impulsive loading is applied normal to the plate surface, localized sliding can occur between layers. To illustrate this phenomenon, the displacements along the x - and y -centerlines of the plate are plotted in Figure 11. A particular noteworthy finding is that localized sliding can arise at locations that are often overlooked (see the rectangle annotations in Figure 11 The displacements along the x- and y-centerlines of the 0.22/0.36 [0°/0°]/0.22 plate under the impulsive load of Z = 1.46  m / kg 1 / 3 . ) especially if the final permanent deflection is small or even close to the zero. The localized sliding due to inconsistency in in-plane displacements might lead to fiber kinking and several other failure types (Greenhalgh, 2009).

Lastly, in order to clearly understand the damage of FML during the dynamic response, the energy dissipated by damage is plotted as shown in Figure 12. In the case of impulsive load Z = 1.82   m / kg 1 / 3 , the energy dissipation due to damage is very close in value for the 0.31/0.18 [0°]/0.31 and 0.22/0.36 [0°/0°]/0.22 FML plates although the configuration of layer thickness is different; the FML plate with two plies of unidirectional fibers laid orthogonally, on the other hand, induces less damage energy. Under the condition of a higher impulsive load Z = 1.82   m / kg 1 / 3 , the energy dissipated by damage increases to twice that of the previous load. Additionally, the difference in energy between 0.31/0.18 [0°]/0.31 and 0.22/0.36 [0°/0°]/0.22 FML plates can be noticed, suggesting a better protective performance of thicker composite layer.