Non-Linear Thermomechanical Modeling of Materials

Plasticity and Fracture in Materials Under High Impact Loading

Under high-speed impact loading, materials experience significant stress and strain which can result in permanent deformation. This work investigates two dominant material failure mechanisms:

Plasticity, used to describe deformation in ductile materials
Fracture, which characterizes crack formation and propagation

Using MOOSE’s finite element framework and thermomechanical models, these failure modes are analyzed in response to dynamic loading to better understand how materials absorb and dissipate energy, and how failure localizes under impact.

While commonly idealized, material microstructures introduce heterogeneity that plays a critical role in the localization of energy which leads to either plasticity or fracture nucleation depending on the failure mechanism of focus. At the meso- and macro-scales, microstructure features such as voids, grain boundaries, and multi-phase inclusions introduce structural heterogeneity [2]. To capture the effects of a microstructure on the material response, I map spatial variation to material parameters (using scaling factor k), with a particular focus on elastic stiffness (K0) and critical surface fracture energy (Gc).

High resolution microstructural models provide greater accuracy but come at a high experimental and computational cost. To improve efficient, I reduce the resolution of these features and evaluate how this impacts the simulation accuracy and overall predictive capabilities of the material response.

MATERIAL Models

The ability to model material behavior is critical when developing simulations as it greatly impacts the quality of the results. This section highlights the solid and thermomechanical models that capture the nature of these materials responses. These material models are a combination of systems seen in material science literature, and they have been calibrated with both experimental and atomistic simulations (molecular dynamics).

The constitutive model decomposes strain into elastic and plastic parts, includes anisotropic elastic stiffness, uses an artificial viscosity to increase numerical stability, and relates the pressure to the temperature and change in volume using an equation of state.

Johnson Cook plasticity is a temperature and strain rate dependent model that includes the effects of strain hardening. As the strain rate increases, the yield stress increases, while increases in temperature results in a lower yield stress.

The phase field damage model (fracture model) compares the strain energy density to the material’s critical surface fracture energy to determine the rate of damage to the material. A brittle fracture approach is used with only strain energy from tensile and shear stress contributing to material damage.

The time rate of change of temperature is described by the thermal transport model, where plastic dissipation, artificial viscosity, and thermal elastic heating sources as well as heat conduction contribute to the evolution in temperature.

We used an infinitesimal strain model that assumes small deformation and can be expressed in terms of position, x, and displacement, u [1].

The strain, ε, can be additively decomposed into elastic, εp, and plastic, εe, parts [2].

The stress, σ, can be decomposed into volumetric, σvol, and deviatoric-volumetric, σdev, parts.

Where ℂ is the fourth order elasticity tensor and is the bulk modulus, calculated as:

At high impact speeds, two distinct regions form within the material, a disturbed region that has a state of high kinetic energy and pressure, and an undisturbed region. The shock front is the boundary between these regions of different states. An equation of state (EOS) describes the pressure associated with the state of the material as a result of the shock wave.

The Hugoniot curve relates the speed of the shock front, Us, to the particle speed, Up, within the material [4]. In this work, this relationship is represented by:

where c0 is the speed of sound in the material.

Computational simulations of shocks introduce numerical instabilities at high shock speeds. An artificial viscosity stress, σav, is introduced to address shock discontinuity and the instabilities that arise from extreme pressure differences at the shock front [5]. The artificial viscosity stress takes the form:

Where hmax is the element size, C0 is the Neuman coefficient used to stabilize the wave front, C1 is the Landshoff coefficient used to reduce oscillations, and J is the Jacobian [4]. The artificial viscosity is only non-zero under the condition that the material is under volumetric compression (J > 1).

The volumetric component of stress is then computed as:

The fourth order elasticity tensor represents the anisotropic stiffness of the material. The elasticity tensor is assumed to exhibit symmetry, resulting in ℂijkl = ℂjikl = ℂijlk = ℂklij. This enables the elasticity tensor to be expressed as a second order stiffness tensor in Voignt notation.

The balance of linear momentum is dependent on the stress, density, and acceleration, ü [2].

The Johnson Cook (JC) plasticity model assumes isotropic plasticity where the yielding stress, σy, is dependent on the equivalent plastic strain, εp, strain rate, ε̇p, and temperature, T, of the material [6], [7].

The material constants in the JC plasticity model are the yield strength in standard conditions, A, strain hardening coefficient, B, strain rate sensitivity constant, C, hardening exponent, n, and thermal softening exponent, m [7].

