### 22.6.1. Equation of State

A general material model requires equations that relate stress to deformation and internal energy (or temperature).

In most cases, the stress tensor may be separated into a uniform hydrostatic pressure (all three normal stresses equal) and a stress deviatoric tensor associated with the resistance of the material to shear distortion (as is the case for most materials in Autodyn).

Then the relation between the hydrostatic pressure, the local density (or specific volume) and local specific energy (or temperature) is known as an equation of state.

The equation of state can be determined from knowledge of the thermodynamic properties of the material and ideally should not require dynamic data to build up the relationship. However, in practice, the only practical way of obtaining data on the behavior of the material at high strain rates is to carry out well-characterized dynamic experiments. It is important to recognize that, since the relationship is required for use in a numerical code, an analytic form is chosen to facilitate solution. Such an analytic form is at best an approximation to the true relationship. Further, the equation of state may be given in extensive tabular form and in that case the analytic form chosen can be considered as an interpolation relationship.

This section describes the Equations of State available in Autodyn that are not described in the Equations of State in the ANSYS Mechanical User's Guide.

#### 22.6.1.1. Tillotson EOS

This form of equation of state (Tillotson 1962 [1], Allen 1967 [2]) was derived to provide an accurate description of the material behavior of metallic elements over the very wide ranges of pressure and density met in hypervelocity phenomena.

Not only must such an equation of state describe normal density material and its condition after shock, but also its expansion and change of phase in cases where the shock energy has been sufficient to melt or vaporize the material.

The pressure range can be so large that the "low pressure" regime of this form of equation of state is defined as from 0 to 10 Mbar and "high pressure" from 10 to about 1000 Mbar. Thus any pressure and results from normal laboratory experiments cover only the "low pressure" region.

For the derivation of an equation of state for the "high pressure" region, analytic forms provide best fit interpolations between Thomas-Fermi-Dirac data at high pressures (above 50 Mbar) and experimental data at low pressures.

The formulation is claimed to be accurate to within 5% of the Hugoniot pressure and to within 10% of the isentropic pressures. It is therefore a very useful form of equation of state to use for general hypervelocity impact problems.

##### 22.6.1.1.1. Regions of Interest

The total range of the pressure-volume plane is divided into four regions as shown in the following figure.

The region to the left of the Hugoniot can only be reached by adiabatic (non-shock) compression and is not relevant to impact problems. It is therefore excluded from the formulation.

Region I represents the compressed phase of the material and extends vertically to shock pressures of about 1000 Mbar.

Region II describes material which has been shocked to an energy less than the sublimation energy es and will therefore, on adiabatic release, return to zero pressure as a solid.

There is no provision in this equation of state to describe the material at pressures less than zero.

Region IV is the expansion phase of material which has been shocked to an energy es' sufficiently large to ensure that it will expand as a gas at very large expansions.

For convenience the change of phase is assumed to take place at v = v0.

For large specific volumes, the formulation for Region IV extrapolates to an ideal gas limit.

It is desirable, or even necessary, to ensure that the formulations in each region provide continuous values of the pressure and its first derivatives at the boundaries between contiguous regions.

This is true for the boundary between Regions I and II at v = v0 but if es' is set equal to es difficulties may arise at the boundary between Regions II and IV at volumes slightly greater than v0.

To avoid this problem Region III is defined by setting

 (22–1)

where ev is the vaporization energy at zero pressure determined, like the boiling point energy es, from thermodynamic data.

Region III lies between Regions II and IV where es < e < es'. In this region the pressure is calculated as a mean between that calculated for Regions II and IV.

The constant k is determined empirically to ensure a well behaved solution and in practice satisfactory solutions have been obtained with values of k in the range of zero to one fifth.

##### 22.6.1.1.2. Regional Formulations

In what follows let:

 (22–2)

You are prompted for the reference density ρref (or ρ0 in the equations below) and the constants in the Tillotson equation of state viz. a, b, A, B, α, β, e0, es, and es'.

