Abstract
We used a physically motivated internal state variable plasticity/damage model containing a mathematical length scale to idealize the material response in finite element simulations of a large-scale boundary value problem. The problem consists of a moving striker colliding against a stationary hazmat tank car. The motivations are (1) to reproduce with high fidelity finite deformation and temperature histories, damage, and high rate phenomena that may arise during the impact accident and (2) to address the material postbifurcation regime pathological mesh size issues. We introduce the mathematical length scale in the model by adopting a nonlocal evolution equation for the damage, as suggested by Pijaudier-Cabot and Bazant in the context of concrete. We implement this evolution equation into existing finite element subroutines of the plasticity/failure model. The results of the simulations, carried out with the aid of Abaqus/Explicit finite element code, show that the material model, accounting for temperature histories and nonlocal damage effects, satisfactorily predicts the damage progression during the tank car impact accident and significantly reduces the pathological mesh size effects.
1. Introduction
The design of accident-resistant hazmat tank cars or the improvement of existing ones requires material models that describe the physical mechanisms that occur during the accident. In the case of impact accidents such as collisions, finite deformation, and temperature histories, damage and high rate phenomena are generated in the vicinity of the impact region. Unfortunately, the majority of material models used in the finite element (FE) simulation of hazmat tank car impact scenarios do not account for such physical features (unavoidably, this will under- or overestimate, for instance, the numerical prediction of the puncture resistance of hazmat tank cars’ structural integrity). Furthermore, in the few models that do, a mathematical length scale aimed at solving the postbifurcation problem is absent. As a consequence, when one material point fails in the course of the numerical simulations, the boundary value problem for such material models changes, from a hyperbolic to an elliptical system of differential equations in dynamic problems, and the reverse in statics. In both cases, the boundary value problem becomes ill posed, Muhlhaus [1], Tvergaard and Needleman [2], de Borst [3], and Ramaswamy and Aravas [4], as the boundary and initial conditions for one system of differential equations are not suitable for the other. Consequently, discontinuities in the strain and damage distribution appear, and all the strain and damage have the tendency to localize in a zone of vanishing volumes; this results in zero dissipated energy at failure of the material. Also, bifurcations with an infinite number of bifurcated branches appear, which raises the problem of selecting the relevant one, especially in numerical computations where this drawback manifests itself as a pathological sensitivity of the results to the FE discretization.
Alternatively to the shortcomings encountered in hazmat tank car impacts’ numerical simulations, we propose to use a nonlocal version of the Bammann-Chiesa-Johnson (BCJ) model, Bammann and Aifantis [5] and Bammann et al. [6], a physically motivated internal state variable plasticity and damage model containing a mathematical length scale. The use of internal state variables will enable the prediction of strain rate and temperature history effects; see Figure 1. These effects can be quite substantial and therefore difficult to incorporate into material models, which assumes that the stress is (1) a unique function of the current strain, strain rate, and temperature and (2) independent of the loading path. The effects of damage are included in the BCJ model, however, through a scalar internal state variable which tends to degrade the elastic moduli of the material as well as to concentrate the stress. The mathematical length scale is introduced in the model via the nonlocal damage approach of Pijaudier-Cabot and Bazant [7]. In the context of concrete damage, these authors have suggested a formulation in which only the damage variable is nonlocal, while the strain, the stress, and other variables retain their local definition. Their formulation has been applied to creep problems by Saanouni et al. [8] and extended to plasticity by, among others, Leblond et al. [9], Tvergaard and Needleman [10], and Enakoutsa et al. [11]. Following suggestion Pijaudier-Cabot and Bazant [7], a nonlocal evolution equation for the damage within an otherwise unmodified BCJ model is adopted in the current study. The time derivative of the damage is expressed as the spatial convolution of a “local damage rate” and bell-shaped weighting function. The width of this function introduces a mathematical length scale.
In this study, a dynamic nonlinear FE analysis, carried out with Abaqus FE code, is used to simulate a moving striker colliding against a stationary hazmat tank car. The structure part of this FE model is represented by Lagrangian elements obeying the nonlocal version of the BCJ model, while the fluid part is represented by Eulerian elements. The objective of this study is twofold: (1) use a high fidelity material model to idealize the physics occurring during the impact accident and (2) rectify the computational drawbacks (postbifurcation mesh dependence issues) for this model on a large-scale boundary value problem. The resulting numerical simulations of hazmat tank car impact scenarios which account for nonlocal damage and temperature history effects predict satisfactorily the tank car failure process independently of the element size. Although several other studies have addressed successfully the mesh dependence issues in finite element computations of problems involving dynamic failure of materials, as in for instance, Voyiadjis et al. [12] and Abu Al-Rub and Kim [13], very few of them deal with real-world large-scale boundary value problems. The present work investigates the application of the nonlocal integral approach to predict the failure of a large-scale vehicle under impact loading conditions. The paper is organized as follows. (i)The first section provides a summary of the equations of the BCJ model and its nonlocal extension. (ii)The second section discusses several methods to numerically implement the integral-type nonlocal damage into existing Abaqus FE BCJ model subroutines. The main difficulty encountered in this implementation relates to the double loop over the integration points required by the calculation of several convolution integrals, which might otherwise dismantle the architecture of the entire code. (iii)The last section is devoted to numerical applications of the local and nonlocal versions of the BCJ model on hazmat tank car shell impact accident simulations.
2. The BCJ Model: Local and Nonlocal Versions
2.1. The BCJ Model
The BCJ model is a physicallybased plasticity model coupled with the Cocks and Ashby [14]’s void growth failure model. The BCJ model incorporates load path, strain rate, and temperature history effects, as well as damage through the use of scalar and tensor internal state variables for which the evolution equations are motivated by dislocation mechanics and cast in a hardening minus recovery format. The BCJ model also accounts for deviatoric deformation resulting from the presence of dislocations and dilatational deformation.
The deformation gradient is multiplicatively decomposed into terms that account for the elastic, deviatoric inelastic, dilatational inelastic, and thermal inelastic parts of the motion. For linearized elasticity, the multiplicative decomposition of the deformation gradient results in an additive decomposition of the Eulerian strain rate into elastic, deviatoric inelastic, dilatational inelastic, and thermal inelastic parts. The constitutive equations of the model are written with respect to the intermediate (stress-free) configuration defined by the inelastic deformation such that the current configuration variables are corotated with the elastic spin. The pertinent equations of the BCJ model are expressed as the rate of change of the observable and internal state variables and consist of the following elements.
(i) A hypoelasticity law connecting the elastic strain rate to an objective time derivative Cauchy stress tensor is given by (strictly speaking, additional terms containing the temperature and its derivatives should be added to (1); the influence of these terms, however, is significant only for problems involving very high temperatures, such as welding problems; therefore, we can safely ignore them here) where is the Lam constant, is the shear modulus, and denotes the damage variable. The Cauchy stress is convected with the elastic spin as where, in general, for any arbitrary tensor variable , represents the convective derivative. Note that the rigid body rotation is included in the elastic spin; therefore, the constitutive model is expressed with respect to a set of directors whose directions are defined by the plastic deformation.
(ii) The decomposition of both the skew symmetric and symmetric parts of the velocity gradient into elastic and inelastic parts for the elastic stretching rate and the elastic spin in the absence of elastic-plastic couplings yields Note that for problems in the shock regime, only the deviatoric elastic strain part is linearized enabling prediction of large elastic volume changes.
(iii) Next, the equation for the plastic spin is introduced, in addition to the flow rules for and and the stretching rate due to the unconstrained thermal expansion . From the kinematics, the dilatational inelastic flow rule is given as Assuming isotropic thermal expansion, the unconstrained thermal stretching rate can be expressed by where is a linearized expansion coefficient.
For the plastic flow rule, a deviatoric flow rule [15] is assumed and defined by where is the temperature, the scalar hardening variable, the objective rate of change of , the tensor hardening variable, and the deviatoric Cauchy stress.
There are several choices for the form of . The assumption that allows recovery of the Jaumann stress rate. Alternatively, this function can be described by the Green-Naghdy rate of Cauchy stress. We used the Jaumann rate for the numerical applications in this paper.
(iv) The evolution equations for the kinematic and isotropic hardening internal state variables are given in a hardening minus recovery format by where and are the hardening moduli, and are scalar functions of describing the diffusion-controlled “static” or “thermal” recovery, and and are the functions of describing dynamic recovery.
(v) To describe the inelastic response, the BCJ model introduces nine functions which can be separated into three groups. The first three are the initial yield, the hardening, and the recovery functions, defined as The second group is related to the kinematic hardening process and consists of the following functions: The last group is related to the isotropic hardening process and is composed of In (8), (9), and (10), is some parameter of the model which needs to be determined.
(vi) The evolution equation for damage, credited to Cocks and Ashby [14], is given by Note that this void growth model displays a “” dependency on the triaxiality factor , as well as an additional parameter along with the initial value of the damage required to calculate damage growth.
(vii) The last equation to complete the description of the model is one that computes the temperature change during high strain rate deformations, such as those encountered in high rate impact loadings. For these problems, a nonconducting (adiabatic) temperature change following the assumption that 90% of the plastic work is dissipated as heat is assumed. Taylor and Quinney [16] were the first to measure the energy dissipation from mechanical work as being between 5% and 50% of the total work for various materials and strain levels. Therefore, the rate of the change of the temperature is assumed to follow where and represent the material density and a specific heat, respectively.
The empirical assumption in (12) has permitted nonisothermal solution by FE that is not fully coupled with the energy balance equation (see [6]). Note that the temperature rise will induce a profound effect on the constitutive behavior of the material. Specifically, the temperature increase will lead to thermal softening (adiabatic shear bands), and as a result shear instabilities may arise. The model is also suitable to predict mechanical softening through a gradual increase of the damage. It is well known that practical F applications of constitutive models involving softening, like the BCJ model, are strongly mesh-dependent. According to Rousselier [17], this problem can be obviated by putting a lower limit on the element size. However, this practice is not optimal theoretically. Another, more elaborated, solution consists of including a mathematical length scale in these constitutive models. The following section presents a technique to embed the BCJ model with just such a mathematical length scale.
2.2. Embedding a Length Scale in the BCJ Model
Following suggestion Pijaudier-Cabot and Bazant [7] in the context of concrete damage, we propose to delocalize the variable(s) responsible for softening. In the BCJ model, softening may arise from two mechanisms: a gradual increase of the damage (under isothermal conditions) or a temperature rise (in adiabatic conditions) followed by an increase of the damage. While temperature and damage parameters seem to govern softening in adiabatic conditions, review of the model’s constitutive equations provided in the previous section reveals that these two variables are related. We choose to introduce the length scale on the damage evolution equation. This choice appears quite appealing from the physical point of view. Indeed, in the case of heterogeneous materials, for example, the damage can only be defined by considering “elementary” volumes of size greater than the voids spacing (the Cocks and Ashby [14] void growth model is based on a cylinder containing a spherical void) and is therefore a nonlocal quantity.
The evolution equation of this variable is given by a convolution integral including a bell-weighting function the width of which introduces a mathematical length scale: In this equation, denotes the volume studied, and the bell weight function defined as where is the mathematical length scale (note that the physical length scale available in the model due to rate sensitivity (see [18] for more details) is so small that it is impossible to solve a boundary value problem, since this will require finite elements of size smaller than that physical length scale, which is already very small). The factor and the “local damage rate ” are given by and (11), respectively. The function is indefinitely differentiable and does not introduce any Dirac’s -distribution at the point . This means that the function is not partially local but entirely nonlocal. The function is also isotropic and normalized. The point here is that must be equal to if the latter variable is spatially uniform. This would not be the case near the boundary of in the absence of the normalization factor . The presence of this term allows for the coincidence everywhere.
The new evolution equation for the damage, (13), along with the equation for the temperature rise, (12), should predict satisfactorily the dynamic failure process of the tank car impact independently of the mesh size. Indeed, the local damage rate (13) implicitly depends on the temperature rise through the strain rate sensitivity parameter (see (6) and [19]). When the temperature gradually rises, which is the case in high-velocity impact loadings, the damage rate in this zone quickly climbs to a high value; as a result, the damage grows rapidly to reach the critical failure damage. The convolution integral in (13) enhances the rapid damage increase, since it involves the sum of several positive terms each of which contains a local damage velocity.
Nonetheless, the numerical implementation of the new evolution equation for the damage rate into an existing FE code is not an easy task because of the double loop over integration points required by the calculation of several convolution integrals, which may potentially compromise the entire architecture of the existing code. The following section is devoted to explaining this implementation.
3. Numerical Treatment of Nonlocal Damage Rate
The numerical implementation of the original BCJ model into an FE code such as Abaqus has been extensively addressed in Bammann et al. [6]; consequently, it will not be repeated here. Recall, however, that this implementation is based on radial return method of R. D. Krieg and D. B. Krieg [20] to solve numerically the equations presented in Section 2.1 for the deviatoric stresses, equivalent plastic strain, pressure, temperature, and damage at each new time step. The algorithm is available in both the implicit and explicit versions of ABAQUS.
To implement the nonlocal extension of the BCJ model in the implicit version of the code, we have computed the nonlocal damage increment at elastic-plastic convergence by means of the built-in subroutine URDFIL, which stores all the variables necessary for the convolution operation performs this operation. This method involves all the integration points of the finite element model considered. However, since the Gaussian influence weighting function vanishes for integrations points that are outside the region limited by the characteristic length scale, only the integration points inside that zone matter in the calculation of the convolution integral. The nonlocal damage increment are stored for all the integration points involved in the FE model and are used to calculate the damage at time for the next time-step following the formula: This updated value is an explicit estimation of the damage at and is not used to repeat the whole process of solution between times and . Thus, the algorithm is not fully implicit, but mixed implicit/explicit: implicit with respect to all parameters (displacements, stresses, hardening parameters, etc.) except the damage parameter, this preserves the compatibility with the implicit schemes used for usual nonsoftening models. Enakoutsa et al. [11] used the same numerical technique to implement a nonlocal version of Gurson [21] model following a slightly modified Aravas [22] algorithm, the so-called “projection into the yield surface” algorithm. The explicit nature of Enakoutsa et al.’s [11] algorithm with respect to the damage allowed these authors to prove that the projection problem associated with Gurson’s nonlocal model has a unique solution. The latter property is a direct consequence of the fact that the constitutive equations of Gurson’s nonlocal model belong to the class of generalized standard materials as defined by Halphen and Nguyen [23]. We apply this numerical treatment of nonlocal damage rate on a double-notched edge plane-strain specimen under uniaxial tension problem to illustrate the validity of the method to avoid ill-posed issues in this laboratory-oriented boundary value problem. The tests have been run using two meshes with element sizes 1.5 and 3 mm in the notched region of the specimen, respectively. The value of the characteristic length scale was 9 mm in the two cases considered; it can be considered as a parameter that determined the region in which the interactions between neighboring integration points are nonzero. In this work, the length scale introduced in the model is not a physical length scale, but a mathematical length scale aimed at eradicating the pathological mesh size dependence issues predicted by the local model. Hence, the only requirement it must satisfy is to be greater than the finite element size considered so that interactions of neighboring points can be accounted for. The specimen was made of 304L stainless steel, and the material parameters of which are provided in Horstemeyer et al. [24]. The major remark that must be made about the simulation results for the plate problem is that while the local BCJ model concentrates the damage within a layer of width of one element size (see Figure 2), the presence of the delocalization term spreads the damage region to a zone beyond this layer and allows for the elimination of the mesh size effects see Figure 3.
(a)
(b)
(a)
(b)
To assess the validity of the method to regularize hazmat tank car impact boundary value problems requires the implementation of the nonlocal damage rate in the explicit version of the BCJ model subroutine. One option of doing so consists of branching outside the “nblock loop” in the subroutine VUMAT, which allows the computation of the nonlocal damage rate at each integration point using the coordinates, the local damage velocities, and the weights of “nblock” integration points. This method also successfully avoids the localization problems arising in the numerical simulations of double-notched edge specimen tensile tests and therefore is excepted to reduce, if not completely remove, the mesh sensitivity issues arising in the hazmat tank car impacts numerical simulations.
4. Hazmat Tank Car Shell Impact FE Simulations
The goal of the hazmat tank car shell impact FE simulations section is to predict numerically the hazmat tank car structural failure process independently of the mesh size effects, as encountered in numerical computations based on material models involving softening, just like the BCJ material model. No comparison with experimental results is needed for this purpose. We describe later the impact accident physical problem and the associated FE model in Abaqus, and we present different numerical predictions of the damage created by the impact.
4.1. Problem Definition
The problem considered here is schematized in Figure 5; it consists of a ram car weighing 286,000 pounds and a protruding beam to which an impactor was attached, moving horizontally at 10 m/s into an immobile hazmat tank car. Ninety percent of the tank car is filled with water mixed with clay slurry, which together has the approximate density of liquid chlorine. The tank is subjected to 100 psi internal pressure. The impactor used in the problem has a inch impact face, see Figure 4.
4.2. FE Model of the Problem
The FE model of the problem consists of an Eulerian mesh representing the fluid in the tank car and a Lagrangian mesh idealizing the tank, the jacket, and the impactor. The impactor consists of R3D4 rigid elements, while the tank and the jacket are meshed using solid C3D8R elements. For the sake of simplicity, the space between the tank and the jacket is assumed to be empty. The two legs are rigid and are idealized with squared analytical surfaces on which reference points located below the jacket are assigned. Each reference point is kinematically coupled with a definite set of nodes on the jacket and the tank car. The contact algorithm available in Abaqus/Explicit FE code is used to account for all possible contacts between the different parts in the model.
The Eulerian mesh is based on the volume-of-fluid (VOF) method. The VOF method, widely used in computational fluid dynamics problems, tracks and locates the fluid-free surface; it belongs to the class of Eulerian methods that are characterized by either a stationary or moving mesh. The VOF method suitably captures the change of the fluid interface topology. In this method, the material in each element is tracked as it flows through the mesh using the element volume fraction (EVF), a unique parameter for each element and each material. The material parameters used in the simulations are determined from tension, compression, and torsion tests under different constant strain rate and temperatures. We used the material parameters provided in Horstemeyer et al. [24]. Although these constants are based on isothermal, constant strain rate tests, the model can be applied to arbitrary temperature and strain rate histories.
The stationary part of the problem consists of a tank car surrounded by a jacket. The dimensions of the tank car and the jacket are the same as in the works of Tang et al. [25] and Tang et al. [26]. The tank is a 0.777-inch-thick cylinder and is closed at its two ends with elliptical caps of aspect ratio 2. The tank body material consists of 304L stainless steel; the jacket is a 0.119-inch-thick made one of the same steel. A 4-inch-thick layer separates the jacket from the tank car. This layer is introduced to account for insulation and thermal protection between the tank car and the jacket. The entire assembly (tank/layer/jacket) is supported by two rigid legs, and this assembly is placed with one side against a rigid wall and the other side exposed to impact from the impactor.
4.3. Results and Discussions
This section presents evidence of the integral nonlocal damage method to reduce mesh dependence issues which arise in the FE solution of hazmat tank car impact accident problems. To that end, the calculations were performed on three different meshes, coarse, medium, and fine meshes (see Figure 6), for the local and nonlocal versions of the BCJ model in adiabatic conditions. Under isothermal conditions, the BCJ model difficultly predicts failure, Bammann et al. [6]; therefore, isothermal conditions are not considered here. The element sizes in the impact region are 75 mm, 50 mm, and 35 mm for the coarse, medium, and fine meshes, respectively. The nonlocal damage effects are introduced at each integration point by taking the convolution integral over “nblock” integration points. This introduces in the model some weak nonlocality defined by the interactions between each integration point and “nblock” other integration points not necessary located in the immediate vicinity of the latter. The length scale associated with such nonlocality can be considered of measure “nblock integration points.”
(a)
(b)
(c)
Due to the calculation of the integral nonlocal evolution equation of the damage, the CPU time required by the nonlocal calculations is higher than the case with the local model. It usually takes up to 12 hours to complete the simulation of the impact problem with nonlocal damage using the coarse mesh (59049 integration points for the structure part of the model). The medium (80373 integration points for the structure part of the model) and fine meshes (158037 integration points for the structure part) simulation times are longer, 15 hours for the medium mesh and almost 20 hours for the refine mesh. This long time period is also due to the presence of the fluid, which is modeled by the element volume fraction (EVF) technique available in ABAQUS.
Figure 7 illustrates the repartition of damage on the surface opposite to the impacted surface of the tank for all three meshes at the same time, in the case of the local BCJ model. The results of the simulations using the local BCJ model are mesh dependent; the pattern of the damage is determined by the size of the elements. Therefore, decreasing of FE size will modify the global response of the structure (which depends explicitly on the number of elements). Furthermore, the energy generated during the impact tends to zero when the size of the FE approaches zero. This leads to the meaningless conclusion that the tank car fails during the accident with zero energy dissipated. In fact, the FE solutions depend not only on the size of the elements, but also on their nature, orientation, degree of interpolation function, as well as on the FE approximation space, as presented in Darve et al. [27]. The mesh size effects encountered in the tank car FE simulations with the local BCJ model are more emphasized in Figure 8 which illustrates the damage pattern on the outer surface of the tank car. In this specific case, the damage pattern recovers the whole deformation region of the tank car; this zone shrinks with reducing mesh size, which is not realistic from the physical point of view.
(a)
(b)
(c)
(a)
(b)
Figure 9 is the analogous of Figure 7. Here, the pattern of the damage for the three meshes is relatively similar; the same observation holds for a latter time. The similarity would be more significant if a very refined mesh was used, because as in the case of local BCJ model, the element size determines the damage pattern. Thus, the nonlocal damage rate, and a fortiori the damage itself, has a significant influence on the response of the hazmat tank car to the impact loading.
(a)
(b)
(c)
Also remarkable from Figure 9 are the bands, which do appear in the local simulations; they seem to have almost completely disappeared. In fact, the nonlocal damage has delayed their apparition to latter times where these bands are extended over several elements. This postponement arises because the damage in the bands is much larger at latter stages of the impact loading, due to the progressive development of considerable stress in this zone. Figure 10, analogous of Figure 8, presents the damage pattern obtained with the nonlocal BCJ model on the outer surface of the tank car body in the impact region for two different mesh sizes. Like in Figure 8, the damage patterns predicted with the nonlocal BCJ model for the two mesh sizes considered are quasi similar.
(a)
(b)
5. Summary
With the aid of Abaqus/Explicit FE code, we have presented a dynamic nonlinear FE simulation of a hazmat tank car shell impact. In these simulations, the tank car body material was idealized with a physically motivated internal state variable model containing a mathematical length scale, the nonlocal BCJ model. This model consists of all the constitutive equations of the original BCJ model, except for the evolution equation for the damage, which is modified into a nonlocal one. Numerical methods to implement this equation in existing BCJ model subroutines were also presented.
The results of the simulations show the ability of the damage delocalization technique to significantly reduce the pathological mesh dependency issues while predicting satisfactory hazmat tank car impact failure, provided that nonlocal damage effects are coupled with temperature history effects (adiabatic effects). However, the performance of our approach will fully demonstrate its power when the damage progression it predicts reproduces the real-world hazmat tank car impact damage progression. This comparison with experimental results will be investigated in the near future.
Acknowledgments
The study was supported by the US Department of Transportation, Office of the Secretary, Grant no. DTOS59-08-G-00103. Dr. Esteban Marin and Mrs. Clémence Bouvard are gratefully thanked for their insightful discussions on the topic.