The dimensionless strain rate, ε̇*, is based on the current strain rate, ε̇, and the reference strain rate, ε̇ref. The homologous temperature, T*, is dependent on the current temperature, T, melting temperature, Tm, and reference temperature, Tref.

At each time step, an initial stress, σtrial, is calculated using constitutive model, where the strain increment is assumed to be purely elastic. A radial return algorithm is used in conjunction with a Newton Raphson method to identify the plastic strain that occurred and correct the stress back to the yield surface.

A deviatoric stress tensor, σtrial_dev, is extracted from the initial stress tensor, and its second invariant, J2, is calculated.

The second invariant of deviatoric stress is compared to the yield stress to determine the proximity of the stress state to the yield surface.

If f(J2) ≤ 0, then the new deformation is purely elastic. However, if f(J2) > 0, then the stress and elastic strain is overestimated, and plastic correction is required to return to the yield surface and account for the contribution of plastic strain.

The heat equation describes the change in temperature, T, over time due to the contribution of thermomechanical processes.

The temperature is dependent on the volumetric isobaric specific heat, Cv, density, ρ0, thermal conductivity, k0, and three sources of heat generation, plastic dissipation heating, qp, artificial viscosity heat generation, qav, and thermal elastic heating, qeos. The last term in the heat equation describes the temperature diffusion caused by temperature gradients.

The presence of plasticity introduces heat generation due to plastic dissipation:

Where ε̇p is the plastic strain rate. Tuning parameters βp and βav represent the fraction of energy from plasticity and artificial viscosity which contribute to heating and were found through molecular dynamics (MD) simulations.

Heat generation due to artificial viscosity is described by:

The contribution of this heat source only occurs when the system is in volumetric compression (J < 1).

The thermal elastic heating describes the heat generation due to thermal stresses in the material.

Where the thermal elastic heating has a negative relationship with the volumetric strain rate, resulting in increasing heat generating at high compressive volumetric strain rates.

A phase field damage model is used to model fracture within the material. A continuum field variable, c(x,t), tracks the evolution of damage, with damaged material represented by c(x,t) = 1, and undamaged material as c(x,t) = 0. Below is a phase field fracture model based on the works of Dandekar [8] and Zhang [9].

A fourth order positive projection tensor of eigen-decomposition, ℙ+, is used to find the positive projection stress. The derivation of ℙ+ is based on [10], [11]. The positive and negative stress projections for the deviatoric-volumetric stress tensor is additively decomposed.

The positive part of strain energy density, a+, is:

The negative part of strain energy density, a-, is:

The total strain energy density a0(σ,ε) is:

Where the degradation is dependent on the damage, c, the positive and negative strain energy densities, and kr, a constant used to prevent the total degradation from being equal to zero.

The elastic free energy density, We, can be found for the volume, Ω:

The fracture surface energy, Wf, is defined on the surface, Γ:

where Gc is the critical surface fracture energy. The crack surface, Γ, is based on a diffuse delta-function, γ(c,∇c), describing the transition zone between the damage and undamaged material.

Thus, the fracture surface energy takes the following form:

The kinetic energy of the material is K(u˙), where is the velocity vector. This results in the balance of momentum equation as defined in Equation 1.12.

To prevent healing of the damaged zone and ensure thermodynamic consistency:

This results in the following conditional requirements:

These conditions can be addressed using an extended Lagrangian as define below.

The time rate of change of the phase field damage variable, c, is described below [12].

Where η is a computational coefficient controlling relaxation of the system used in the MOOSE finite element software to control the convergence [8]. The Macaulay bracket function below is used to prevent negative time rate of change in damage.

Boundary Conditions

A controlled impact velocity is applied to one edge of the material to compress the sample, while the left and right boundaries are fixed in the x-direction, and the top boundary is fixed in the y-direction. Our research is focused on high-speed impacts, promoting high strain rates and large strain energy density. To solve the system using a two-dimensional mesh, we assume plain strain.

Impact speeds: 0.1 - 1.0 km/s

Thermal Plasticity Response

To explore the effects of variable microstructures on the thermo-mechanical response of a material, the elasticity tensor was varied as a function of space within a material sample. The resulting deformation is described by anisotropic elasticity, an equation of state, thermal transport, and Johnson Cook plasticity.

Key Take Aways:

Higher temperatures form in regions of lower elastic stiffness due to higher rates of plastic and thermal elastic heating.
Resolution reduction leads to underestimation of high temperatures, and overall decrease in standard deviation of temperature.