For Region I (μ 0) the pressure p1 is given by a Mie-Gruneisen equation of state, but since the formulation is to be valid for a very large range of pressure, the Gruneisen Gamma is a function of both v and e, not just a function of v alone.

The constants fit the low pressure shock data but they are adjusted to fit the asymptotic Thomas-Fermi behavior for the variation of pressure at maximum compressions (like a monatomic gas).

The formulation for Region II is as for Region I with a slight modification to one term to avoid difficulties as μ goes increasingly negative.

In Region IV the formulation is chosen to give the correct behavior both at high pressure/normal density and for very large expansion (where it must converge to an ideal gas behavior).

With these constraints the different formulations are:

Region I (μ 0)

The pressure p1 is given by

 (22–3)

Region II (μ < 0, e es)

The pressure p2 is given by the same formulation as that for Region I except that B = 0, i.e.

 (22–4)

Region III (μ < 0, es < e < es')

The pressure p3 is given by

 (22–5)

Region IV (μ < 0, e es')

The pressure p4 is given by

 (22–6)

where x = 1 - 1 / η.

#### 22.6.1.2. Puff EOS

Like the Tillotson EOS, this equation of state was constructed to cover the behavior of material from cold shocked regions to hot, highly expanded regions.

The equation of state was first used in the PUFF/KO hydrocodes (Brodie & Hormuth [3], 1966: Bakken & Anderson, 1969 [4]) and is based on a Mie-Gruneisen form but with a variable Gruneisen Gamma in the expanded phase to give the required convergence to perfect gas behavior at very large expansions if the energy e (ρ = ρ0) is greater than the sublimation energy es.

Like the Tillotson EOS, the (p, v) plane is divided into separate regions according to whether the material is compressed or expanded, and if expanded whether the energy is less than or greater than the sublimation energy.

##### 22.6.1.2.1. Regions of Interest

The total range of the pressure-volume plane is divided into three regions as shown in the following figure.

The region to the left of the Hugoniot can only be reached by adiabatic (non-shock) compression and is not relevant to impact problems. It is therefore excluded from the formulation.

Region I represents the compressed phase of the material and extends vertically to shock pressures of about 1000 Mbar.

Region II describes material which has been shocked to an energy less than the sublimation energy es and will therefore, on adiabatic release, return to zero pressure as a solid.

There is no provision in this equation of state to describe the material at pressures less than zero.

Region III is the expansion phase of material which has been shocked to an energy es sufficiently large to ensure that it will expand as a gas at very large expansions.

##### 22.6.1.2.2. Regional Formulations

You are prompted for a Gruneisen Gamma, Γ0, an expansion coefficient, H (the effective ideal gas constant at large expansions, γ - 1) the sublimation energy, es, and constants A1, A2, A3, T1, and T2.

Then the formulation is:

Region I (μ 0)

The pressure p1 is given by

 (22–7)

in this equation Γ ρ = Γ0 ρ0.

Region II (μ < 0, e < es)

The pressure p2 is given by

 (22–8)

If T1 is input as zero it is reset to T1 = A1

Region III (μ < 0, e es)

The pressure p3 is given by

 (22–9)

where N = A1 / ρ0 Γ es.

Unlike the Tillotson EOS, there is no interpolative region between regions II and III.

#### 22.6.1.3. Two-Phase EOS

This formulation (Morgan 1973 [5]) only considers material which has expanded from normal density but gives a realistic description of both the single phase (liquid or vapor) and the two-phase (in which liquid and vapor co-exist) behavior of the expanding material, provided that the conditions along the saturation curve for this material can be tabulated.

The pressure-volume (p, v) plane for v > v0 is divided into two regions by the saturation curve as shown in the following figure.

Above the saturation curve the material exists in a single phase state, either liquid or vapor, while below the saturation curve both phases coexist.

On the saturation curve the values of p, v, e, Γ, and temperature T are tabulated at suitable intervals to ensure accurate interpolations in the required calculations.

This equation of state should be used with some other form which will cover the region of compressed states (e.g. Shock or Polynomial). In this way the whole (p, v) plane will be covered.

##### 22.6.1.3.1. Single Phase Region

In the single phase region the behavior of the material is described by a Gruneisen equation of state of the form

 (22–10)

where in this case pr (v), er (v) and Γr (v) are the tabulated functions along the saturation curve and are determined for a specific value of v by interpolation within the tabular entries.

It is assumed that on the liquid side of the saturation curve Tr, er and Γr vary quadratically with vr along the curve and that pr may be obtained from a relationship of the form

 (22–11)

On the vapor side of the saturation curve er and Γr are taken to vary linearly with Tr, and vr and Tr are assumed to be related by

 (22–12)

and

 (22–13)
##### 22.6.1.3.2. Two-Phase Region

In the two-phase region below the saturation curve the specific volume v of the liquid-vapor mixture is written as

 (22–14)

and the internal energy e as

 (22–15)

where α is the mass fraction of vapor in the liquid-vapor mixture and vg, vl, eg and el are the specific volumes and internal energies of the saturated liquid and vapor at the saturated vapor pressure p and temperature T, i.e. they are conditions on the saturation curve.

Eliminating α from the two equations above gives

 (22–16)

This is an implicit equation defining p (or T) in terms of v and e.

Data on the saturation curve for water has been extracted from Morgan’s paper and included as a standard option in Autodyn.

For any other material you will have to supply the required saturation curve data using subroutine EXTAB (the comments in EXTAB explain how to use this subroutine).

Note:  Energy deposition using two phase EOS in Euler may cause instabilities.

#### 22.6.1.4. Sesame EOS (Tables)

The Sesame equation of state tables have been included in Autodyn in the form of a Sesame material library.

This is an extensive library of tables of thermodynamic properties developed and maintained by the Equation of State and Opacity Group of the Theoretical Division of Los Alamos National Laboratory in the USA.

The library currently contains data for over 200 materials including metals, minerals, polymers and mixtures. Most of the tables have data for very wide ranges of density and energy, and are typically used for applications where these wide ranges are required, for example when materials undergo phase changes.

The library format for this EOS is different and you cannot add to or modify the tables, only retrieve data.

In general, our experience suggests that the Sesame tables cannot be used as a black box. The tables contain data for specific ranges in pressure, density, temperature space for each material. Results should be interrogated thoroughly to ensure that simulated material states are consistent with those expected for a particular application.

If you want to use the Sesame tables in Autodyn, you must obtain the data base from LASL first.

The procedure below describes how to use the SESAME database from Los Alamos with Autodyn:

• Obtain the sesconv.exe program from ANSYS and run this to convert the ASCII SESAME file to a binary format suitable for use with Autodyn.

• Put the resulting seslib file in the application folder C:\Documents and Settings\<usr name>\Application Data\Ansys\v160\AUTODYN.

#### 22.6.1.5. Explosive EOS

The approach required to model the initiation and detonation of explosives and the expansion of the detonation products varies greatly depending on the type of explosive being modeled.

The numerical modeling of the detonation of High Explosives (HE) is well established and precise. In this case, the detonation wave travels at the local sound velocity, ensuring that all the available chemical energy is released before the explosive reacts with its environment. This means that the amount of energy released during the detonation process is fixed and not dependent on the environment (e.g. whether the HE is detonated in air or within a closed container).

This is not the case for other types of explosives (e.g. pyrotechnics, propellants, low density explosives). These explosives burn much slower and consequently, the amount of energy they release depends critically on their environment. This type of burning is referred to as deflagration rather than detonation.

Even with HE, the question sometimes arises as to whether or not the HE is initiated strongly enough for full detonation to ensue.

Autodyn has material models that allow you to simulate all these effects to some degree.

##### 22.6.1.5.1. Lee-Tarver Model

Using the JWL EOS with burn on time provides a good simulation of the detonation process if it is clear that the explosive does indeed fully detonate (release all its chemical energy).

Early experimental work on the initiation of detonation indicated that as a shock front proceeds through a heterogeneous explosive it interacts with local density discontinuities, producing numerous local hot spots that explode but do not propagate.

However energy is released which strengthens the shock so that when it interacts with further non-homogeneities, higher temperatures are produced at these hot spots and more of the explosive is decomposed.

There is thus a run distance for the build-up to detonation and it has been found (Ramsay & Popolato 1965 [6]) that each explosive has a characteristic run to detonation versus shock pressure—the so-called "Pop plots".

Making the further assumption that the explosive passes through the same pressure, distance and state of partial decomposition, regardless of the initial conditions, i.e. that it follows a unique path in distance, time and state space, Forest (Forest 1976 [7]) developed the Forest Fire model to be consistent with the Pop plots by suitable choice of arbitrary coefficients.

This so-called ignition and growth model was not, however, so successful in matching pressure-time data obtained from gauges embedded in the explosive.

A more recent ignition and growth model, the Lee-Tarver model is included in Autodyn.

This model, as with Forest Fire model, is based on the assumption that ignition starts at local hot spots and grows outward from these sites.

The initial version of this model (Lee & Tarver 1980 [8]) described a two-part reaction rate model with a term for ignition of the explosive and a term describing the growth. The reaction rate in this model was given by:

 (22–17)

where the compression μ is, as usual,

 (22–18)

F was the reaction ratio (ratio of the mass of the gaseous explosive to the total mass of the explosive), p was the pressure in the explosive and I, x, r, G, y, and z were constants. In this model r = 4, so the ignition term depended on the fourth power of the compression of the explosive.

In the more recent model included in Autodyn (Tarver, Hallquist & Erickson 1985 [9]), an additional term is added to overcome the deficiencies of the earlier model in matching data for very short shock pulse initiation.

In this new model, there is a three-fold process of ignition, growth and completion and the reaction rate is given by:

 (22–19)

where F, μ and p are defined as above and I, b, a, x, G1, c, d, y, G2, e, g and z are constants.

For the explosive PBX-9404, x = 20, so the ignition term depends upon the 20th power of the compression.

The model described by the above equation gives a very rapid pressure spike on ignition, followed by a slow growth of the reaction that accelerates when the regions around the hot spots coalesce.

In Autodyn, the Lee-Tarver equation of state is used to model both the detonation and expansion of high explosives.

The solid inert explosive is modeled using either a JWL EOS or a shock EOS, while the gaseous expansion products are modeled with the JWL EOS with the reaction rate given by the equation above.

For a more physical post-burn behavior for medium and large expansions, the strength of the material is dependent on the burn fraction, F:

 (22–20)

where σ is the current yield strength and σ0 is the original yield strength of the unburned material.

##### 22.6.1.5.2. Powder Burn EOS

In the past, the Powder Burn model has also been known in Autodyn as the Slow Burn material model.

Two different formulations of the powder burn model have been implemented, dependent on the material behavior of the reacted (burned) material:

• A JWL equation of state with constant burn front velocity

• An Exponential equation of state with a burn front velocity that is dependent on the pressure.

The powder burn model is available for the 2D and 3D structural Lagrange, SPH, and Euler multi-material solver types. For Euler, only the formulation using an Exponential reacted equation of state can be used.

There is a need to numerically simulate the combustion of materials where the dominant physical characteristic of the burning is deflagration rather than detonation of the combustible material. The complicated chemical reactions that dictate the burning of the energetic material are extremely difficult to model from first principles. Therefore, a mathematical model based on the averaged quantities of the governing properties has been proposed [10]. The model was developed using detailed studies of a variety of experiments involving energetic materials where the convective burning front moves behind the pressure front of a shockwave. This type of burning occurs in many incendiary devices and munitions.

Theory

The powder burn model is a multi-phase model where gas and solid can be present in one cell at the same time step. The total mass within each cell is found by added together the mass of the gas and of the solid. The volume taken by the gas and solid are both known and therefore the density, compression, etc., of the material within the cell can be calculated. The burning of the material is dependent on the burn fraction F(t) where,

 (22–21)

and

 (22–22)

The burning velocity, c(t), is given by a user supplied function, H, dependent on the pressure of the gas,

 (22–23)

It has been shown [10] that the burn fraction can be represented as a function of H and a function, G(F(t)), which expresses the size and shape of the burning surface, such that:

 (22–24)

The functions G and H are user supplied and are entered in the material data menu, and a typical burn fraction variation in time is shown below:

The burn fraction should only be calculated on the amount of solid material within each cell which has been ignited. Therefore the total burn fraction within each cell, Ft(t) is found by multiplying the burn fraction, F(t), multiplied by the ignition fraction, I(t), such that

 (22–25)

The ignition fraction, I(t) is dependent on the velocity of the burn front, V(t), which is also known as ignition velocity and which depends on the equation of state selected by the user for the reacted gaseous material.

When the JWL EOS is selected for the Reacted EOS the ignition velocity is taken constant

 (22–26)

and when the Exponential EOS is selected the ignition velocity depends on pressure and density according to:

 (22–27)

where C1, and C2 are user defined constants, H and γ are user defined functions and ρs is the solid density averaged over the cell volume.

The total burn fraction then dictates the amount of solid which is burnt and has changed into a gaseous state at each time step.

When an Exponential reacted EOS is chosen the pressure of the gas, Pg, is calculated as:

 (22–28)

where ρg is the gas density, eg the gas internal energy and D is a user defined constant.

The pressure of the solid is evaluated using a separate equation of state, which is user defined and can be either a linear or a compaction EOS. As the solid material burns it changes state such that when the cell is fully burned only gas is present within the cell. A typical density variation in time is shown below,

Material Data Input

Autodyn requires the user to input the values of the burn fraction growth parameters, G and c. These growth parameters are used to evaluate the change in the burn fraction where,

 (22–29)

When the JWL EOS is chosen for the reacted gaseous material the constant burn front velocity must be entered as well as the relationship between H and Pg. Because the burn velocity is constant, the ignition time of the elements needs to be defined using the on-time detonation logic.

When the Exponential EOS is chosen for the reacted gaseous materials the constants C1, C2, and the relationship between H and Pg must be entered in addition to the variation between density and γ in the relationship below:

 (22–30)

Because the burn velocity for this reacted EOS is not constant but varies depending on gas pressure and density, the ignition times for the elements cannot be defined using the on-time detonation logic, and to start the ignition of elements the option of manually setting the detonation (ignition) times needs to be used.

The unreacted solid can be represented using a compaction or linear EOS.

Example

 Example application of Powder Burn model used in the Multi-material solver to propel a sabot and projectile inside an experimental gun chamber

#### 22.6.1.6. Rigid

A rigid material is available under the Equation of State option. Regions of the model filled with such a material will act as a single rigid body, irrespective of whether the individual elements filled are connected to each other. Multiple rigid bodies can be represented by filling different regions of the model with different rigid body materials.

A rigid body defines an un-deformable object that has mass, inertia and momentum. If the mass and inertia parameters are set to zero in the material definition, these quantities will be derived automatically from the material density and the volume of the elements filled. Alternatively, the user can directly specify the body’s mass and inertia using the first 10 input parameters.

By default, a kinematic rigid body is defined in Autodyn and its motion will depend on the resultant forces and moments applied to it through interaction with other Parts of the model. Elements filled with rigid materials can interact with other regions via contact or Euler coupling.

A constraint or imposed motion can be applied to a rigid body by selecting the “Constrain rigid body option” in the material definition. Alternatively, a boundary condition can be applied to all the nodes belonging to the rigid body.

Note:  Parts containing any of the new unstructured solvers can be filled with rigid materials. Rigid materials cannot be used in the structured solvers.