The temperature fields were found for each simulation by capturing the state of the material when the shock has approximately reached the end of the domain. The resulting temperature fields show resemblance to the stiffness fields for each resolution, with softer regions demonstrating higher temperatures. Regions of lower stiffness experience greater elastic and plastic strain, increasing thermal elastic and plastic heating. Reducing the resolution of the microstructure causes a decrease in the standard deviation of stiffness and thus a smaller temperature range is exhibited. The elastic stiffness is expressed as a bulk modulus, K0.

The effect of impact velocities on the temperature field was also explored. The figure below shows a comparison of the temperature fields at all three resolutions and impact velocities of 0.1, 0.4, and 1.0 km/s. Lower impact velocities result in lower strain rates and thus a decrease in the overall temperature field. At all impact velocities the effects of homogenization due to resolution reduction are evident, where lower resolutions underestimate the temperature of a region.

Fracture Response

Both the critical surface fracture energy and elastic stiffnesses were independently examined to determine their influence on fracture initiation under impact and thermal loading. The deformation is simulated in MOOSE and is described by the fracture model which includes anisotropic elasticity, an equation of state, and the phase field damage model.

Key Findings:

Fracture nucleation occurs in regions of high critical surface fracture energy and stiff regions with significant surrounding stiffness gradients.
With high strain energy density, underlying microstructures have less influence over fracture propagation paths, and branching can occur.

In shock-based loading, fracture initiation occurs primarily along the left boundary edge where the impact velocity is applied. At the early stages of impact, there is high localized pressure near the boundary, resulting in significant strain energy that leads to fracture nucleation.

To better understand the effect that a heterogeneous microstructure has on material fracture, I independently analyze spatial variance of critical surface fracture energy and elastic stiffness.

Spatial Variation of Critical Surface Fracture Energy

For a microstructure with spatial variation of Gc, regions of lower will accumulate damage at a higher rate, acting as nucleation sites for fracture. The test below highlights the role that Gc plays in the development of the fracture network. Fracture nucleation occurs at the interfaces of regions with dissimilar Gc, with propagation continuing at an angle. Regions of high Gc, demonstrate a resistance to fracture, resulting in fractures propagating away from the region. Branching is observed in both cases, with nucleation sites developing into two fractures which propagate at an angle. At high fracture propagation speeds, there is sufficient energy to generate greater fracture surface area, thus fracture bifurcation can occur [13].

Shock simulations were conducted at impact velocities of 1.0 km/s and 0.4 km/s to assess the influence of impact speed on the fracture behavior of a microstructure. The time evolution of the fracture response highlights a dependence of fracture nucleation on the availability of strain energy, with lower impact speeds resulting in lower strain rates and thus fewer fractures.

Due to high strain energy density, the Gc field has limited influence on the direction of fracture propagation, with fractures moving through regions regardless of the Gc value. Fracture coalescence can be observed, with in-field nucleation sites favoring regions of low Gc. Over time, the new nucleation site will propagate both forwards and backwards, eventually connecting to a primary fracture.

Spatial Variation of Elastic Stiffness

For a microstructure with spatial variation of elastic stiffness, regions of higher elastic stiffness will generate greater pressure due to the material’s higher resistance to volumetric compression. This leads to a pressure gradient which causes pressure diffusion into adjacent regions of lower pressure, resulting in higher strain energy concentrations at material interfaces, particularly where the elastic stiffness has high contrast. In the cases below, fracture nucleation occurs at the interface between stiffness regions in the stiffer material.

The Laplacian of the elastic stiffness field can identify regions of high concavity in elastic stiffness, signifying regions with the greatest dissimilarity in stiffnesses. Regions with a large negative magnitude of the Laplacian of stiffness exhibit higher accumulation of strain energy due to the significant localized pressure and pressure gradients. These locations act as nucleation sites for fracture.

Shock simulations were conducted at multiple impact velocities of 1.0 km/s, 0.4 km/s and 0.1 km/s to assess the influence of impact speed on the fracture behavior. As shown below, higher impact speeds lead to greater strain energy densities, resulting in early fracture nucleation, greater fracture branching, and greater accumulated damage. The nucleation behavior shows consistency with early findings, with higher stiffness regions with strong surrounding stiffness gradients acting as primary nucleation sites. After nucleation, the fracture propagation is less influenced by the stiffness field due to high strain energy availability.