- Email: [email protected]
A High-Order Eulerian Godunov Method for Elastic–Plastic Flow in Solids
Journal of Computational Physics 167, 131–176 (2001) doi:10.1006/jcph.2000.6665, available online at http://www.idealibrary.com on
A High-Order Eulerian Godunov Method for Elastic–Plastic Flow in Solids1 G. H. Miller†, ‡ and P. Colella† †Applied Numerical Algorithms Group, Lawrence Berkeley National Laboratory, Berkeley, California 94720; and ‡Department of Applied Science, University of California, Davis, One Shields Avenue, Davis, California 95616 E-mail: [email protected]; [email protected] Received May 2, 2000; revised November 10, 2000
We present an explicit second-order-accurate Godunov finite difference method for the solution of the equations of solid mechanics in one, two, and three spatial dimensions. The solid mechanics equations are solved in nonconservation form, with the novel application of a diffusion-like correction to enforce the gauge condition that the deformation tensor be the gradient of a vector. Physically conserved flow variables (e.g., mass, momentum, and energy) are strictly conserved; only the deformation gradient field is not. Verification examples demonstrate the accurate capturing of plastic and elastic shock waves across approximately five computational cells. 2D and 3D results are obtained without spatial operator splitting. °c 2001 Academic Press Key Words: solid mechanics; shock waves; Godunov method; elasticity; plasticity.
1. INTRODUCTION
In this work, we present a high-order Godunov method for computing in Eulerian coordinates the multidimensional dynamics of elastic–plastic solids undergoing large deformations. Our approach is based on a new formulation of the equations of solid mechanics as a first-order system of hyperbolic partial differential equations (PDEs), a modification of that used by Trangenstein and Colella [25]. In [25], the usual conservation laws for mass, momentum, and energy, plus a constitutive model, are augmented by a form of equality of mixed partial derivatives that yields conservation equations for the entries of the inverse 1
Work at the Lawrence Berkeley National Laboratory was sponsored by the U.S. Department of Energy (DOE) Mathematical, Information, and Computing Sciences Division under Contract DE-AC03-76SF00098. Other work was supported by a subcontract from the Caltech Center for the Simulation of Dynamic Response in Materials, which in turn is supported by the Academic Strategic Alliances Program of the Accelerated Strategic Computing Initiative (ASCI/ASAP) under Subcontract B341492 of DOE Contract W-7405-ENG-48. 131 0021-9991/01 $35.00 c 2001 by Academic Press Copyright ° All rights of reproduction in any form reserved.
132
MILLER AND COLELLA
deformation gradient. This leads to equations of the form ∂U + ∇ · F(U ) = S(U ). ∂t
(1)
Here S(U ) contains source terms associated with the treatment of plasticity. These equations by themselves are not sufficient to specify the problem. In addition, we must impose linear constraints on the solution to guarantee that the inverse of the deformation gradient is in fact a gradient, i.e., that the curl of the rows of the deformation gradient vanish. These constraints can be written in the form L C (U ) = 0.
(2)
Here L C is a system of linear differential operators with constant coefficients. The constraint equation is an initial-value constraint: if (1) is satisfied, and L C (U ) is identically zero at some initial time, then L C (U ) vanishes identically for all later times. The constraint (2) plays an essential role in the analysis of the characteristic structure of the system (1). To get the physically correct eigenvectors and eigenvalues from the quasilinear form of the equations, one must use the constraint to replace some of the spatial derivatives. In general, solutions to (1), without imposing (2), give rise to unphysical wave propagation properties, even for linearized waves as was observed in [25]. A difficulty arises when one attempts to compute solutions to (1) and (2) using a conservative finite difference method. To the extent that a modified equation analysis is valid, we expect the behavior of the numerical solution to behave very similarly to the solution to the following system of PDEs: ∂U Mod + ∇ · F(U Mod ) = S(U Mod ) + τU (U Mod ) ∂t L C (U Mod ) = τC (U Mod ).
(3)
Here τU and τC are truncation error terms, which are nonzero. In general, these terms, and in particular τC , cannot be eliminated. The practice of enforcing a discretized form of the constraint (2) at the end of each time step using a Hodge projection would guarantee that a discretized form of (2) is satisfied identically. However, that will change the form of τC , but not set it to zero. The observation that the truncation error terms are a small perturbation to the equations is not sufficient to guarantee that U Mod is close to U . There is much less known about the well-posedness of systems of equations that are combinations of evolution equations and constraints than there is about pure evolution equations, and unexpected pathologies are known to occur [17]. The approach we want to take on this problem starts with an analysis due to Godunov [8, 9]. Numerical methods based on this approach have been recently investigated for the MHD equations in [20], the case for which Godunov first applied this analysis. Godunov modifies (1) in the following way: ∂U + ∇ · F(U ) = S(U ) + ξ L C (U ). ∂t
(4)
Here ξ = ξ(U ) can be chosen so that the system has the physically correct linearized eigenstructure, independent of whether L C vanishes. In addition, L C satisfies a transport
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
133
equation such that if L C (U ) is identically zero at some time, then it remains so for all later times. The numerical method we present here is based on the form of the equations given by (4). Thus we are discretizing a well-posed initial value problem without constraints, independent of whether the constraint (2) is satisfied. This gives us a high degree of confidence that a stable and consistent method can be developed. Of course, the extent to which we compute a solution to the original physical problem (1) and (2) depends strongly on whether the constraint is satisfied, but now that is purely an accuracy issue, without any impact on the stability of the method. In fact, we will investigate the use of various methods of limiting discrete measures of L C , similar to filtering methods developed for incompressible flow [11, 21]. Examples that demonstrate the method for elastic–plastic deformations of a homogeneous solid domain are presented. Aside from some simple one-dimensional problems that involve free surfaces, methods for handling material interfaces (contact discontinuities) are not described here. To address more general engineering problems, including those concerned with fluid–solid coupling, our intent is to combine this solid mechanics method with a volume-of-fluid interface treatment analogous to [15]. 2. GOVERNING EQUATIONS
The mechanical behavior of solids is described by observable variables (e.g., density ρ, momentum ρv, internal energy E, and the deformation F with respect to a chosen reference state) and also unobservable internal parameters which describe the response of the material to deviatoric stress. One constitutive representation of this behavior is through the multiplicative decomposition of the total deformation into elastic and inelastic components [12], F = F eF p.
(5)
Here F is the Lagrangian coordinate deformation which relates the spatial coordinate frame x = x(a, t) to the material coordinate frame a: Fαβ =
∂ xα . ∂aβ
(6)
We refer to F p as the plastic deformation tensor, although the numerical scheme we will present applies to more general inelastic deformations. According to (5), F p is a fictitious state of total deformation in which there is no elastic deformation: given an initial total deformation F, and a purely elastic relaxation path F e → I , the total observable deformation will evolve to F p , F → F p . The state F p is a function of the deformation history of the material. We represent this history through a single scalar parameter κ, a work-hardening measure, and constitutive flow rules F˙ p = h(ρ, g, F p , E, κ) κ˙ = K(ρ, g, F p , E, κ), which depend on the state variables but not their gradients.
(7) (8)
134
MILLER AND COLELLA
The equations of solid mechanics are then given by 0 ρvα ρ ρf ρvvα − σ eα ρv ρ(8 + v· f) ρ Ev − v σ α β βα ρE (v × (∇ × g T ))T e ge x gvδxα x T T (v × (∇ × g )) e ge gvδ yα y y ∂ ∂ = , gez + T T ∂t gvδz α (v × (∇ × g )) ez ∂ xα p ρF ex ρF p ex vα ρhex ρF p e p y ρF e v y α ρhe y p p ρF ez ρF ez vα ρhez ρκ
ρκvα
(9)
ρK
where ex , e y , and ez are the Cartesian unit vectors, and E is the sum of internal energy and kinetic energy (E = E + 12 v · v). For generality, we include a heat source term 8 and a body force vector f . The system of equations (9) is abbreviated ∂U ∂ Fα (U ) = S(U ), + ∂t ∂ xα
(10)
where U is the vector of quasi-conservation-form variables (ρ, ρv, ρ E, etc. . . .), Fα (U ) is the flux in direction eα , and where S(U ) is the vector of source terms. Here we follow the treatment in [25] and use the inverse deformation gradient g = F −1 as dependent variables. However, we introduce additional nonconservative terms in the evolution equations of g. We will show below that the addition of these terms leads to a well-behaved hyperbolic structure for the equations, independent of whether the curl of g T vanishes. However, we note here that G = ∇ × g T satisfies the following evolution equation: ∂G + ∇ · (vG − Gv) = 0. ∂t
(11)
In particular, if G vanishes identically at time t = 0, it vanishes at all later times. To solve these equations we adopt a predictor–corrector strategy. For each time step, we first solve the conservative flux differencing left-hand side of (10) using fluxes derived (by solution to Riemann problems) from edge- and time-centered variables that include timecentered contributions from the source terms. The solution obtained by flux differencing is then modified by addition of the source terms, evaluated using time-centered and cellcentered variables, and acting over the full time step 1t. The solution to the flux differencing equations is based upon the standard high-order Godunov strategy. This strategy begins with a characteristic analysis of the equations, which makes use of the linearized 1D equations in direction eα , ρ
v E gex ge y ∂ gez + A ∂ ∂t F p ex ∂ xα p F ey p F ez κ σ eα
ρ
0
v f E 8 gex 0 ge y 0 ge z = 0 , p he F ex x p he F ey y p F ez hez κ σ eα
K bα
(12)
135
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
where
vα
0 0 0 0 A= 0 0 0 0 0 0
ρeαT
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
−(σ eα )T /ρ vα
0
0
0
0
0
0
0
gδxα
0
vα I
0
0
0
0
0
0
gδ yα
0
0
vα I
0
0
0
0
0
gδzα
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα
−Aαα
0
0
0
0
0
0
0
0
0
−I /ρ 0 0 0 0 (13) 0 0 0 0 vα I
with Aαβ = −
∂σ eα g ∂geβ
(14)
and bα =
∂σ eα ∂σ eα ∂σ eα :h+ K+ 8. ∂F p ∂κ ∂E
(15)
The eigenvalue decomposition of A uses the technique of eigenvalue deflation and hinges upon recognition of the matrices Aαα as being acoustic wave propagation tensors for waves traveling in direction eα , ρ u¨ β = (Aαα )γ δ
∂ 2uδ , ∂ xβ ∂ xγ
(16)
where u is the displacement vector. The matrices Aαα are positive definite as a requirement of thermodynamic stability. This is made clear by writing Aαα in terms of gradients of the spatial displacements uˆ defined relative to the current configuration, (Aαα )βγ = ρ
∂ 2E . ∂ uˆ βα ∂ uˆ γ α
(17)
Here, uˆ βα is related to the deformation tensor Fβα with the reference coordinate frame {a} chosen to correspond to the current spatial frame {x}: uˆ βα = Fβα |{a}={x} − δβα .
(18)
Aαα is therefore a component of the Hessian of E, which is positive definite for a
136
MILLER AND COLELLA
thermodynamically stable material, and consequently Aαα has positive real eigenvalues and three linearly independent eigenvectors. Recognizing Aαα as being the acoustic wave propagation tensor suggests the wave equation solution 2 , Aαα X ac = ρ X ac 3ac
(19)
where 3ac is the diagonal matrix of acoustic wave speeds c, 3ac = diag(c1 , c2 , c3 ), and X ac are the acoustic displacement vectors. The linearized 1D matrix A then has eigenvalue decomposition A = X 3X −1 ,
(20)
with X , the matrix of right eigenvectors, given by
1 0 0 0 0 0 0 0 0
0 0 0 0 X = 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 I
0 0 0 0 0 0
0 0 I
0 0 0 0 0
0 0 0 I
0 0 0 0
0 0 0 0 I
0 0 0
0 0 0 0 0 I
0 0
0 0 0 0 0 0 I
0
0 0 0 0 0 0 0 1
−ρeαT X ac
−ρeαT X ac
−X ac 3ac T T (σ eα ) X ac /ρ (σ eα ) X ac /ρ −g X ac δxα −g X ac δxα −g X ac δ yα −g X ac δ yα −g X ac δzα −g X ac δzα , 0 0 0 0 0 0 0 0 X ac 3ac
2 ρ X ac 3ac
0 0 0 0 0 0 0 0 0
(21)
2 X ac 3ac ρ
and 3, the diagonal matrix of eigenvalues, given by
0 vα 0 0 vα 0 0 0 vα I 0 0 0 0 0 0 3= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 vα I 0 0 0 0 0 0 0
0 0 0 0 vα I 0 0 0 0 0 0
0 0 0 0 0 vα I 0 0 0 0 0
0 0 0 0 0 0 vα I 0 0 0 0
0 0 0 0 0 0 0 vα I 0 0 0
0 0 0 0 0 0 0 0 vα 0 vα I 0
0 0 0 0 0 0 0 0 0 − 3ac 0 vα I
0 0 0 0 0 0 . (22) 0 0 0 0 + 3ac
The wave speeds are Galilean invariant and properly analogous to the Lagrangian
137
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
representation, with 3 − waves with velocities vα − cγ , 3 + waves with velocities vα + cγ , and 21 material waves with speeds vα . X −1 , the inverse of X , is given by
1
0
0 0 0 0 0 0 0 0 −1 X = 0 0 0 0 0 0 0 0 1 −1 −1 0 2 3ac X ac −1 −1 0 − 12 3ac X ac
0 0 0 0 0 0 0 0
−2 −1 eαT X ac 3ac X ac
−2 −1 1 0 0 0 0 0 0 0 −(σ eα )T X ac 3ac X ac /ρ 2 −2 −1 0 I 0 0 0 0 0 0 g X ac 3ac X ac δxα /ρ −2 −1 X ac δ yα /ρ 0 0 I 0 0 0 0 0 g X ac 3ac −2 −1 0 0 0 I 0 0 0 0 g X ac 3ac X ac δzα /ρ . (23) 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 1 −2 −1 0 0 0 0 0 0 0 0 3 X 2ρ ac ac 1 −2 −1 0 0 0 0 0 0 0 0 3 X 2ρ ac ac 3. NUMERICAL METHOD: 1D
In 1D we discretize space into cells, indexed with subscript i, with width 1xi . Time is discretized in steps of 1t with integer superscript index n; t n+1 − t n = 1t. The generalization to 2D and 3D is similar, with indices j and k used for the second and third dimensions, respectively. Half-integral subscript indices represent edge-centered quantities. Lowercase Greek subscripts are used to denote vector and tensor indices. We begin by evaluation of the equation of state in each cell to determine the Cauchy stress σ , the acoustic wave propagation tensor Aαα , and the thermodynamic derivatives ∂σ/∂E|g,F p ,κ , ∂σ/∂g|E,F p ,κ , ∂σ/∂F p |E,g,κ , and ∂σ/∂κ|E,g,F p . Next, we evaluate the 1D slopes dq/d xα of the 27 primitive cell-centered variables q, q = (ρ, v, E, g, F p , κ, σ eα ).
(24)
We construct these slopes beginning with the van Leer slope in cell i which uses the monotonized limiter [26]: µ
∂q ∂x
¶vL i
µ
¶ 2|qi+1 − qi−1 | 2|qi − qi−1 | 2|qi+1 − qi | = sign(qi+1 − qi−1 ) min , , . 1xi−1 + 21xi + 1xi+1 1xi 1xi (25)
A fourth-order-accurate slope is then constructed as [5] µ
∂q ∂x
¶4th = i
³h ¡ ¢vL i h ¡ ¢vL i´ qi+1 − 14 1xi+1 ∂q − qi−1 + 14 1xi−1 ∂q ∂ x i+1 ∂ x i−1 1 1xi−1 4
+ 1xi + 14 1xi+1
.
(26)
138
MILLER AND COLELLA
To prevent overshoot and ringing, dissipation at strong shocks may be introduced via a “flattening parameter” χ, 0 ≤ χ ≤ 1, whence [5–7] µ
∂q ∂x
¶
µ = χi
i
∂q ∂x
¶4th .
(27)
i
The determination of this flattening parameter is described in a later section. These limited slopes are used to construct time-centered, edge-valued estimates of the primitive variables. The exact solution of the linearized equations, which we abbreviate as ∂q ∂q +A = s, ∂t ∂x
(28)
gives time-centered edge values n+1/2 q R,i−1/2 n+1/2
q L ,i+1/2
µ µ ¶ ¶ 1t 1xi 1t −1 ∂q Xi si = − 3i + I X i + 2 1xi ∂x i 2 µ µ ¶ ¶ 1t 1xi 1t −1 ∂q n = qi − 3i − I X i + Xi si . 2 1xi ∂x i 2 qin
(29a) (29b)
However, this construction uses both upwind and downwind characteristics. We make the method strictly upwind by filtering out the downwind characteristics: n+1/2 q R,i−1/2 n+1/2
q L ,i+1/2
µ µ ¶ ¶ 1t 1xi 1t −1 ∂q X i P− si = − 3i + I X i + 2 1xi ∂x i 2 µ µ ¶ ¶ 1t ∂q 1xi 1t X i P+ si = qin − 3i − I X i−1 + 2 1xi ∂x i 2 qin
(30a) (30b)
with projection operators P± defined as µ
µ P−
µ
µ P+
1t 3+ I 1x 1t 3− I 1x
(¡ 1t
¶¶ αβ
=
¶¶ αβ
1x
0 (¡ 1t
=
1x
¢ 3αα + 1 δαβ ¢ 3αα − 1 δαβ
3αα ≤ 0 3αα > 0 3αα ≥ 0 3αα < 0.
0
(31a)
(31b)
At each cell edge (i + 1/2), time-centered values are thus obtained from the left (i) and right (i + 1) neighboring cells. These edge values are then used to pose a Riemann problem: n+1/2 an initial value problem with constant left and right initial states given by q L ,i+1/2 and n+1/2 q R,i+1/2 , respectively. We approximate the solution to the Riemann problem by decomposing n+1/2 n+1/2 the jump q R,i+1/2 − q L ,i+1/2 in terms of the eigenvectors X of the linearized coefficients A. Specifically, n+1/2
n+1/2
q R,i+1/2 − q L ,i+1/2 =
27 X
ϕγ X γ ,i+1/2 ,
(32)
γ =1
where eigenvector column X γ ,i+1/2 is evaluated with certain L or cell-i properties if 3γ γ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
139
is a member of the − family (i.e., of the form veα − c), or with certain R or cell-(i + 1) properties if 3γ γ is a member of the + family (of the form veα + c), as given by the discretization
1 0 0 0 0 0 0 0 0
X i+1/2
0 0 0 0 = 0 0 0 0 0 0
−ρi eαT X ac,i
0 0 0 0 0 0 0 0
X ac,i 3ac,i
1 0 0 0 0 0 0 0
(σi+1/2,L eα )T X ac,i ρi
0
I
0 0 0 0 0 0 −gi+1/2,L X ac,i δxα
0 0
I
0 0 0 0 0 −gi+1/2,L X ac,i δ yα
0 0 0 0 0 0
0 0 0 0 0 0
I 0 0 0 0 0
0 0 0 0 0 0
0 I 0 0 0 0
0 0 I 0 0 0
0 0 0 I 0 0
0 0 0 0 1 0
−gi+1/2,L X ac,i δzα 0 0 0 0 2 ρi X ac,i 3ac,i
−ρi+1 eαT X ac,i+1
−X ac,i+1 3ac,i+1 (σi+1/2,R eα )T X ac,i+1 ρi+1 −gi+1/2,R X ac,i+1 δxα −gi+1/2,R X ac,i+1 δ yα . −gi+1/2,R X ac,i+1 δzα 0 0 0 0 2 X ac,i+1 3ac,i+1 ρi+1 (33)
In this expression, the density ρ, and the components (X ac , 3ac ) of the acoustic propagation tensor, are evaluated at the cell centers to avoid multiple evaluations of the equation of state. From the coefficients ϕγ of the jump decomposition, the material velocity v ∗ · eα at the cell edge is determined by adding to the L state the contributions of the − family, or by subtracting from the R state the contributions of the + family, that is, ∗ · eα = v ∗L ≡ vi+1/2,L · eα + ϕ6 X 6β,i+1/2 + ϕ7 X 7β,i+1/2 + ϕ8 X 8β,i+1/2 vi+1/2
(34)
or ∗ vi+1/2 · eα = v ∗R ≡ vi+1/2,R · eα − ϕ9 X 9β,i+1/2 − ϕ10 X 10β,i+1/2 − ϕ11 X 11β,i+1/2 , (35)
where β = 2, 3, 4 for directions eα equal to ex , e y , ez , respectively. We average the results of these calculations to determine the normal-direction edge velocity v ∗ · eα , ∗ · eα = vi+1/2
1 ∗L (v + v ∗R ). 2
(36)
For other properties to be evaluated at the cell edge as solutions to the Riemann problem we do not average the values evaluated from the L and the R states as above. Instead, we evaluate from the L state if v ∗ · eα is positive, or from the R state if v ∗ · eα is negative. Only if v ∗ · eα is approximately zero do we average these estimates. The evaluations include only
140
MILLER AND COLELLA
upwind characteristics by writing
∗ qi+1/2
P qi+1/2,L + γ w L ,γ ϕγ X γ P = qi+1/2,R + γ w R,γ ϕγ X γ ¢ P 1 ¡q i+1/2,L + qi+1/2,R + γ (w L ,γ + w R,γ )ϕγ X γ 2
∗ vi+1/2 · eα > ² ∗ vi+1/2 · eα < −² (37) ¯ ¯ ∗ ¯ ¯v i+1/2 · eα ≤ ²
with ( w L ,γ = ( w R,γ =
1
∗ 3γ γ ,i − vi · eα + vi+1/2 · eα < −²
0
otherwise
1
∗ · eα > ² 3γ γ ,i+1 − vi+1 · eα + vi+1/2
0
otherwise.
(38a)
(38b)
The value of w L ,γ is 1 when eigenvalue γ , estimated using the ∗ value of the material velocity together with the i cell-centered acoustic wave speeds, is negative and w L ,γ is 0 otherwise; w R,γ is 1 when the approximated value of eigenvalue γ is positive and 0 otherwise. In our computations presented below, we use a value of ² = 10−9 . By this procedure, we obtain the edge ∗ value solutions of the Riemann problem, ρ ∗ , v ∗ , E ∗ , g ∗ , F p∗ , κ ∗ , and (σ e j )∗ . These are then used to compute edge-valued fluxes (cf. (9, 10)). For example, in direction ex ,
Fx,i+1/2
ρvx
∗
ρvx2 − σx x ρv y vx − σ yx ρvz vx − σzx ρ Ev − v σ − v σ − v σ x x xx y yx z zx vx gex + v y ge y + vz gez = . 0 0 p ρvx F ex p ρvx F e y ρvx F p ez ρvx κ i+1/2
(39)
In 1D we obtain a preliminary update U˜ of the variables U by conservatively differencing the fluxes: ¢ 1t ¡ ∗ ∗ Fi+1/2 − Fi−1/2 . U˜ in+1 = Uin − 1xi
(40)
The final time-(n + 1) value of the variables U is obtained from the preliminary values U˜ by addition of the source terms S: Uin+1 = U˜ in+1 + 1t Si .
(41)
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
141
We discretize this in the general 3D (Cartesian) case as
0
n+1/2 n+1/2 ρi fi ¢ n+1/2 ¡ n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 ρi + vx,i f x,i + v y,i f y,i + vz,i f z,i 8 ∗ ∗ ∗ ∗ ¡ ¢ © gi, j+1/2,k ex − gi, j−1/2,k ex n+1/2 gi+1/2, j,k e y − gi−1/2, j,k e y − v y,i jk 1xi 1y j ∗ ∗ ∗ ∗ ¡ ¢ª g e − g e g e − g e n+1/2 i, j,k+1/2 x i, j,k−1/2 x i+1/2, j,k z i−1/2, j,k z − − vz,i jk 1z k 1xi ∗ ∗ ¢ © n+1/2 ¡ gi,∗ j+1/2,k ez − gi,∗ j−1/2,k ez gi, j,k+1/2 e y − gi, j,k−1/2 e y − v z,i jk 1y j 1z k ∗ ∗ ∗ ∗ ¡ ¢ª gi, j+1/2,k ex − gi, j−1/2,k ex n+1/2 gi+1/2, j,k e y − gi−1/2, j,k e y n+1 n+1 . ˜ − − v Ui = Ui + 1t x,i jk 1xi 1y j ∗ ∗ ¢ © n+1/2 ¡ gi,∗ j,k+1/2 ex − gi,∗ j,k−1/2 ex gi+1/2, j,k ez − gi−1/2, j,k ez − vx,i jk 1z k 1xi ∗ ∗ gi,∗ j,k+1/2 e y − gi,∗ j,k−1/2 e y ¢ª n+1/2 ¡ gi, j+1/2,k ez − gi, j−1/2,k ez − − v y,i jk 1y j 1z k n+1/2 n+1/2 ρi jk h i jk ex n+1/2 n+1/2 ρi jk h i jk e y n+1/2 n+1/2 ρi jk h i jk ez n+1/2
ρi jk
n+1/2
K i jk
(42) In 1D we use the 3D discretization above but retain only terms in ∂/∂ x and ∂ 2 /∂ x 2 and omit derivatives in all transverse directions. n+1/2 In the above expression, time-centered terms (e.g., ρi jk ) are estimated with n+1/2
qi jk
≈
¢ 1¡ n q + q˜ in+1 jk , 2 i jk
(43)
except for the g ∗ s appearing in the (v × G) terms. These are obtained at the half time step and cell edges as components of the Riemann problem solutions.
4. NUMERICAL METHOD: 2D AND 3D
To extend the 1D method described above to multiple spatial dimensions, we use a spatially unsplit fully corner-coupled second-order-accurate scheme after [6, 22]. In 2D, this predictor–corrector approach begins by estimating the 1D x and y fluxes at each cell edge, using the higher order 1D approach described in the previous section. These predictor fluxes, F˜ x and F˜ y , are given schematically as solutions to the Riemann problem R as ¡ ¡ n+1/2 ¢¢ x n+1/2 F˜ i+1/2, j = Fx R qx L ,i+1/2, j , qx R,i+1/2, j ¡ ¡ n+1/2 ¢¢ y n+1/2 F˜ i, j+1/2 = Fy R q y L ,i, j+1/2 , q y R,i, j+1/2 .
(44a) (44b)
142
MILLER AND COLELLA
The predictor fluxes are used to pose a corrector problem, wherein the edge values are augmented by transverse predictor fluxes. Schematically, ¡ ¡ n+1/2 ¢¢ n+1/2 Fx,i+1/2, j = Fx R q 0x L ,i+1/2, j , q 0x R,i+1/2, j ¡ ¡ n+1/2 ¢¢ n+1/2 Fy,i, j+1/2 = Fy R q 0y L ,i, j+1/2 , q 0y R,i, j+1/2
(45a) (45b)
with, for example, ¢ 1t ¡ ˜ y y F i, j+1/2 − F˜ i, j−1/2 21y j X ¡ y ¡ n ¢T ¢T 1t y + 0(γ , δ). v × e y × g˜ i, j+1/2 − g˜ i, j−1/2 γδ 21y j γ δ i, j
q 0 x L ,i+1/2, j = qx L ,i+1/2, j − n+1/2
n+1/2
(46)
The vector 0 is introduced to align the elements of the matrix (v × ∇ × g T )T with the appropriate elements of the vector q: 0(r, s) = (0, 0, 0, 0, 0, δγ 1 δδ1 , δγ 2 δδ1 , δγ 3 δδ1 , δγ 1 δδ2 , δγ 2 δδ2 , δγ 3 δδ2 , δγ 1 δδ3 , δγ 2 δδ3 , δγ 3 δδ3 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)T .
(47)
In 2D there are therefore four Riemann problems solved per cell: two in the predictor and two in the corrector steps. In setting up the corrector step, the components ρ, v, E, g, F p , and κ of the vectors q 0 are updated as indicated in Eq. (46). Our 1D Riemann solver also requires time-centered edge values of the stresses, (σ eα ) L/R in direction eα , and these components of q 0 are calculated by updating the (σ eα ) L/R components of q L/R with the change in stress accompanying the 0 − q L/R in E, g, F p , and κ using cell-centered thermodynamic derivatives. changes q L/R For example, (σ 0 ex ) L ,i+1/2, j = (σ ex ) L ,i+1/2, j +
µ
¯ ¶n ¡ 0 ¢ ∂σ ex ¯¯ E L ,i+1/2, j − E L ,i+1/2, j ¯ ∂E g,F p ,κ i j
¯ ¶n ¡ 0 ¢ ∂σ ex ¯¯ κ L ,i+1/2, j − κ L ,i+1,2, j + ¯ ∂κ g,F p ,E i j ¯ ¶n X µ ∂σ ex ¯ ¡ ¢ ¯ (gγ δ )0L ,i+1/2, j − (gγ δ ) L ,i+1/2, j + ¯ ∂gγ δ g6=gγ δ ,F p ,κ,E i j γδ µ ¶n X ∂σ ex ¯¯ ¡ p¢ ¡¡ p ¢0 ¢ ¯ + Fγ δ L ,i+1/2, j − Fγ δ L ,i+1/2, j . (48) p ¯ p ∂F γ δ g,F p 6=Fγ δ ,κ,E i j γδ µ
By employing this approximation we require only one equation-of-state evaluation per time step per cell for problems involving only elasticity. In problems that also include plasticity, additional equation-of-state evaluations are required for the computation of plastic source terms. In 3D there are two corrector steps: first, ¡ ¡ n+1/2 ¢¢ x n+1/2 F˜ i+1/2, j,k = Fx R qx L ,i+1/2, j,k , qx R,i+1/2, j,k
(49a)
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
¡ ¡ n+1/2 ¢¢ y n+1/2 F˜ i, j+1/2,k = Fy R q y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ n+1/2 ¢¢ z n+1/2 F˜ i, j,k+1/2 = Fz R qz L ,i, j,k+1/2 , qz R,i, j,k+1/2 ,
143
(49b) (49c)
and then ¡ ¡ ¢¢ n+1/2 x|y 0(y) n+1/2 F˜ i+1/2, j,k = Fx R q 0(y) x L ,i+1/2, j,k , q x R,i+1/2, j,k ¡ ¡ ¢¢ n+1/2 x|z 0(z) n+1/2 F˜ i+1/2, j,k = Fx R q x0(z) L ,i+1/2, j,k , q x R,i+1/2, j,k ¡ ¡ ¢¢ n+1/2 y|x 0(x) n+1/2 F˜ i, j+1/2,k = Fy R q 0(x) y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ ¢¢ n+1/2 y|z 0(z) n+1/2 F˜ i, j+1/2,k = Fy R q 0(z) y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ ¢¢ n+1/2 z|x 0(x) n+1/2 F˜ i, j,k+1/2 = Fz R q 0(x) z L ,i, j,k+1/2 , q z R,i, j,k+1/2 ¡ ¡ ¢¢ n+1/2 z|y 0(y) n+1/2 F˜ i, j,k+1/2 = Fz R q 0(y) z L ,i, j,k+1/2 , q z R,i, j,k+1/2 ,
(50a) (50b) (50c) (50d) (50e) (50f)
with, e.g., ¢ 1t ¡ ˜ y y F i, j+1/2,k − F˜ i, j−1/2,k 31y j X ¡ ¢T ¢T ¡ y 1t y + vi,n j,k × e y × g˜ i, j+1/2,k − g˜ i, j−1/2,k 0(γ , δ). (51) γδ 31y j γ δ
n+1/2
q 0(y) x L ,i+1/2, j,k = q x L ,i+1/2, j,k − n+1/2
The final fluxes, which enter the conservative differencing step of the integration, are then computed as ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 (52a) Fx,i+1/2, j,k = Fx R qx L ,i+1/2, j,k , qx R,i+1/2, j,k ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 Fy,i, j+1/2,k = Fy R q y L ,i, j+1/2,k , q y R,i, j+1/2,k (52b) ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 Fz,i, j,k+1/2 = Fz R qz L ,i, j,k+1/2 , qz R,i, j,k+1/2 (52c) with, e.g., 00 n+1/2
¢ ¢ 1t ¡ ˜ y|z 1t ¡ ˜ z|y y|z z|y F i, j+1/2,k − F˜ i, j−1/2,k − F i, j,k+1/2 − F˜ i, j,k−1/2 21y j 21z k ¡ y|z ¢T ¢T 1t X¡ n y|z + 0(γ , δ) vi, j,k × e y × g˜ i, j+1/2,k − g˜ i, j−1/2,k γδ 21y j γ δ n+1/2
qx L ,i+1/2, j,k = qx L ,i+1/2, j,k −
+
¡ z|y ¢T ¢T 1t X¡ n z|y × ez × g˜ i, j,k+1/2 − g˜ i, j,k−1/2 0(γ , δ). v γδ 21z k γ δ i, j,k
(53)
There are a total of 12 Riemann solves per unit cell in 3D: 9 in the predictor steps and 3 in the corrector step. The σ components of the vectors q 0 and q 00 are computed as in the 2D case. 5. PLASTIC SOURCE TERMS
We present here an associated plasticity evolution equation for the rate of change of the plastic deformation tensor F p with time. The more common approach (e.g., [19, 23])
144
MILLER AND COLELLA
is to consider evolution equations for the plastic strain η p = 12 (F pT F p − I ), the plastic Green tensor C p = F pT F p , or the plastic Finger tensor b p = F p F pT . We choose instead to evolve the full nine-component plastic deformation tensor F p . This choice is necessary to be capable of modeling arbitrary crystal systems (see, e.g., [24]). For example, the elastic response of the lowest symmetry crystal system (triclinic) depends upon all six components of the elastic Green tensor. If one were to specify the total inverse deformation g, and either η p , C p , or b p , then all six components of C e could not be determined. Although our examples will make use of isotropic equation of state models (whose elastic invariants may be determined using g and C p ), our goal is to construct a framework of more general applicability. To motivate our choice of evolution equations for F p we begin by postulating the existence of a hyperelastic equation of state, ˆ F p , κ, S), E = E(g,
(54)
where S is the specific entropy. The material derivative of E is ∂E ∂E ˙ 1 ∂E ˙ p ∂E E˙ = κ˙ + S = − σβγ Fγ α g˙ αβ + 8, g˙ αβ + p Fαβ + ∂gαβ ∂κ ∂S ρ ∂Fαβ
(55)
where the second equality equates energy change with the sum of work and heat. Solving for entropy production (dissipation) we have 0 ≤ S˙ = − =
1 T
µ
¶ 1 ∂E ˙ p 1 ∂E σβγ Fγ α ∂E ∂S − + κ˙ + 8 g˙ αβ − p F ρ ∂gαβ T ∂Fαβ αβ T ∂κ ∂E
1 (9plast + 9therm ). ρT
(56)
Here ∂E/∂S = T is the temperature, and we have introduced the specific power of thermal dissipation, 9therm = ρ8,
(57)
and the specific power of plastic dissipation,
9plast = −ρ
∂E ˙ p ∂E − ρ κ˙ p F ∂κ ∂Fαβ αβ
e ˙ Fαβ − ϑ κ˙ = gβγ σγ δ Fδα p
p ˙p Fαβ − ϑ κ˙ = gβγ σγ δ Fδν gνα
= 6 : L p − ϑ κ˙
(58)
p ˙ ˙ Fαβ = (F p )−1 with ϑ = ρ∂E/∂κ being the work-hardening modulus. L νβ = gνα να Fαβ is the p
p
p
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
145
plastic distortion rate [14], and 6βν = gβγ σγ δ Fδν is the thermodynamic force conjugate to L p . The dependence of S˙ on g˙ αβ vanishes because σαβ = −ρ
∂E gγ α . ∂gγβ
(59)
In evaluating (58) we have assumed that E depends on g and F p only through the elastic deformation F e = F g p = (F p g)−1 , e.g., ˜ e , S, κ), E = E(F
(60)
¯ ∂E ¯¯ 1 e = − gβγ σγ δ Fδα . p ρ ∂Fαβ ¯S,κ
(61)
whence
Thermodynamics requires that the internal energy depends upon the volume, and we assume by (60) that this energy dependence is carried by the tensor F e , e.g., V = V0 det F e . For this to be true, it is necessary that det F p = 1 at all times (i.e., V = V0 det F = V0 det F e det F p ; V = V0 det F e iff det F p = 1.) Therefore, (60) assumes that plastic flow is volume-preserving. We postulate a plastic yield surface f = 0, which we represent for illustrative purposes with a Mises–Huber constitutive model written in terms of the Cauchy stress σ , a constant yield stress parameter σY , and the work-hardening modulus ϑ: r f (σ, ϑ) = kdev σ k −
2 (σY + ϑ). 3
(62)
Here, dev σ = σ − 13 (tr σ )I is the stress deviator, and kAk is the Schur norm of A, kAk2 = Aαβ Aαβ = tr(A T A). The flow model we adopt is derived from (62) by the postulate of maximum plastic dissipation [10, 13]. The plastic dissipation (58) is considered as a function of the variables 6 ˙ 9plast = 9plast (6, ϑ; L p , κ). ˙ The plastic dissipation and ϑ, with fixed parameters L p and κ; is then maximized with respect to 6 and ϑ, subject to the constraint that f = 0 during plastic flow. The resulting flow laws are dev(σ ) F˙ p = ζ F p g F kdev(σ )k r 2 κ˙ = ζ 3
(63) (64)
with ζ a parameter chosen to satisfy the Kuhn–Tucker complementarity conditions and the “consistency condition” [23] f =0
(65)
ζ ≥0
(66)
ζf = 0
(67)
ζ f˙ = 0 (if f = 0) .
(68)
146
MILLER AND COLELLA
The flow model (63) is consistent with the assumption that plastic flow is volumepreserving, p p (det˙F p ) = (det F p )gνδ F˙ δν µ ¶ (dev σ )βγ p p = ζ (det F p )gνδ Fδα gαβ Fγ ν kdev σ k
(dev σ )αα kdev σ k = 0 because tr(dev σ ) = 0,
= ζ (det F p )
(69)
and is therefore compatible with the assumption made in evaluating 9plast (58). As an example, we use a modified Mooney–Rivlin equation of state, ρ0 E(C e , S) =
λ(S) √ µ(S) µ(S) (ln det C e )2 + tr C e − log det C e 2 2 2 µ ¶ 1 ρ0 ϑ0 κ + e−ϑ1 κ + ρ ϑ1
(70)
where C e is the elastic Green tensor, C e = F eT F e .
(71)
This equation of state gives a work-hardening modulus, ϑ(κ) = ρ
∂E = ϑ0 (1 − e−ϑ1 κ ), ∂κ
(72)
in terms of two parameters: ϑ0 is the ultimate, asymptotic value of the work-hardening modulus, and ϑ1 dictates the rate of approach of the asymptotic limit. The combined elastic–plastic evolution problem is solved with a predictor–corrector strategy. The inverse total deformation g is advanced in accordance with the equations of motion, with the plastic deformation F p being conservatively advected. This step may predict a coordinate in state space that lies outside the convex manifold of permissible states f (σ, ϑ) ≤ 0, in which case a plastic corrector step is used to bring state back to the yield surface. The algorithmic approach is a return mapping algorithm [23], modified to require only one equation-of-state evaluation. Begin the iteration sequence with iteration index m = 0, F p(0) = F p,n+1 κ (0) = κ n+1 σ (0) = σ (g n+1 , F p,n+1 , κ n+1 )
(73)
ϑ (0) = ϑ(κ n+1 ) There is one equation-of-state evaluation at the beginning of the iteration in which the Cauchy stress σ , work-hardening modulus ϑ, and the derivatives ∂σ/∂F p |E,g,κ , ∂σ/ ∂κ|E,g,F p , and ∂ϑ/∂κ, are calculated. Next, evaluate the yield criterion f (m) = f (σ (m) , ϑ (m) ).
(74)
147
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
If m = 0 and f (m) ≤ ², then the state point is interior to the yield surface, and no plastic flow occurs. If f (0) > ² and | f (m) | ≤ ², then F p,n+1 ← F p(m)
(75)
κ n+1 ← κ (m) and stop. Otherwise, calculate 1ζ (m) = ζ (m+1) − ζ (m) using Newton’s method 1ζ (m) = − f (m)
µ
d f (m) dζ
¶−1 (76)
with d f /dζ estimated from d f (m) ≈ dζ
¯ ¶ ·µ ¶µ ˙ p ¶(m) µ ¯ ¶µ ¶(m) ¸ ∂ f (m) ∂σ ¯¯ ∂σ ¯¯ dF d κ˙ + ¯ ¯ p ∂σ ∂F κ,E,g dζ ∂κ F p ,E,g dζ µ ¶(m) µ ¶µ ¶(m) d κ˙ ∂ϑ ∂f . + ∂ϑ ∂κ dζ µ
(77)
Next, calculate revised estimates µ F˜ p = F p(m) +
∂ F˙ p ∂ζ
¶(m)
1ζ (m)
F p(m+1) = (det F˜ p )−1/3 F˜ p µ ¶(m) ∂ κ˙ (m+1) (m) κ =κ + 1ζ (m) ∂ζ ¯ µ ¶ µ ¯ ¡ p(m+1) ¢ ∂σ ¯¯ ∂σ ¯¯ (m+1) (0) p(0) =σ + −F F σ + ¯ p ∂F ∂κ ¯ E,g,κ
ϑ
(m+1)
= ϑ(κ
(m+1)
(78) ¶
¡
E,g,F p
κ (m+1) − κ (0)
¢
),
set m ← m + 1, and retest the stopping criterion. In this procedure we evaluate the equation of state once to determine σ and the thermodynamic derivatives ∂σ/∂F p |g,E,κ and ∂σ/∂κ|g,E,F p . The stress σ (m) , for m > 0, is approximated by first-order Taylor expansion about the initial m = 0 value. The method converges in 1 or 2 iterations, with ² = 10−6 , in each of the test problems involving plasticity described below. The framework described by Eq. (9) calls for rates of plastic deformation h and rates of work hardening K . In the example above, which is rate-independent, we use τ h = 1F p τ K = 1κ,
(79a) (79b)
where τ = 1t/2 in the predictor step of the method (Eqs. (30a) and (30b), and τ = 1t in the corrector (Eq. (41)). A generalization of this approach to rate-dependent plasticity is described in [19].
148
MILLER AND COLELLA
6. DISSIPATION
In certain problems in hydrodynamics it has been found that the higher order Godunov strategy we adapted here will give rise to spurious post-shock oscillations (e.g., [5]). A solution that rectifies this problem is the addition of a small amount of dissipation at strong shocks. This dissipation is added by introducing an additional slope limiter via a “flattening” parameter χ (see Eq. (27)). A variety of flattening strategies have been proposed. Perhaps the simplest variant, employed by Miller and Puckett [15], uses the divergence of the velocity to detect potential shocks, and uses a simple measure of shock strength, the ratio of pressure jump across a cell to the isentropic bulk modulus, |1P|/K S where K S = ∂ P/∂ log ρ|S , to compute a flattening measure. This introduces additional dissipation in regions where the pressure change is large compared to the bulk modulus—where linearization of the equation of state is expected to become error-prone. This strategy may introduce extra dissipation in regions that do not require it, however, as when a shock is spread over a large (>5 or 6) number of grid cells. It is therefore desirable to also include measures of the shock structure to minimize application of this dissipation mechanism. Elaborate strategies for computing χ are described by Colella and Woodward [7]. One of their strategies is to restrict the use of this dissipative mechanism to regions where the detected shock is steep. In our solid mechanics computations we found this strategy to be useful, and in conjunction with a measure of shock strength it provides judicious, adequate additional dissipation. We detect a strong shock by measuring in 1D the divergence of the velocity field, and calculating a normalized jump in stress. We define zi =
k(σ eα )i+1 − (σ eα )i−1 k∞ (det Aαα,i )1/3
(80)
as a measure of shock strength in the neighborhood of cell i in direction eα . The numerator is the maximum of the absolute value of the jump in those stress components that may change in direction eα 1D purely elastic flow, and the denominator is a mean modulus of the acoustic propagation tensor in direction eα . Following Colella and Woodward, we discriminate between steep and broad shocks by the ratio βi =
k(σ eα )i+1 − (σ eα )i−1 k∞ . k(σ eα )i+2 − (σ eα )i−2 k∞
(81)
In the limit βi = 12 , stress is approximately linear across five grid cells, and so a shock discontinuity is not being captured. When βi ≈ 1 the discontinuity is captured in three cells: the shock may be overly steep and postshock oscillations are expected. Accordingly, the minimum value that our flattening parameter χ should have, based upon shock steepness, is µ µ ¶¶ a1 − βi , , (82) χmin i = max 0, min 1, a1 − a0 where a0 and a1 are numerical constants. We use the values a0 = 0.75 and a1 = 0.85 in the computations presented here.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
149
A local shock-strength-sensitive flattening parameter χ˜ i , χmin i ≤ χ˜ i ≤ 1, is thus ¡ ¡ − zi ¢¢ min 1, max zz11 − , χmin i z0
(v · e j )i+1 < (v · e j )i−1
1
otherwise.
( χ˜ i =
(83)
In our example calculations we use the numerical values z 0 = 0.25 and z 1 = 0.75. In 1D we limit the slopes by the minimum over nearest neighbor cells of the local flattening parameter, χi = min(χ˜ i−1 , χ˜ i , χ˜ i+1 ).
(84)
In higher dimensions, we employ the same 1D local flattening parameters—measured separately in each direction. All slopes (∂q/∂ x, ∂q/∂ y, and ∂q/∂z) are limited by the same cell-valued flattening parameter, which is given by the minimum of the directional local measures. In 2D, χi j = min(χ˜ x,i−1, j , χ˜ x,i, j , χ˜ x,i+1, j , χ˜ y,i, j−1 , χ˜ y,i, j , χ˜ y,i, j+1 ),
(85)
and in 3D, χi jk = min(χ˜ x,i−1, j,k , χ˜ x,i, j,k , χ˜ x,i+1, j,k , χ˜ y,i, j−1,k , χ˜ y,i, j,k , χ˜ y,i, j+1,k , χ˜ z,i, j,k−1 , χ˜ z,i, j,k , χ˜ z,i, j,k+1 ).
(86)
7. ACCURACY
The term (v × ∇ × g T )T was introduced to the evolution equations of the inverse deformation gradient g to make the system of equations stable and well-posed when the gauge constraint ∇ × g T = 0 fails to be satisfied. Although the partial differential equations show that when satisfied initially, it will be satisfied for all times, numerical errors cause the constraint to be violated to some degree. We propose a modification of (9) to control inaccuracy that may arise from violation of the gauge constraint. The conservation law (11) indicates that G will be created by numerical errors as dipoles. Thus, a numerical strategy that will control this truncation error is to diffuse G, ∂G + ∇ · (vG − Gv) = D(∇ 2 G) ∂t
(87)
∂geα ∂ (gv) = (v × G)T eα − D(∇ × G)T eα . + ∂t ∂ xα
(88)
or, equivalently,
g is also related to the density via ρ = ρ0 det g,
(89)
150
MILLER AND COLELLA
where ρ0 is the mass density in the reference state F = g = I . Multiplying the g equations by ρ0 det(g)g −T , and summing over the nine components of g gives a conservation law for ρˆ ≡ ρ0 det(g): ∂ ρˆ + ∇ · (ρv) ˆ = 0. ∂t
(90)
Thus the continuity equation is embodied in the g equations as well. However, because of discretization errors the equivalence of ρˆ and the mass density ρ cannot be assured. To make the method strictly conservative, we keep ρ as a redundant variable, and we invoke a relaxation mechanism on g to enforce the condition ρˆ = ρ. This relaxation alone (not including the diffusion modification) is accomplished by writing ∂ ∂geα (gv) = (v × G)T eα + η + ∂t ∂ xα
µ
¶ ρ − 1 geα . ρˆ
(91)
The “continuity” equation for ρˆ is then Dgαβ D ρˆ ˆ = ρF ˆ βα = −ρ∇ ˆ · v + ρF ˆ βα (v × G)βα + 3η (ρ − ρ) Dt Dt ∂ ρˆ when G = 0. ˆ + ∇ · (ρv) ˆ = 3η (ρ − ρ) ∂t
(92) (93)
Including the diffusion and relaxation terms, the system of equations we will solve is ρvα ρ ρvv − σ e α α ρv ρ Evα − vβ σβα ρE gvδxα gex gvδ yα ∂ ge y ∂ gez + ∂ xα ∂t gvδ zα ρF p ex ρF p ex vα ρF p e y ρF p e v y α p ρF ez ρF p ez vα ρκ ρκvα 0 ρf ρ(8 + v · f ) T −D(∇ × G)T ex e (v × G) 1 (v × G)T e −D(∇ × G)T e 2 y = + (v × G)T e3 T −D(∇ × G) ez ρhe x ρhe y ρhez ρK
0 0 0 +η +η +η 0 0 0 0
¡ ¡ ¡
ρ ρ0 det g ρ ρ0 det g ρ ρ0 det g
¢ − 1 gex ¢ − 1 ge y . ¢ − 1 gez
(94)
151
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
Our discretization of the diffusion and relaxation terms takes the form ¡ ρ˜ n+1 ¢ n+1 − 1 g˜ ex −D(∇ × G)T,n ex + η ρ0 det n+1 g˜ n+1 n+1 ˜ x ge gex ¡ ¢ ρ˜ n+1 T,n n+1 ˜ −D(∇ × G) e + η − 1 g e ge y = ge y y ˜ y + 1t ρ0 det g˜ n+1 ¡ ρ˜ n+1 ¢ n+1 ˜ z i jk ge gez i jk −D(∇ × G)T,n ez + η ρ0 det − 1 g˜ ez n+1 g˜
,
(95)
i jk
where g˜ n+1 denotes g after flux differencing and evaluation of source terms in (10) (cf. (40)). The second derivatives of g appearing in the G diffusion term, £ 2 £ ∂ 2 gx x £ ∂ 2 gx y 2 2 ∂2g ∂ gx z − ∂∂zg2x x − ∂ x 2x y − ∂∂ yg2x z ∂ x∂z ∂ y∂ x ∂z∂ y ¤ ¤ 2 ∂ 2 gx y ¤ ∂ 2 gx z ∂ 2 gx x ∂ 2 gx z + ∂ gx y − ∂ 2 g2x x + ∂ y∂z − ∂z 2 + ∂z∂ x − ∂ x 2 ∂ x∂ y ∂y £ £ £ 2 ∂ 2 g yz ∂ 2 g yx ∂ 2 g yy ∂ 2 g yy ∂ 2 g yz ∂ x∂z − ∂∂zg2yx − − ∂ x∂ y ∂x2 ∂ y∂z ∂ y2 T (∇ × G) = , (96) ∂ 2 g yz ∂ 2 g yy ¤ ∂ 2 g yx ∂ 2 g yz ¤ + ∂ 2 g yy − ∂ 2 g yx ¤ − − + + 2 2 2 ∂ x∂ y ∂y ∂ y∂z ∂z ∂ x∂z ∂x £ ∂ 2 g £ £ 2 2 ∂ gzy ∂ gzy ∂ 2 gzx ∂ 2 gzx ∂ 2 gzz zz − − − ∂ x∂z ∂z 2 ∂ x∂ y ∂x2 ∂ y∂z ∂ y2 ¤ ¤ ¤ 2 2 2 2 2 2 ∂ g ∂ g gzz gzx + ∂ x∂zyy − ∂∂ yg2zx − ∂z 2zy − ∂∂ xg2zz + ∂∂ y∂z + ∂∂ x∂z are computed using time-n cell-centered values of g, with a standard three-point stencil for homogeneous second derivatives, e.g., ¢ ¢¶ ¡ µ ¡ n µ 2 ¶ n 2 gi+1, jk − ginjk 2 ginjk − gi−1, 1 ∂ g jk , (97) = − ∂ x 2 i jk 1xi 1xi+1 + 1xi 1xi + 1xi−1 and heterogeneous derivatives are computed with a four-point stencil, e.g., µ 2 ¶ n n n n gi+1, ∂ g j+1,k − gi+1, j−1,k − gi−1, j+1,k + gi−1, j−1,k . =4 ∂ x∂ y i jk (1xi−1 + 21xi + 1xi+1 )(1y j−1 + 21y j + 1y j+1 )
(98)
A von Neumann stability analysis of the diffusion update in (94), considered independently of other source terms or the basic solid mechanics equations, gives a bound on the diffusion coefficient: 2 h in 1D, 21t D≤ (99) 2 h in 2D or 3D. 41t This suggests an approximate overall Courant–Friedrichs–Lewy stability criterion of ( 1t (|v| + cmax ) CFL = CFL < 1,
1x
+
2D1t (1x)2
in 1D,
1t (|v| + cmax ) 1x
+
4D1t (1x)2
in 2D or 3D,
(100) (101)
where here it is assumed that 1x = 1y = 1z, a constant. The more rigorous CFL condition 4D1t CFL = max( 1t (|v|1x+ cmax ) , (1x) 2 ) (in 2D or 3D) would hold if the mechanics equations and G diffusion steps were performed sequentially.
152
MILLER AND COLELLA
The optimal damping conditions for the diffusion of G are obtained by choosing the empirical diffusion constant D to satisfy D=
(1x)2 , 4d1t
(102)
where d = 1, 2, 3 is the dimensionality of the problem. Similarly, optimal relaxation is obtained by choosing the empirical relaxation parameter η to satisfy η=
1 61t
(103)
for all dimensions. According to the approximate CFL condition (101), the optimal value of D will contribute 1/2 to the CFL value in 1D and 2D, and 1/3 in 3D, limiting the overall step size 1t by factors of 1/2 and 2/3 (respectively) relative to the D = 0 value. Thus, some of the examples presented below use smaller values of D than indicated by (102). In some cases, however, we find that values of CFL > 1 provide stable solutions [consistent with (101) being only an approximation]. Assumptions underlying our plastic yield model require that det F p be constant. The differential equations describing our plastic flow model F˙ p preserves det F p , but again numerical errors will lead to some violation of this constraint. To remedy this problem we renormalize the plastic deformation tensor at the end of each time step, F p ← (det F p )−1/3 F p .
(104)
8. EXAMPLES
8.1. Convergence: Elasticity To demonstrate the convergence properties of the algorithm we model in 1D the smooth flow resulting from an initial Gaussian-shaped disturbance. For these computations we use a hyperthermoelastic model of the Mooney–Rivlin variety, ρ0 E(C e , S) =
√ λ(S) µ(S) µ(S) (log det C e )2 + tr C e − log det C e 2 2 2 µ ¶ 1 −ϑ1 κ ρ0 ϑ0 κ+ e , + ρ ϑ1
(105)
where S is the entropy. Entropy dependence is introduced by supposing λ(S) = λ0 + λS f (S)
(106a)
µ(S) = µ0 + µS f (S)
(106b)
where f (S) is an unspecified function of the entropy. From this equation we evaluate σ (E, g, F p , κ), and other derivatives including the acoustic propagation tensors, by first solving this equation of state for f (S) and then differentiating E with respect to the elements of C e while holding f (S) constant. We use values ρ0 = 1, µ0 = λ0 = 0.6, and µS = λS =
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
153
p
0.01, with initial values g0 = 1.1I , F0 = I , κ = 0, and v = 0. The initial disturbance is generated by distributing internal energy from E 0 to E 00 , ¡ ¢ p E 0 = E g0 , F0 , f (S) = 0 00
E = 10E0 .
(107a) (107b)
These limiting values are used to construct a Gaussian initial profile, via Ei = E 0 ωi + (1 − ωi )E 00
(108)
¸ · 1 r2 ωi = √ exp − i 2 , 2a a 2π
(109)
with
where a 2 = 100 is the variance of the distribution, and where ri is the coordinate of the center of cell i, in the domain [0, 40]. Boundary conditions are reflecting at r = 0 and r = 40. We pick the time step 1t to satisfy the Courant–Friedrichs–Lewy constraint (101), with CFL = 0.8. This problem was chosen to give a nontrivial shockless flow, with initial conditions that strictly obey G = 0. Plasticity is not incorporated into this test problem, and no flattening is required. Figure 1 shows the initial and final conditions of this test problem in Cartesian geometry. At this scale, the difference between results at 40, 80, and 160 Cartesian points is not resolvable. A comparison of results using 40, 80, and 160 grid points is used to estimate the L 1 , L 2 , and L ∞ (max) norm rates of convergence using the volume-weighted variables (Table I). In Cartesian geometry the method exhibits approximately third-order convergence: as the number of grid cells is doubled, the error diminishes by a factor of 23 . Slightly lower rates of convergence are seen in cylindrical and spherical geometries, but in all cases the order exceeds 2.
FIG. 1. Initial conditions and computed results using 160 cells in Cartesian geometry for (a) density, ρ; and (b) stress, σx x .
154
MILLER AND COLELLA
TABLE I Convergence Test: Pure Elasticity Geometry
Field
L1
L2
L∞
Cartesian
ρ vx σx x σ yy , σzz
3.33 3.02 3.30 2.95
3.24 2.97 3.47 2.91
3.06 2.89 3.80 2.69
Cylindrical
ρ vr σrr σθθ σzz
2.51 2.84 2.82 3.12 3.22
2.68 2.70 2.78 3.24 3.17
2.80 2.53 2.68 3.30 2.97
Spherical
ρ vr σrr σθθ , σφφ
2.42 2.77 2.85 3.31
2.53 2.64 2.78 3.34
2.66 2.55 2.79 3.37
8.2. Convergence: Plasticity To assess the rate of convergence in a plasticity-dominated flow we pose a model problem similar to the purely elastic problem presented in (Fig. 2). A Gaussian distribution with width
FIG. 2. Initial conditions and computed results using 160 cells in Cartesian geometry for (a) density, ρ; (b) stress, σx x ; (c) plastic deformation Fxpx ; and (d) work hardening parameter κ.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
155
TABLE II Convergence Test: Elastic–Plastic Flow Geometry
Field
Cartesian
ρ vx σx x σ yy σzz Fxpx F yyp Fzzp κ
L1
L2
L∞
2.14 2.07 1.99 1.62 1.89 2.01 2.18 2.36 2.05
2.03 2.08 1.93 1.55 1.76 1.80 2.08 2.32 1.98
1.84 2.36 2.12 1.21 1.50 1.60 2.10 2.71 1.91
5 is used to vary g yy and gzz as functions of coordinate x according to gx x,i = 1.1 g yy,i = (1 + 9ωi )1.1
(110)
gzz,i = 1.1/(1 + 9ωi ) with homogeneous initial density, internal energy, and zero velocity. We use the equation of state (105) with yield model (62) and flow rates (63, 64). The equation-of-state parameters are as used in the purely elastic convergence test, and the plastic constitutive parameters are σY = 0.1, ϑ0 = 0.1, and ϑ1 = 10.0. The flow field in this problem is C 0 , which lowers the overall order of convergence. Density converges at greater than second order (Table II), but the tangential stress components converge only at first order. 8.3. Blake’s Problem Blake [3] presented an analytical solution to the problem of an unbounded solid medium characterized by an isotropic linear elastic equation of state, ρ0 E =
1 1 λ[tr(C e − I )]2 + µ tr(C eT C e − 2C e − I ), 8 4
(111)
loaded by a prescribed pressure boundary condition on the interior of a spherical cavity of initial radius a. We present a numerical solution to this problem in 1D spherical coordinates (see Appendix), with slight modification of the code to accommodate the moving boundary with prescribed flux (Neumann) boundary conditions. This problem is selected to verify the behavior of the elastic algorithm in the weak shock limit. The cavity wall represents a material interface across which the mass flux will be zero. Accordingly, the flux at this boundary is given by F(U B ) − sU B , where U B is the vector of conserved quantities at the boundary, s is the velocity of the boundary, and F(U B ) is the radial flux vector evaluated at the boundary. Blake’s solution provides u(r, t), the displacement of a mass element in the radial direction. In spherical coordinates, this gives
156
MILLER AND COLELLA
rise to an inverse deformation tensor (1 + ∂u/∂r )−1 0 g(r, t) = 0 (1 + u/r )−1 0
0
0
.
0
(112)
−1
(1 + u/r )
The velocity of the material interface is s = ∂u/∂t|r =a . Cell 1, whose left boundary is r = a at t = 0, and whose right boundary is fixed at a + 1r , has a volume which varies with time. Applying Gauss’s divergence theorem to this cell gives ¡ n+1/2 n+1/2 V1n+1 U1n+1 = V n U1n + 1t A1/2 F1/2 − A3/2 F3/2 ¡ n+1/2 n+1/2 ¢¢ n+1/2 + A¯ 1 H1/2 − H3/2 , (113) + 1t V¯ 1 G 1 where F denotes the radial flux component that enters as (1/r 2 )∂(r 2 F)/∂r , H denotes the radial flux component that enters as ∂ H/∂r (see Appendix), A¯ 1 is the average area (r 2 ) n+1/2 over r in [a − u(a, t), a + 1r ], V¯ 1 is the time-averaged cell volume, and G 1 is the cell-centered vector of (geometric) source terms, which we time-center with a predictor– corrector strategy. In general (see Wilkins’s problem below), an algebraic solution of this discretization is unstable. In the particular case of our discretization of Blake’s problem, however, |u(a, t)| ¿ 1r and so V1 does not vary appreciably with time and in particular is of order a 2 1r . Our solution to Blake’s problem therefore uses (113) as written. It is also necessary to modify the algorithm to account for the absence of cell values at i − 1 and i − 2. The gradient ∂q/∂r at i = 1 is obtained by first order forward finite difference with a van Leer limiter. The flattening parameter χ operates on a stencil that requires cell values at 0 and −1. However, for this weak problem additional flattening is never required, so the algorithm is modified by omission of the flattening computation (χ = 1). Following Trangenstein and Colella [25] we use parameters a = 0.1 m, ρ0 = 3000 kg/m3 , λ = 2.36 × 1010 Pa, and µ = 2.78 × 1010 Pa. The pressure inside the spherical cavity is 106 Pa, and the solution is plotted at time 1.6 × 104 s. We compare in Figs. 3–6 our computed results for radial stress, σrr = (λ + 2µ)(∂u/∂r ) + 2λ(u/r ),
(114)
σθ θ = σφφ = λ(∂u/∂r ) + 2(λ + µ)(u/r ),
(115)
¶ µ 1 2 P = − σii = − λ + µ [(∂u/∂r ) + 2(u/r )] , 3 3
(116)
hoop stress
pressure
and radial velocity vr against Blake’s analytical results. These results verify the method in the case of weak (linear) waves. The leading shock is captured in approximately five grid cells. A single stress undershoot precedes the shock, and a corresponding overshoot follows it, but the wave speed and amplitude are correctly modeled.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
157
FIG. 3. σrr calculated from Blake’s analytical solution at time 1.6 × 104 , compared with values computed using 200 grid cells.
FIG. 4. Hoop stress, σθθ = σφφ .
FIG. 5. Pressure, (σrr + σθ θ + σφφ )/3.
158
MILLER AND COLELLA
FIG. 6. Material velocity vr for Blake’s problem.
8.4. Wilkins’s Problem Wilkins’s flying plate problem [27] involves a 5-mm-thick aluminum plate impacting an initially stationary aluminum halfspace. The rear (left) surface of the flying plate is a free surface (vacuum). Initially, left- and right-traveling shocks propagate outward from the point of contact of the plate with the halfspace. When the left-traveling shock reaches the free surface, a right-traveling rarefaction is created, which ultimately overtakes the right-traveling shock. This problem incorporates plasticity. To model this problem, we modify our 1D algorithm to allow for the moving free-surface boundary. This is an example of volume-of-fluid front reconstruction applied to multi-fluid modeling, and details will be described in a future correspondence. Briefly, we modify the approach adopted for Blake’s problem using the flux redistribution ideas of Chern and Colella [4]. Application of this approach to stationary incompressible boundaries is described in [16], and to reaction front tracking in [1, 18]. Our implementation is similar, but the free-surface boundary moves at a velocity determined by the solid-vacuum Riemann problem. This problem is solved as described above for the solid–solid case but uses only the 3 × 3 stress component of the eigenvectors. This interface velocity, and the surrounding material velocities, are used with a volume-pushing algorithm (after [2]) to update the fractional occupancy of the interface cells. We construct a hyperelastic model of aluminum in close correspondence to Wilkins’s (rate model) description, with ¶ µ µZ µ ¶2/3 ¶ µ0 P(ρ 0 ) 0 ρ0 e + tr C , (117) dρ − 3 E(g, F p ) = 02 ρ 2ρ0 ρ where P(ρ) is the hydrostatic pressure (in GPa) P(ρ) = 72(ρ/ρ0 − 1) + 172(ρ/ρ0 − 1)2 + 40(ρ/ρ0 − 1)3 ,
(118)
with ρ0 = 2.7 kg/m3 . The shear modulus is µ0 = 24.8 GPa. The problem is perfectly plastic (no work hardening), and uses the von Mises yield surface function r 2 (119) σY f (σ ) = kdev σ k − 3 with constant flow stress σY = 0.2976 GPa.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
159
FIG. 7. Longitudinal stress σx x for Wilkins’s problem with impact velocity 0.8 km/s. Time in µs.
Computations with impact velocities of 0.8 km/s (Figs. 7 and 8) and 2.0 km/s (Figs. 9 and 10) were obtained with CFL = 0.80 and 500 Cartesian grid points. At 0.8 km/s, a plastic shock trails a leading elastic shock precursor. When the left-facing shocks reach the free surface, right-traveling elastic and trailing plastic rarefaction waves begin to overtake the initial right-facing shocks. The shock stress at 2.0 km/s is above the elastic limit, so only plastic shocks are formed. On rarefaction from the left free surface, a leading rightfacing elastic rarefaction is formed, followed by the plastic wave. These results are in good quantitative agreement with those of Wilkins. 8.5. A Test in 2D This test problem compares a 1D cylindrical coordinate computation against a 2D Cartesian result, for a problem with cylindrical symmetry. We use the modified Mooney– Rivlin model presented in Eq. (70) with initial conditions ρ0 = 1, g = 1.1I , F p = I , and κ = 0. The plasticity parameters are σY = 0.1, ϑ0 = 0.1, and ϑ1 = 10. All boundary conditions are reflecting. The material is initially at rest, except for a cylindrical shell r ∈ [5, 15], which moves toward the axis with a velocity of −1. This generates a diverging rarefaction, and a convergent shock, which reflects off the axis of symmetry.
160
MILLER AND COLELLA
FIG. 8. Mass density ρ for Wilkins’s problem with impact velocity 0.8 km/s. Time in µs.
In Figs. 11–15 we compare results from a 1D cylindrical calculation (500 cells, CFL = 0.8), and an equivalent 2D Cartesian calculation using 250 × 250 cells, also at CFL = 0.8. The 2D results are presented as 1D scatter plots in order to demonstrate the accurate preservation of cylindrical symmetry obtained with the spatially unsplit 2D method. The high-resolution 1D results and lower resolution 2D results are in good agreement, although there is some discrepancy in κ and σrr near the axis. Using this same 2D test we demonstrate the errors associated with the gauge constraints ρ − ρ0 det(g) = 0
(120)
G = ∇ × g T = 0.
(121)
and
These conditions are enforced in the computation by way of a relaxation term to satisfy (120) and a diffusion-like term to satisfy (121). In Fig. 16 we plot the left-hand side of
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
161
FIG. 9. Longitudinal stress σx x for Wilkins’s problem with impact velocity 2.0 km/s. Time in µs.
(120) comparing results from the computation presented in Figs. 11–15, and results from a similar computation in which, however, neither a relaxation nor a diffusion correction was applied. In Fig. 17 we plot the L 2 norm of the tensor ∇ × g T , comparing results from the computation with relaxation and diffusion to results from a computation using neither correction. These figures demonstrate over an order of magnitude reduction in density error is achieved by the relaxation mechanism. Approximately a factor of 2 reduction of ||G||2 is achieved by the diffusion mechanism. 8.6. A Test in 3D This test problem compares a 1D spherical coordinate computation against a 3D Cartesian result, for a problem with spherical symmetry. The equation of state is identical to the 2D test above, and the initial conditions are similar: a spherical shell r ∈ [5, 15] is given an initial velocity of −1. This computation, with 100 × 100 × 100 cells at CFL = 0.8, is underresolved. Nevertheless, there is good agreement between the 3D Cartesian results and the 1D spherical calculation and excellent preservation of spherical symmetry (Figs. 18–20).
162
MILLER AND COLELLA
FIG. 10. Mass density ρ for Wilkins’s problem with impact velocity 2.0 km/s. Time in µs.
FIG. 11. Density at time 4.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
FIG. 12. Radial velocity at time 4.
FIG. 13. σrr at time 4.
FIG. 14. F p rr component of plastic deformation tensor.
163
164
MILLER AND COLELLA
FIG. 15. Work-hardening parameter κ.
FIG. 16. Error in density constraint, ρ − ρ0 det(g).
FIG. 17. Error in curl constraint,
p
Gαβ Gαβ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
FIG. 18. Density at time 4.
FIG. 19. Radial velocity at time 4.
FIG. 20. σrr at time 4.
165
166
MILLER AND COLELLA
9. CONCLUSIONS
We have presented a new method for the solution to equations of solid mechanics in one, two, and three spatial dimensions on Eulerian grids. Our method addresses the problem of gauge constraints (∇ × g T = 0) by adopting a nonconservation approach first proposed by [8] for the equations of magnetohydrodynamics. We write the partial differential equations of solid mechanics in such a way that the constraint, if applicable initially, holds true for all time. The constraint is violated by the truncation error of the method, and reinforced with an explicit diffusion term which annihilates the dipolar field of ∇ × g T . Another constraint of the system, a correspondence between density variation and the deformation field (ρ = ρ0 det g) is also satisfied for all times by the PDEs, if satisfied in the initial conditions. Truncation errors of the method are compensated with an explicit relaxation term. The method presented here does not incorporate artificial viscosity, but its solutions are sensitive to six adjustable parameters: D and η control accuracy of the gauge constraints, and a0 , a1 , z 0 , and z 1 in Eqs. (82) and (83) govern the introduction of dissipation near strong shocks to prevent overshoot and ringing by locally reducing the high-order Godunov method to first order. Our strategy for damping modes violating the curl gauge constraint, g T : = g T − λ∇ × ∇ × g T , = g T + λ(∇ 2 g T − ∇(∇ · g T )),
(122)
= g + λ(∇ g − ∇ Q(g )), T
2 T
2
T
(λ = 1tD), uses a single central difference operator acting on cell-centered variables. Here, Q(g T ) = ∇ −2 ∇(∇ · g T ) is the projection onto the curl-free part of g T . Defining P(x) = 1 − Q(x) as the projection onto the divergence-free part of x, and noting P Q = Q P = 0, we have P(g T ) := P(g T ) + λ∇ 2 P(g T );
(123)
thus we are diffusing the divergence-free part of g T without modifying the curl-free part. A similar scheme may be used to modify a vector field B subject to a divergence-free constraint, B := B + λ∇(∇ · B), Q(B) := Q(B) + λ∇ 2 Q(B)
(124)
with a single matrix-valued central difference operator for the projection ∇(∇ · B). This will directly target odd–even and checkerboard short-wavelength modes of ∇ · B by diffusing the curl-free part of B. The application of this extension to magnetohydrodynamics, where B is the magnetic field subject to gauge constraint divB = 0, is currently being investigated (R. Crockett, personal communication). APPENDIX: CYLINDRICAL AND SPHERICAL COORDINATES
The equations of solid mechanics in cylindrical and spherical coordinates (like those of gas dynamics) differ from the Cartesian equations by the existence of both spatial and
167
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
volumetric spatial derivatives, and by the introduction of “geometric source terms.” The coordinate transformation is accomplished by rotating the Cartesian basis vectors into the curved coordinate frame via the rotation matrices
Rcyl
0 0 1
cos θ −sin θ = sin θ cos θ 0 0
sin θ cos φ
cos θ cos φ
cos θ
cos θ sin φ −sin θ
Rsph = sin θ sin φ
−sin φ
(A.1)
cos φ , 0
(A.2)
where we adopt the standard curved coordinate notation x = r cos θ y = r sin θ
(A.3)
z=z in cylindrical coordinates and x = r sin θ cos φ y = r sin θ sin φ
(A.4)
z = r cos θ in spherical coordinates. R is the matrix of inner products of unit vectors in the curved coordinate system, eα0 , and the Cartesian system eβ ; Rαβ = eα0 · eβ . These rotation matrices transform the Cartesian tensors F, g, and σ transform as σcyl = R T σCart R, etc. and transform the velocity vector v as vcyl = R T vCart . In cylindrical coordinates, the system of transformed equations may be written (cf. Eq. 9) as
0 ! à ! ρvθ σrr ρ ρvv − σ0 rθ ρvvθ − σ eθ r ρv − 0 σr z ρE ρ Evθ − v T σ eθ 0 ge 0 ρ Evr − v T σ er 0 r ∂ 1 ∂ 1 ∂ ∂ gv gv ge 0 θ + + + r r ∂r 0 0 ge ∂t ∂r r ∂θ z 0 p p ρF er 0 ρF er vθ 0 p ρF p e v p 0 ρF eθ ρF e v θ θ r r 0 ρF p e ρF p e v p z θ r ρF e v z θ 0 ρF p ez vr ρκ ρκv ρvà r
ρκvr
0
θ
168
MILLER AND COLELLA
0 ¡ ¢ − ρv 2 + σrr − σθ θ (ρvθ v − σ ) r θ rθ 0 0 0 −g v − g v − g v θr r θθ θ θz z ρvz grr vr + gr θ vθ + gr z vz ρvvz − σ ez 0 ρ Ev − v T σ e z z 0 0 p p −Fr θ − Fθr ∂ 1 0 + = − p gv ∂z r ρvθ Frrp − Fθ θ p ρF er vz p −F zθ p ρF eθ vz ρF p e v p p z z Frr − Fθ θ ρκvz p p F + F ρv θ rθ θr p Fzr p −Fθ z Frpz ρv θ 0
0 0 ρf ρ(8 + v · f ) rr − vz ∂g∂zrr − vθ ∂g + r ∂θ
v ∂gr θ + v ∂gr z vθ (gr θ + gθr ) θ ∂r z ∂r r ∂gθθ ∂gθ z ∂gθr ∂gθr vθ (gθθ − grr ) vθ ∂r + vz ∂r − vθ r ∂θ − vz ∂z + r zr vθ ∂g∂rzθ + vz ∂g∂rzz − vθ ∂g − vz ∂g∂zzr + vθ rgzθ r ∂θ v ∂grr + v ∂gr z − v ∂gr θ − v ∂gr θ − vr (gr θ + gθr ) + vz gθ z r r ∂θ z r ∂θ r ∂r z ∂z r v ∂gθr + v ∂gθ z − v ∂gθθ − v ∂gθ θ + vr (grr − gθ θ ) + vz gr z r r ∂θ z r ∂θ r ∂r z ∂z r . + vr gzθ ∂gzr ∂gzz ∂gzθ ∂gzθ vr r ∂θ + vz r ∂θ − vr ∂r − vz z − r ∂grr ∂gr θ ∂gr z ∂gr z vθ gθ z vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ + r vθ gr z ∂gθr ∂gθθ ∂gθ z ∂gθ z vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ − r ∂gzr ∂gzθ ∂gzz ∂gzz vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ ρher ρheθ ρhez
(A.5)
ρK
This equation does not include the det g relaxation term, whose representation is unaffected by the change in variables, nor does it include the G diffusion correction, which will be described separately below. There is some latitude in the partitioning of terms between the LHS and the RHS geometric source vector. This is particularly evident in the stress terms appearing in the momentum equations. The choice of representations described here was chosen in order that some cancellation between σrr and σθ θ occur in the r -momentum source term.
169
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
The linearized equations of solid mechanics [cf. (28)], used in the construction of L and R edge states, also has a geometric source vector. Expressed in terms of the primitive variables q, but omitting stress components which are described later, we have −ρvr vθ2 + σrr /ρ − σθθ /ρ −v v + 2σ /ρ r θ rθ σr z /ρ (vr σθθ − vθ σr θ )/ρ Ã ! 0 0 0 Ã v g +v g +v g ! r θr θ θθ z θz −vr grr − vθ gr θ − vz gr z 0 Ã ! 0 1 . s= (A.6) 0 r 0 p Fr θ + Fθrp p p ρvθ Fθ θ − Frr p Fzθ p Fθ θ − Frrp ρvθ −Frpθ − Fθrp −Fzrp p Fθ z ρvθ −Frpz 0 0
The stress evolution equations, used in the predictor steps of the method, are (plastic source terms omitted) σ er ∂ σ eθ ∂t σ ez
+ vr
Arr = Aθr Azr
σ er ∂ σ eθ ∂r σ ez
vr · ∂ vθ ∂r vz
+ vθ
σ er ∂ σ eθ r ∂θ σ ez
Ar θ + Aθ θ Azθ
+ vz
vr · ∂ vθ r ∂θ vz
2vθ σr θ + vr (Ar θ )r θ − vθ (Ar θ )rr
σ er ∂ σ eθ ∂z σ ez
Ar z
where the tensors A are defined by (14).
vr + A θ z · ∂ vθ ∂z vz Azz
− vθ (σrr − σθθ ) + vr (Ar θ )θ θ − vθ (Ar θ )θr vθ σθ z + vr (Ar θ )r θ − vθ (Ar θ )zr −v (σ − σ ) + v (A ) − v (A ) θ rr θθ r θθ rθ θ θ θ rr 1 , + −2vθ σr θ + vr (Aθ θ )θ θ − vθ (Aθ θ )θr r −vθ σr z + vr (Aθθ )zθ − vθ (Aθθ )zr vθ σθ z + vr (Azθ )r θ − vθ (Azθ )rr −vθ σr z + vr (Azθ )θθ − vθ (Azθ )θr +vr (Azθ )zθ − vθ (Azθ )zr
(A.7)
170
MILLER AND COLELLA
The g relaxation term,−D(∇ ×∇ ×g T )T , transforms in cylindrical coordinates as −D × £ ∂ 2 g rz ∂r ∂z
∂ 2 gr θ r ∂r ∂θ ∂ 2 grr + grrr2 − grθθ2 r 2 ∂θ 2 ∂gr θ θr θθ + 2∂g − ∂g r 2 ∂θ r 2 ∂θ r ∂r
−
∂ 2 grr ∂z 2
+
£ ∂ 2 grr
2 2 rz − ∂∂rg2r θ + r∂∂θg∂z r ∂r ∂θ 2 rr − ∂∂zg2r θ − r∂g2 ∂θ + grr2θ + grθr2 ¤ rθ θr θz − ∂g − ∂g − ∂g r ∂r r ∂r r ∂z
£ ∂ 2 gr θ
2 gr z ∂ 2 grr − r∂2 ∂θ 2 + ∂r ∂z r ∂θ ∂z 2 θz − ∂∂rg2r z + grr2z + 2∂g r 2 ∂θ ¤ rr rz θθ − ∂g − ∂g + ∂g r ∂z r ∂r r ∂z
− ¤ + £ ∂ 2 g £ £ 2 2 2 2 2 2 2 2 ∂ gθr ∂ gθθ ∂ gθr ∂ gθ θ ∂ gθ z ∂ gθ θ ∂ gθ z ∂ gθr θz − + − + − + ∂r ∂z ∂z 2 r ∂r ∂θ r ∂r ∂θ ∂r 2 r ∂θ ∂z r ∂θ ∂z r 2 ∂θ 2 ∂r ∂z 2 2 2 ∂gθθ ∂ gθ θ ∂gθr ∂grr ∂ gθ z 2∂gr z ∂gr θ rr − ∂2 gθr2 − 2∂g + r 2 ∂θ − ∂z 2 − r 2 ∂θ + r ∂r − ∂r 2 − r 2 ∂θ + r ∂z r ∂θ r 2 ∂θ ¤ ¤ ¤ ∂g g g ∂gθ θ grr gθθ ∂gθr ∂gθ z gθ z rz − − + − + + ∂g + + r ∂rr θ + rr2θ + rθr2 2 2 2 r ∂z r ∂r r r r ∂z r ∂r r £ ∂ 2 g £ £ ∂ 2 gzr ∂ 2 gzθ ∂ 2 gzr ∂ 2 gzθ ∂ 2 gzz ∂ 2 gzθ ∂ 2 gzz ∂ 2 gzr zz − + − + − + ∂r ∂z ∂z 2 r ∂r ∂θ r ∂r ∂θ ∂r 2 r ∂θ ∂z r ∂θ ∂z r 2 ∂θ 2 ∂r ∂z ¤ ¤ ¤ 2 2 2 ∂ gzr ∂gzθ ∂ gzθ ∂gzr ∂gzθ gzθ ∂ gzz ∂gzr ∂gzz − ∂z 2 − r 2 ∂θ − r ∂r + r 2 − ∂r 2 + r ∂z − r ∂r − r 2 ∂θ 2 + r 2 ∂θ (A.8)
The transformed system of equations in spherical coordinates may be written as ρ ρvr ρv 0 ρvvr −σ er ρE ρ Evr − v T σ er ger 0 0 ∂ geθ 1 ∂ 2 ∂ gv 0 + r + ∂r 0 ∂t geφ r 2 ∂r 0 ρF p er ρF p er vr 0 ρF p e ρF p e v 0 θ θ r p p 0 ρF eφ ρF eφ vr ρκ ρκvr
ρvθ ρvvθ T ρ Evθ − v σ eθ 0 1 ∂ 0 ∂ 1 sin θ + + r ∂θ 0 r sin θ ∂θ ρF p er vθ ρF p e v θ θ ρF p eφ vθ ρκvθ
ρvφ ρvv − σ e 0 φ φ T −σ eθ ρ Ev − v σ eφ φ 0 0 0 0 ∂ + 1 gv r sin θ ∂φ gv ρF p er vφ 0 ρF p e v 0 θ φ p 0 ρF eφ vφ ρκvφ
171
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
0 −ρ vθ2 + vφ − 2σrr + σθ θ + σφφ − cot θ σr θ ¢ ¡ ρ vr vθ − cot θvφ2 − 3σr θ + cot θ (σφφ − σθ θ )
¡
¢ 2
ρ(vr vφ + vθ vφ cot θ ) − 3σr φ − 2 cot θ σθ φ 0 0 −gθr vr − gθ θ vθ − gθ φ vφ grr vr + gr θ vθ + gr φ vφ 0 −gφr vr − gφθ vθ − gφφ vφ −cot θ(gφr vr + gφθ vθ + gφφ vφ ) (grr vr + gr θ vθ + gr φ vφ ) + cot θ (gθr vr + gθ θ vθ + gθ φ vφ ) 1 =− ¡ p ¡ p p¢ p ¢ r −vθ Fr θ + Fθr − vφ Fr φ + Fφr ¡ p p¢ p p ρ vθ Frr − Fθ θ − vφ Fθ φ − cot θ vφ Fφr ¡ ¢ p p p p −vθ Fφθ + vφ Frr − Fφφ + cot θ vφ Fθr ¡ p p¢ p p vθ Frr − Fθ θ − vφ Fφθ − cot θ vφ Fr φ ¡ ¢ ¡ ¢ p p p p ρ vθ Fr θ + Fθr − cot θ vφ Fφθ + Fθ φ ¡ ¢ p p p p F + v F + cot θ v − F v F θ φ φ φr rθ θθ φφ ¡ p ¢ p p p −vθ Fθ φ + vφ Frr − Fφφ + cot θ vφ Fr θ ¡ ¢ p ρ vθ Frpφ + vφ Fθrp + cot θ vφ Fθpθ − Fφφ ¡ p ¡ p p ¢ p ¢ vφ Fφr + Fr φ + cot θ vφ Fφθ + Fθ φ 0 0 ρf ρ(8 + v · f ) ¡ ∂gr θ ∂grr ¢ ¡ ∂grr ¢ ∂g v (g + g ) + v (g + g ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂rr φ + θ r θ θr r φ r φ φr ¡ ∂gθr ¡ ∂gθθ ∂gθr ¢ ∂gθφ ¢ vθ (gθθ − grr ) + vφ (gθ φ + cot θgφr ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂r + r ¡ ¢ ¡ ¢ ∂gφθ ∂gφr ∂gφr ∂gφφ vθ gφθ + vφ (gφφ − grr − cot θgθr ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂r + r ¡ ∂g ¢ ¡ ¢ −vr (gr θ + gθr ) + vφ (gφθ − gθ φ + cot θgr φ ) ∂gr θ ∂gr θ ∂grr rφ vφ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r ¡ ∂g ¢ ¡ ∂gθ θ ∂gθr ¢ vr (grr − gθ θ ) + vφ (gr φ + cot θ (gθφ + gφθ )) ∂gθθ θφ v φ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r + . ¡ ∂g ¢ ¡ ¢ ∂gφθ ∂gφθ ∂gφr −vr gφθ + vφ (−gr θ + cot θ (gφφ − gθθ )) φφ vφ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r ¡ ∂grr ¢ ¡ ¢ ∂gr φ ∂gr φ −vr (gφr + gr φ ) + vθ (gθ φ − gφθ − cot θgr φ ) ∂gr θ vr r sin θ ∂φ − ∂r − vθ r ∂θ − r sin θ ∂φ + r ¡ ∂gθr ¢ ¡ ¢ ∂gθφ ∂gθφ −vr (gθφ + cot θgφr ) − vθ (gr φ + cot θ (gφθ + gθφ )) ∂gθ θ vr − − − v + θ r ∂θ r sin θ ∂φ ∂r r sin θ ∂φ r ¢ ¡ ¢ ¡ ∂gφr ∂gφφ ∂gφφ ∂gφθ vr (grr − gφφ + cot θgθr ) + vθ (gr θ + cot θ (gθ θ − gφφ )) vr − ∂r − vθ r ∂θ − r sin θ ∂φ + r sin θ ∂φ r ρhe r ρhe θ ρheφ ρK
[A.9]
172
MILLER AND COLELLA
Again, the det g source term, not included above, is unaffected by the transformation of variables. The G diffusion term is described separately below. The geometric source terms in the vector s [cf. (28)] corresponding to the primitive variables q, but omitting the direction-dependent stress terms, are
−ρ(2vr + cot θ vθ ) 2 2 vθ + vφ + (2σrr − σθ θ − σφφ + cot θ σr θ )/ρ
−vr vθ + cot θ v 2 + (3σr θ + cot θ (σθ θ − σφφ ))/ρ φ −vr vφ − cot θ vθ vφ + (3σr φ + 2 cot θ σθ φ )/ρ (vr (σθ θ + σφφ ) + vθ (−σθr + cot θ σφφ ) − vφ (σφr + cot θ σθ φ ))/ρ 0 0 0 vr gθr + vθ gθ θ + vφ gθ φ −vr grr − vθ gr θ − vφ gr φ 0 vr gφr + vθ gφθ + vφ gφφ cot θ (vr gφr + vθ gφθ + vφ gφφ ) 1 . s= −vr grr − vθ gr θ − vφ gr φ − cot θ (vr gθr + vθ gθ θ + vφ gθ φ ) r ¡ p ¢ ¡ ¢ p p p vθ Fr θ + Fθr + vφ Fr φ + Fφr ¡ p ¢ vθ F − F p + vφ F p + cot θ vφ F p ρ θθ θφ φr rr ¡ ¢ p p p p vθ Fφθ + vφ Fφφ − Frr − cot θ vφ Fθr ¡ p ¢ p p p vθ Fθ θ − Frr + vφ Fφθ + cot θ vφ Fr φ ¡ p ¢ ¡ ¢ −vθ F + F p + cot θ vφ F p + F p ρ rθ θr φθ θφ ¡ p ¢ p p p −vθ Fφr − vφ Fr θ + cot θ vφ Fφφ − Fθ θ ¡ p ¢ p p p vθ Fθ φ + vφ Fφφ − Frr − cot θ vφ Fr θ ¡ p ¢ p p p −vθ F − vφ F + cot θ vφ F − F ρ rφ θr φφ θθ ¡ p ¡ p p ¢ p ¢ −vφ Fφr + Fr φ − cot θ vφ Fφθ + Fθ φ 0 (A.10) The stress evolution equations, used in the predictor steps of the method, are (nongeometric source terms omitted)
σ er σ er σ er σ er ∂ ∂ ∂ ∂ σ eθ σ eθ + vr σ eθ + vθ σ eθ + vφ ∂t ∂r r ∂θ r sin θ ∂φ σ eφ σ eφ σ eφ σ eφ vr vr vr Ar θ Ar φ Arr ∂ ∂ ∂ vθ = Aθr · vθ + Aθ θ · vθ + A θ φ · ∂r r ∂θ r sin θ ∂φ Aφr vφ Aφθ vφ Aφφ vφ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
173
[2(vθ σr θ + vφ σr φ ) + vr (Ar θ )r θ − vθ (Ar θ )rr + (vr + vθ cot θ )(Ar φ )r φ − vφ (Ar φ )rr − cot θ vφ (Ar φ )r θ ] [−vθ (σrr − σθ θ ) + vφ (σθ φ + cot θ σr φ ) + vr (Ar θ )θ θ − vθ (Ar θ )θr + (vr + vθ cot θ)(Ar φ )θ φ − vφ (Ar φ )θr − cot θ vφ (Ar φ )θ θ ] σ − v (σ − σ + cot θ σ ) + v (A ) − v (A ) [v θ θφ φ rr φφ rθ r rθ rθ θ r θ φr + (vr + vθ cot θ )(Ar φ )φφ − vφ (Ar φ )rr − cot θ vφ (Ar φ )φθ ] [−vθ (σrr − σθ θ ) + vφ (σθ φ + cot θ σr φ ) + vr (Aθ θ )r θ − vθ (Aθ θ )rr + (vr + vθ cot θ)(Aθ φ )r φ − vφ (Aθ φ )rr − cot θ vφ (Aθ φ )r θ ] σ − v cot θσ ) + v (A ) − v (A ) + (v + v cot θ ) [−2(v 1 θ rθ φ θφ r θθ θθ θ θ θ θr r θ + r × (Aθ φ )θ φ − vφ (Aθ φ )θr − cot θ vφ (Aθ φ )θ θ ] [−vθ σr φ − vφ (σr θ + cot θ (σθ θ − σφφ )) + vr (Aθ θ )φθ − vθ (Aθ θ )φr + (vr + vθ cot θ)(Aθ φ )φφ − vφ (Aθ φ )φr − cot θ vφ (Aθ φ )φθ ] [vθ σθ φ − vφ (σrr − σφφ + cot θ σr θ ) + vr (Aφθ )r θ − vθ (Aφθ )rr + (vr + vθ cot θ )(Aφφ )r φ − vφ (Aφφ )rr − cot θ vφ (Aφφ )r θ ] [−vθ σr φ − vφ (σr θ + cot θ (σθ θ − σφφ )) + vr (Aφθ )θ θ − vθ (Aφθ )θr + (vr + vθ cot θ)(Aφφ )θ φ − vφ (Aφφ )θr − cot θ vφ (Aφφ )θ θ ] [−2vφ (σr φ + cot θσθφ ) + vr (Aφθ )φθ − vθ (Aφθ )φr + (vr + vθ cot θ ) × (Aφφ )φφ − vφ (Aφφ )φr − cot θ vφ (Aφφ )φθ ]
The g relaxation term −D(∇ ×∇ ×g T )T transforms in spherical coordinates as ((∇ × ∇ × g T )T )rr =
((∇ × ∇ × g T )T )r θ =
((∇ × ∇ × g T )T )r φ =
∂ 2 gr φ ∂ 2 gr θ ∂gr θ ∂ 2 grr ∂ 2 grr 2∂gθr + − 2 2 − 2 2+ 2 + 2 2 r sin θ∂r ∂φ r sin θ ∂φ r ∂r ∂θ r ∂θ r ∂θ r ∂θ ∂gθ θ 2∂gφr ∂gr φ cot θ ∂grr ∂gφφ − + 2 + − − r ∂r r ∂r r sin θ ∂φ r 2 sin θ ∂φ r 2 ∂θ 2grr gθ θ cot θ gr θ 2cot θ gθr gφφ cot θ∂gr θ + − 2 + 2 − 2 + + r ∂r r r r r2 r2 (A.12a) ∂gφφ ∂ 2 gr θ ∂ 2 gr θ 2∂gr θ ∂ 2 grr ∂ 2 gr φ − − 2 2 2 − 2 − + 2 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂ φ r ∂θ r ∂r cot θ ∂gr φ 2∂gφθ ∂gθ φ gr θ ∂gθr + 2 + − + 2 − r ∂r r sin θ ∂φ r 2 sin θ ∂φ r 2 sin θ ∂φ r cot θ gφφ cot θ gθ θ − + (A.12b) 2 r r2 ∂ 2 gr φ ∂ 2 gr φ 2∂gθ φ ∂ 2 gr θ ∂ 2 grr ∂gφθ − − + − 2 + 2 r 2 sin θ∂θ∂φ r 2 ∂θ 2 r sin θ ∂r ∂φ ∂r 2 r ∂θ r ∂θ 2∂gr φ cot θ ∂gr θ cot θ ∂gr φ ∂gθ θ gr φ ∂gφr − − 2 − − 2 + 2 − r ∂r r ∂r r sin θ ∂φ r 2 ∂θ r sin θ ∂φ r gr φ cot θ gθ φ cot θ gφθ + 2 2 + + (A.12c) 2 r r2 r sin θ
174
MILLER AND COLELLA
((∇ × ∇ × g T )T )θr =
((∇ × ∇ × g T )T )θ θ =
((∇ × ∇ × g T )T )θ φ =
((∇ × ∇ × g T )T )φr =
((∇ × ∇ × g T )T )φθ =
((∇ × ∇ × g T )T )φφ =
∂ 2 gθr ∂ 2 gθr ∂gθ θ ∂ 2 gθ φ ∂ 2 gθ θ 2∂grr − 2 2 − + 2 + − 2 r sin θ∂r ∂φ r sin θ ∂φ 2 r ∂r ∂θ r 2 ∂θ 2 r ∂θ r ∂θ ∂gr θ 2cot θ ∂gφr ∂gθ φ cot θ ∂gθr cot θ ∂gφφ + + 2 + 2 − − 2 r ∂r r sin θ ∂φ r sin θ ∂φ r ∂θ r ∂r cot θ gφφ gr θ gθr cot θ gθ θ cot θ∂gθ θ + 2 + 2 2 − + (A.12d) + r ∂r r r2 r2 r sin θ ∂ 2 gθr ∂ 2 gθ φ ∂grr ∂ 2 gθ θ ∂ 2 gθ θ 2∂gθ θ + 2 + − − 2 2 − 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂φ 2 r ∂r r ∂r cot θ ∂gθ φ cot θ ∂gφφ ∂gr φ 2cot θ∂gφθ + 2 − + 2 + 2 r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r sin θ ∂φ gθ θ cot θ gr θ gθ θ gφφ gφφ (A.12e) + 2 − 2 − 2 2 + 2 2 + r r r2 r sin θ r sin θ ∂ 2 gθ φ ∂ 2 gθ φ 2∂gθ φ ∂ 2 gθ θ ∂ 2 gθr 2∂gr φ − − − + − 2 r 2 sin θ∂θ∂φ r 2 ∂θ 2 r sin θ ∂r ∂φ ∂r 2 r ∂θ r ∂r cot θ ∂gφθ ∂gr θ cot θ ∂gθ φ cot θ ∂gφr cot θ∂gθ θ − + 2 − − − 2 2 2 r sin θ∂φ r ∂θ r sin θ ∂φ r ∂θ r ∂r gθ φ gφθ gθ φ cot θ gr φ gφθ (A.12f) − 2 + 2 + 2 2 + 2 2 − r r r2 r sin θ r sin θ ∂ 2 gφφ ∂ 2 gφθ ∂gφθ ∂ 2 gφr ∂ 2 gφr ∂gr φ + + 2 + − 2 2 − 2 2 2 r sin θ∂r ∂φ r sin θ ∂φ r ∂r ∂θ r ∂θ r ∂θ r ∂r ∂gφφ cot θ ∂gφr 2∂grr cot θ ∂gφθ 2cot θ∂gθr + 2 − − 2 + − 2 r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r sin θ ∂φ r ∂r gφr cot θ gθ φ gr φ cot θ∂gθ φ cot θ gφθ + (A.12g) + 2 + 2 2 + + r ∂r r r2 r2 r sin θ ∂ 2 gφr ∂ 2 gφφ ∂gr φ ∂ 2 gφθ ∂ 2 gφθ 2∂gφθ + + 2 − − − 2 2 2 2 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂φ r ∂θ r ∂r 2cot θ ∂gθ θ 2∂gr θ cot θ ∂gθ φ 2gθ φ cot θ∂gφφ − 2 − 2 + − 2 + 2 r sin θ∂φ r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r gφθ gθ φ 2cot θ gr φ + 2 2 + 2 2 + (A.12h) r2 r sin θ r sin θ ∂ 2 gφθ 2 r sin θ∂θ∂φ
∂ 2 gφr ∂gr θ ∂ 2 gφφ ∂ 2 gφφ + + 2 − 2 2 r ∂θ r sin θ ∂r ∂φ ∂r 2 r ∂θ ∂grr cot θ ∂gφθ cot θ ∂gφφ cot θ ∂gθ θ 2∂gφφ + − 2 − + − r ∂r r ∂r r sin θ ∂φ r 2 ∂θ r 2 ∂θ gφφ cot θ∂gθr gθ θ + 2 2 − 2 2 + r ∂r r sin θ r sin θ −
(A.12i)
We have implemented these equations in 1D, direction r , with only slight modifications to the strategy described for Cartesian geometry. Schematically, we represent the overall system of equations in the form ∂ AF(U ) ∂ H (U ) ∂U + + = G(q, r ) + S(q). ∂t ∂V ∂r
(A.13)
175
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
Here we distinguish between the area-weighted volumetric flux terms, AF, and the spatial flux terms H . The geometric source terms are represented by G(q, r ), and the plastic source terms are S(q). Note that in strict 1D-r flow, there is no angular dependence to any flow variable, and therefore terms proportional to cot θ (for example) vanish identically. As in the Cartesian case, we solve the time-centered edge Riemann problems to deduce ∗,n+1/2 single-valued time-centered edge states Ui+1/2 . These edge states are then used to construct the flux terms F and H , which are used to compute a preliminary update U˜ n+1 via the difference scheme ¢ ¢ 1t ¡ 1t ¡ ∗ ∗ ∗ ∗ − Ai−1/2 Fi−1/2 Ai+1/2 Fi+1/2 − Hi+1/2 − Hi−1/2 . U˜ in+1 = Uin − Vi 1ri
(A.14)
Next, we modify the preliminary update by inclusion of the geometric source terms. This is made second-order using a predictor-corrector strategy, ¡ ¢ 0 U˜ i = U˜ in+1 + 1t G qin , σin µ ¶n ∂σ (q˜ 0 − q n ) σ0 = σn + ∂q ¢ 1t ¡ ¡ n n ¢ 00 U˜ i = U˜ in+1 + G qi , σi + G(q˜ i0 , σi0 ) . 2
(A.15)
The plastic source terms are then evaluated at the half time step, giving the final result ¢ 1¡ n qi + q˜ i00 2 = U˜ i0 + 1t S(q¯ i ).
q¯ i = Uin+1
(A.16)
ACKNOWLEDGMENTS G.M. benefitted from correspondences with B. Plohr (SUNY at Stony Brook), V. Barcilon (University of Chicago), and M. Ortiz (Caltech).
REFERENCES 1. J. B. Bell, P. Colella, and M. L. Welcome, Conservative front-tracking for inviscid compressible flow, in AIAA 10th Comput. Fluid Dynamics Conf., Honolulu, Hawaii, 1991. 2. J. B. Bell, C. N. Dawson, and G. R. Shubin, An unsplit, higher order Godunov method for scalar conservation laws in multiple dimensions, J. Comput. Phys. 74, 1 (1988). 3. F. G. Blake, Spherical wave propagation in solid media, J. Acoust. Sci. Am. 24, 211 (1952). 4. I.-L. Chern and P. Colella, A Conservative Front Tracking Method for Hyperbolic Conservation Laws, Technical Report (LLNL, 1987), UCRL-97200. 5. P. Colella, A direct Eulerian MUSCL scheme for gas dynamics, SIAM J. Sci. Stat. Comput. 6, 104 (1985). 6. P. Colella, Multidimensional upwind methods for hyperbolic conservation laws, J. Comput. Phys. 87, 171 (1990). 7. P. Colella and P. R. Woodward, The piecewise parabolic method (PPM) for gas dynamical simulations, J. Comput. Phys. 54, 174 (1984). 8. S. K. Godunov, Symmetric form of magnetohydrodynamics equations, Chislennye Metody Mekhaniki Sploshnoi Sredy, Novosibirsk 3, 26 (1972), in Russian.
176
MILLER AND COLELLA
9. S. K. Godunov and E. I. Romensky, Thermodynamics, conservation laws, and symmetric forms of differential equations in mechanics of continuous media, in Computational Fluid Dynamics Review 1995 (Wiley, New York, 1995), pp. 19–31. 10. R. Hill, The Mathematical Theory of Plasticity (Oxford University Press, Oxford, U.K., 1950). 11. Mindy Fruchtman Lai, A Projection Method for Reacting Flow in the Zero Mach Number Limit, Ph.D. thesis (Dept. of Mechanical Engineering, Univ. of California, Berkeley, 1993). 12. E. H. Lee and D. T. Liu, Finite-strain elastic-plastic theory with application to plane-wave analysis, J. Appl. Phys. 38, 19 (1967). 13. J. Lubliner, Plasticity Theory (McMillan, New York, 1990). 14. J. Mandel, Plasticit´e Classique et Viscoplasticit´e (Springer-Verlag, New York, 1972). 15. G. H. Miller and E. G. Puckett, A high-order Godunov method for multiple condensed phases, J. Comput. phys. 128, 134 (1996). 16. R. B. Pember, J. B. Bell, P. Colella, W. Y. Crutchfield, and M. L. Welcome, An adaptive Cartesian grid method for unsteady compressible flow in irregular regions, J. Comput. Phys. 120, 278 (1995). 17. L. Petzold, Differential/algebraic equations are not ODEs, SIAM J. Sci. Stat. Comput. 3, 367 (1982). 18. J. E. Pilliod, A Second-Order Unsplit Method for Modeling Flames in Two-Dimensional Compressible Flow, Ph.D. thesis (University of California, Davis, 1996). 19. B. J. Plohr and D. H. Sharp, A conservative formulation for plasticity, Adv. Appl. Math. 13, 462 (1992). 20. K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. DeZeeuw, A solution-adaptive upwind scheme for ideal magnetohydrodynamics, J. Comput. Phys. 154, 284 (1999). 21. W. J. Rider, Filtering non-solenoidal modes in numerical simulations of incompressible flow, Int. J. Numer. Meth. Fluids 28, 789 (1998). 22. J. Saltzman, An unsplit 3D upwind method for hyperbolic conservation laws, J. Comput. Phys. 115, 53 (1994). 23. J. C. Simo and T. J. R. Hughes, Computational Inelasticity (Springer-Verlag, New York, 1997). 24. G. F. Smith and R. S. Rivlin, The strain-energy function for anisotropic elastic materials. Trans. Am. Math. Soc. 88, 175 (1958). 25. J. A. Trangenstein and P. Colella, A higher-order Godunov method for modeling finite deformation in elasticplastic solids, Comm. Pure Appl. Math. 44, 41 (1991). 26. B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys. 32, 101 (1979). 27. M. L. Wilkins, Calculation of elastic-plastic flow, in Methods in Computational Physics, edited by B. Alder, S. Fernbach, and M. Rotenberg (Academic, New York, 1964), Vol. 3, p. 211.
A High-Order Eulerian Godunov Method for Elastic–Plastic Flow in Solids1 G. H. Miller†, ‡ and P. Colella† †Applied Numerical Algorithms Group, Lawrence Berkeley National Laboratory, Berkeley, California 94720; and ‡Department of Applied Science, University of California, Davis, One Shields Avenue, Davis, California 95616 E-mail: [email protected]; [email protected] Received May 2, 2000; revised November 10, 2000
We present an explicit second-order-accurate Godunov finite difference method for the solution of the equations of solid mechanics in one, two, and three spatial dimensions. The solid mechanics equations are solved in nonconservation form, with the novel application of a diffusion-like correction to enforce the gauge condition that the deformation tensor be the gradient of a vector. Physically conserved flow variables (e.g., mass, momentum, and energy) are strictly conserved; only the deformation gradient field is not. Verification examples demonstrate the accurate capturing of plastic and elastic shock waves across approximately five computational cells. 2D and 3D results are obtained without spatial operator splitting. °c 2001 Academic Press Key Words: solid mechanics; shock waves; Godunov method; elasticity; plasticity.
1. INTRODUCTION
In this work, we present a high-order Godunov method for computing in Eulerian coordinates the multidimensional dynamics of elastic–plastic solids undergoing large deformations. Our approach is based on a new formulation of the equations of solid mechanics as a first-order system of hyperbolic partial differential equations (PDEs), a modification of that used by Trangenstein and Colella [25]. In [25], the usual conservation laws for mass, momentum, and energy, plus a constitutive model, are augmented by a form of equality of mixed partial derivatives that yields conservation equations for the entries of the inverse 1
Work at the Lawrence Berkeley National Laboratory was sponsored by the U.S. Department of Energy (DOE) Mathematical, Information, and Computing Sciences Division under Contract DE-AC03-76SF00098. Other work was supported by a subcontract from the Caltech Center for the Simulation of Dynamic Response in Materials, which in turn is supported by the Academic Strategic Alliances Program of the Accelerated Strategic Computing Initiative (ASCI/ASAP) under Subcontract B341492 of DOE Contract W-7405-ENG-48. 131 0021-9991/01 $35.00 c 2001 by Academic Press Copyright ° All rights of reproduction in any form reserved.
132
MILLER AND COLELLA
deformation gradient. This leads to equations of the form ∂U + ∇ · F(U ) = S(U ). ∂t
(1)
Here S(U ) contains source terms associated with the treatment of plasticity. These equations by themselves are not sufficient to specify the problem. In addition, we must impose linear constraints on the solution to guarantee that the inverse of the deformation gradient is in fact a gradient, i.e., that the curl of the rows of the deformation gradient vanish. These constraints can be written in the form L C (U ) = 0.
(2)
Here L C is a system of linear differential operators with constant coefficients. The constraint equation is an initial-value constraint: if (1) is satisfied, and L C (U ) is identically zero at some initial time, then L C (U ) vanishes identically for all later times. The constraint (2) plays an essential role in the analysis of the characteristic structure of the system (1). To get the physically correct eigenvectors and eigenvalues from the quasilinear form of the equations, one must use the constraint to replace some of the spatial derivatives. In general, solutions to (1), without imposing (2), give rise to unphysical wave propagation properties, even for linearized waves as was observed in [25]. A difficulty arises when one attempts to compute solutions to (1) and (2) using a conservative finite difference method. To the extent that a modified equation analysis is valid, we expect the behavior of the numerical solution to behave very similarly to the solution to the following system of PDEs: ∂U Mod + ∇ · F(U Mod ) = S(U Mod ) + τU (U Mod ) ∂t L C (U Mod ) = τC (U Mod ).
(3)
Here τU and τC are truncation error terms, which are nonzero. In general, these terms, and in particular τC , cannot be eliminated. The practice of enforcing a discretized form of the constraint (2) at the end of each time step using a Hodge projection would guarantee that a discretized form of (2) is satisfied identically. However, that will change the form of τC , but not set it to zero. The observation that the truncation error terms are a small perturbation to the equations is not sufficient to guarantee that U Mod is close to U . There is much less known about the well-posedness of systems of equations that are combinations of evolution equations and constraints than there is about pure evolution equations, and unexpected pathologies are known to occur [17]. The approach we want to take on this problem starts with an analysis due to Godunov [8, 9]. Numerical methods based on this approach have been recently investigated for the MHD equations in [20], the case for which Godunov first applied this analysis. Godunov modifies (1) in the following way: ∂U + ∇ · F(U ) = S(U ) + ξ L C (U ). ∂t
(4)
Here ξ = ξ(U ) can be chosen so that the system has the physically correct linearized eigenstructure, independent of whether L C vanishes. In addition, L C satisfies a transport
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
133
equation such that if L C (U ) is identically zero at some time, then it remains so for all later times. The numerical method we present here is based on the form of the equations given by (4). Thus we are discretizing a well-posed initial value problem without constraints, independent of whether the constraint (2) is satisfied. This gives us a high degree of confidence that a stable and consistent method can be developed. Of course, the extent to which we compute a solution to the original physical problem (1) and (2) depends strongly on whether the constraint is satisfied, but now that is purely an accuracy issue, without any impact on the stability of the method. In fact, we will investigate the use of various methods of limiting discrete measures of L C , similar to filtering methods developed for incompressible flow [11, 21]. Examples that demonstrate the method for elastic–plastic deformations of a homogeneous solid domain are presented. Aside from some simple one-dimensional problems that involve free surfaces, methods for handling material interfaces (contact discontinuities) are not described here. To address more general engineering problems, including those concerned with fluid–solid coupling, our intent is to combine this solid mechanics method with a volume-of-fluid interface treatment analogous to [15]. 2. GOVERNING EQUATIONS
The mechanical behavior of solids is described by observable variables (e.g., density ρ, momentum ρv, internal energy E, and the deformation F with respect to a chosen reference state) and also unobservable internal parameters which describe the response of the material to deviatoric stress. One constitutive representation of this behavior is through the multiplicative decomposition of the total deformation into elastic and inelastic components [12], F = F eF p.
(5)
Here F is the Lagrangian coordinate deformation which relates the spatial coordinate frame x = x(a, t) to the material coordinate frame a: Fαβ =
∂ xα . ∂aβ
(6)
We refer to F p as the plastic deformation tensor, although the numerical scheme we will present applies to more general inelastic deformations. According to (5), F p is a fictitious state of total deformation in which there is no elastic deformation: given an initial total deformation F, and a purely elastic relaxation path F e → I , the total observable deformation will evolve to F p , F → F p . The state F p is a function of the deformation history of the material. We represent this history through a single scalar parameter κ, a work-hardening measure, and constitutive flow rules F˙ p = h(ρ, g, F p , E, κ) κ˙ = K(ρ, g, F p , E, κ), which depend on the state variables but not their gradients.
(7) (8)
134
MILLER AND COLELLA
The equations of solid mechanics are then given by 0 ρvα ρ ρf ρvvα − σ eα ρv ρ(8 + v· f) ρ Ev − v σ α β βα ρE (v × (∇ × g T ))T e ge x gvδxα x T T (v × (∇ × g )) e ge gvδ yα y y ∂ ∂ = , gez + T T ∂t gvδz α (v × (∇ × g )) ez ∂ xα p ρF ex ρF p ex vα ρhex ρF p e p y ρF e v y α ρhe y p p ρF ez ρF ez vα ρhez ρκ
ρκvα
(9)
ρK
where ex , e y , and ez are the Cartesian unit vectors, and E is the sum of internal energy and kinetic energy (E = E + 12 v · v). For generality, we include a heat source term 8 and a body force vector f . The system of equations (9) is abbreviated ∂U ∂ Fα (U ) = S(U ), + ∂t ∂ xα
(10)
where U is the vector of quasi-conservation-form variables (ρ, ρv, ρ E, etc. . . .), Fα (U ) is the flux in direction eα , and where S(U ) is the vector of source terms. Here we follow the treatment in [25] and use the inverse deformation gradient g = F −1 as dependent variables. However, we introduce additional nonconservative terms in the evolution equations of g. We will show below that the addition of these terms leads to a well-behaved hyperbolic structure for the equations, independent of whether the curl of g T vanishes. However, we note here that G = ∇ × g T satisfies the following evolution equation: ∂G + ∇ · (vG − Gv) = 0. ∂t
(11)
In particular, if G vanishes identically at time t = 0, it vanishes at all later times. To solve these equations we adopt a predictor–corrector strategy. For each time step, we first solve the conservative flux differencing left-hand side of (10) using fluxes derived (by solution to Riemann problems) from edge- and time-centered variables that include timecentered contributions from the source terms. The solution obtained by flux differencing is then modified by addition of the source terms, evaluated using time-centered and cellcentered variables, and acting over the full time step 1t. The solution to the flux differencing equations is based upon the standard high-order Godunov strategy. This strategy begins with a characteristic analysis of the equations, which makes use of the linearized 1D equations in direction eα , ρ
v E gex ge y ∂ gez + A ∂ ∂t F p ex ∂ xα p F ey p F ez κ σ eα
ρ
0
v f E 8 gex 0 ge y 0 ge z = 0 , p he F ex x p he F ey y p F ez hez κ σ eα
K bα
(12)
135
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
where
vα
0 0 0 0 A= 0 0 0 0 0 0
ρeαT
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
−(σ eα )T /ρ vα
0
0
0
0
0
0
0
gδxα
0
vα I
0
0
0
0
0
0
gδ yα
0
0
vα I
0
0
0
0
0
gδzα
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα I
0
0
0
0
0
0
0
0
0
vα
−Aαα
0
0
0
0
0
0
0
0
0
−I /ρ 0 0 0 0 (13) 0 0 0 0 vα I
with Aαβ = −
∂σ eα g ∂geβ
(14)
and bα =
∂σ eα ∂σ eα ∂σ eα :h+ K+ 8. ∂F p ∂κ ∂E
(15)
The eigenvalue decomposition of A uses the technique of eigenvalue deflation and hinges upon recognition of the matrices Aαα as being acoustic wave propagation tensors for waves traveling in direction eα , ρ u¨ β = (Aαα )γ δ
∂ 2uδ , ∂ xβ ∂ xγ
(16)
where u is the displacement vector. The matrices Aαα are positive definite as a requirement of thermodynamic stability. This is made clear by writing Aαα in terms of gradients of the spatial displacements uˆ defined relative to the current configuration, (Aαα )βγ = ρ
∂ 2E . ∂ uˆ βα ∂ uˆ γ α
(17)
Here, uˆ βα is related to the deformation tensor Fβα with the reference coordinate frame {a} chosen to correspond to the current spatial frame {x}: uˆ βα = Fβα |{a}={x} − δβα .
(18)
Aαα is therefore a component of the Hessian of E, which is positive definite for a
136
MILLER AND COLELLA
thermodynamically stable material, and consequently Aαα has positive real eigenvalues and three linearly independent eigenvectors. Recognizing Aαα as being the acoustic wave propagation tensor suggests the wave equation solution 2 , Aαα X ac = ρ X ac 3ac
(19)
where 3ac is the diagonal matrix of acoustic wave speeds c, 3ac = diag(c1 , c2 , c3 ), and X ac are the acoustic displacement vectors. The linearized 1D matrix A then has eigenvalue decomposition A = X 3X −1 ,
(20)
with X , the matrix of right eigenvectors, given by
1 0 0 0 0 0 0 0 0
0 0 0 0 X = 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 I
0 0 0 0 0 0
0 0 I
0 0 0 0 0
0 0 0 I
0 0 0 0
0 0 0 0 I
0 0 0
0 0 0 0 0 I
0 0
0 0 0 0 0 0 I
0
0 0 0 0 0 0 0 1
−ρeαT X ac
−ρeαT X ac
−X ac 3ac T T (σ eα ) X ac /ρ (σ eα ) X ac /ρ −g X ac δxα −g X ac δxα −g X ac δ yα −g X ac δ yα −g X ac δzα −g X ac δzα , 0 0 0 0 0 0 0 0 X ac 3ac
2 ρ X ac 3ac
0 0 0 0 0 0 0 0 0
(21)
2 X ac 3ac ρ
and 3, the diagonal matrix of eigenvalues, given by
0 vα 0 0 vα 0 0 0 vα I 0 0 0 0 0 0 3= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 vα I 0 0 0 0 0 0 0
0 0 0 0 vα I 0 0 0 0 0 0
0 0 0 0 0 vα I 0 0 0 0 0
0 0 0 0 0 0 vα I 0 0 0 0
0 0 0 0 0 0 0 vα I 0 0 0
0 0 0 0 0 0 0 0 vα 0 vα I 0
0 0 0 0 0 0 0 0 0 − 3ac 0 vα I
0 0 0 0 0 0 . (22) 0 0 0 0 + 3ac
The wave speeds are Galilean invariant and properly analogous to the Lagrangian
137
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
representation, with 3 − waves with velocities vα − cγ , 3 + waves with velocities vα + cγ , and 21 material waves with speeds vα . X −1 , the inverse of X , is given by
1
0
0 0 0 0 0 0 0 0 −1 X = 0 0 0 0 0 0 0 0 1 −1 −1 0 2 3ac X ac −1 −1 0 − 12 3ac X ac
0 0 0 0 0 0 0 0
−2 −1 eαT X ac 3ac X ac
−2 −1 1 0 0 0 0 0 0 0 −(σ eα )T X ac 3ac X ac /ρ 2 −2 −1 0 I 0 0 0 0 0 0 g X ac 3ac X ac δxα /ρ −2 −1 X ac δ yα /ρ 0 0 I 0 0 0 0 0 g X ac 3ac −2 −1 0 0 0 I 0 0 0 0 g X ac 3ac X ac δzα /ρ . (23) 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 1 −2 −1 0 0 0 0 0 0 0 0 3 X 2ρ ac ac 1 −2 −1 0 0 0 0 0 0 0 0 3 X 2ρ ac ac 3. NUMERICAL METHOD: 1D
In 1D we discretize space into cells, indexed with subscript i, with width 1xi . Time is discretized in steps of 1t with integer superscript index n; t n+1 − t n = 1t. The generalization to 2D and 3D is similar, with indices j and k used for the second and third dimensions, respectively. Half-integral subscript indices represent edge-centered quantities. Lowercase Greek subscripts are used to denote vector and tensor indices. We begin by evaluation of the equation of state in each cell to determine the Cauchy stress σ , the acoustic wave propagation tensor Aαα , and the thermodynamic derivatives ∂σ/∂E|g,F p ,κ , ∂σ/∂g|E,F p ,κ , ∂σ/∂F p |E,g,κ , and ∂σ/∂κ|E,g,F p . Next, we evaluate the 1D slopes dq/d xα of the 27 primitive cell-centered variables q, q = (ρ, v, E, g, F p , κ, σ eα ).
(24)
We construct these slopes beginning with the van Leer slope in cell i which uses the monotonized limiter [26]: µ
∂q ∂x
¶vL i
µ
¶ 2|qi+1 − qi−1 | 2|qi − qi−1 | 2|qi+1 − qi | = sign(qi+1 − qi−1 ) min , , . 1xi−1 + 21xi + 1xi+1 1xi 1xi (25)
A fourth-order-accurate slope is then constructed as [5] µ
∂q ∂x
¶4th = i
³h ¡ ¢vL i h ¡ ¢vL i´ qi+1 − 14 1xi+1 ∂q − qi−1 + 14 1xi−1 ∂q ∂ x i+1 ∂ x i−1 1 1xi−1 4
+ 1xi + 14 1xi+1
.
(26)
138
MILLER AND COLELLA
To prevent overshoot and ringing, dissipation at strong shocks may be introduced via a “flattening parameter” χ, 0 ≤ χ ≤ 1, whence [5–7] µ
∂q ∂x
¶
µ = χi
i
∂q ∂x
¶4th .
(27)
i
The determination of this flattening parameter is described in a later section. These limited slopes are used to construct time-centered, edge-valued estimates of the primitive variables. The exact solution of the linearized equations, which we abbreviate as ∂q ∂q +A = s, ∂t ∂x
(28)
gives time-centered edge values n+1/2 q R,i−1/2 n+1/2
q L ,i+1/2
µ µ ¶ ¶ 1t 1xi 1t −1 ∂q Xi si = − 3i + I X i + 2 1xi ∂x i 2 µ µ ¶ ¶ 1t 1xi 1t −1 ∂q n = qi − 3i − I X i + Xi si . 2 1xi ∂x i 2 qin
(29a) (29b)
However, this construction uses both upwind and downwind characteristics. We make the method strictly upwind by filtering out the downwind characteristics: n+1/2 q R,i−1/2 n+1/2
q L ,i+1/2
µ µ ¶ ¶ 1t 1xi 1t −1 ∂q X i P− si = − 3i + I X i + 2 1xi ∂x i 2 µ µ ¶ ¶ 1t ∂q 1xi 1t X i P+ si = qin − 3i − I X i−1 + 2 1xi ∂x i 2 qin
(30a) (30b)
with projection operators P± defined as µ
µ P−
µ
µ P+
1t 3+ I 1x 1t 3− I 1x
(¡ 1t
¶¶ αβ
=
¶¶ αβ
1x
0 (¡ 1t
=
1x
¢ 3αα + 1 δαβ ¢ 3αα − 1 δαβ
3αα ≤ 0 3αα > 0 3αα ≥ 0 3αα < 0.
0
(31a)
(31b)
At each cell edge (i + 1/2), time-centered values are thus obtained from the left (i) and right (i + 1) neighboring cells. These edge values are then used to pose a Riemann problem: n+1/2 an initial value problem with constant left and right initial states given by q L ,i+1/2 and n+1/2 q R,i+1/2 , respectively. We approximate the solution to the Riemann problem by decomposing n+1/2 n+1/2 the jump q R,i+1/2 − q L ,i+1/2 in terms of the eigenvectors X of the linearized coefficients A. Specifically, n+1/2
n+1/2
q R,i+1/2 − q L ,i+1/2 =
27 X
ϕγ X γ ,i+1/2 ,
(32)
γ =1
where eigenvector column X γ ,i+1/2 is evaluated with certain L or cell-i properties if 3γ γ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
139
is a member of the − family (i.e., of the form veα − c), or with certain R or cell-(i + 1) properties if 3γ γ is a member of the + family (of the form veα + c), as given by the discretization
1 0 0 0 0 0 0 0 0
X i+1/2
0 0 0 0 = 0 0 0 0 0 0
−ρi eαT X ac,i
0 0 0 0 0 0 0 0
X ac,i 3ac,i
1 0 0 0 0 0 0 0
(σi+1/2,L eα )T X ac,i ρi
0
I
0 0 0 0 0 0 −gi+1/2,L X ac,i δxα
0 0
I
0 0 0 0 0 −gi+1/2,L X ac,i δ yα
0 0 0 0 0 0
0 0 0 0 0 0
I 0 0 0 0 0
0 0 0 0 0 0
0 I 0 0 0 0
0 0 I 0 0 0
0 0 0 I 0 0
0 0 0 0 1 0
−gi+1/2,L X ac,i δzα 0 0 0 0 2 ρi X ac,i 3ac,i
−ρi+1 eαT X ac,i+1
−X ac,i+1 3ac,i+1 (σi+1/2,R eα )T X ac,i+1 ρi+1 −gi+1/2,R X ac,i+1 δxα −gi+1/2,R X ac,i+1 δ yα . −gi+1/2,R X ac,i+1 δzα 0 0 0 0 2 X ac,i+1 3ac,i+1 ρi+1 (33)
In this expression, the density ρ, and the components (X ac , 3ac ) of the acoustic propagation tensor, are evaluated at the cell centers to avoid multiple evaluations of the equation of state. From the coefficients ϕγ of the jump decomposition, the material velocity v ∗ · eα at the cell edge is determined by adding to the L state the contributions of the − family, or by subtracting from the R state the contributions of the + family, that is, ∗ · eα = v ∗L ≡ vi+1/2,L · eα + ϕ6 X 6β,i+1/2 + ϕ7 X 7β,i+1/2 + ϕ8 X 8β,i+1/2 vi+1/2
(34)
or ∗ vi+1/2 · eα = v ∗R ≡ vi+1/2,R · eα − ϕ9 X 9β,i+1/2 − ϕ10 X 10β,i+1/2 − ϕ11 X 11β,i+1/2 , (35)
where β = 2, 3, 4 for directions eα equal to ex , e y , ez , respectively. We average the results of these calculations to determine the normal-direction edge velocity v ∗ · eα , ∗ · eα = vi+1/2
1 ∗L (v + v ∗R ). 2
(36)
For other properties to be evaluated at the cell edge as solutions to the Riemann problem we do not average the values evaluated from the L and the R states as above. Instead, we evaluate from the L state if v ∗ · eα is positive, or from the R state if v ∗ · eα is negative. Only if v ∗ · eα is approximately zero do we average these estimates. The evaluations include only
140
MILLER AND COLELLA
upwind characteristics by writing
∗ qi+1/2
P qi+1/2,L + γ w L ,γ ϕγ X γ P = qi+1/2,R + γ w R,γ ϕγ X γ ¢ P 1 ¡q i+1/2,L + qi+1/2,R + γ (w L ,γ + w R,γ )ϕγ X γ 2
∗ vi+1/2 · eα > ² ∗ vi+1/2 · eα < −² (37) ¯ ¯ ∗ ¯ ¯v i+1/2 · eα ≤ ²
with ( w L ,γ = ( w R,γ =
1
∗ 3γ γ ,i − vi · eα + vi+1/2 · eα < −²
0
otherwise
1
∗ · eα > ² 3γ γ ,i+1 − vi+1 · eα + vi+1/2
0
otherwise.
(38a)
(38b)
The value of w L ,γ is 1 when eigenvalue γ , estimated using the ∗ value of the material velocity together with the i cell-centered acoustic wave speeds, is negative and w L ,γ is 0 otherwise; w R,γ is 1 when the approximated value of eigenvalue γ is positive and 0 otherwise. In our computations presented below, we use a value of ² = 10−9 . By this procedure, we obtain the edge ∗ value solutions of the Riemann problem, ρ ∗ , v ∗ , E ∗ , g ∗ , F p∗ , κ ∗ , and (σ e j )∗ . These are then used to compute edge-valued fluxes (cf. (9, 10)). For example, in direction ex ,
Fx,i+1/2
ρvx
∗
ρvx2 − σx x ρv y vx − σ yx ρvz vx − σzx ρ Ev − v σ − v σ − v σ x x xx y yx z zx vx gex + v y ge y + vz gez = . 0 0 p ρvx F ex p ρvx F e y ρvx F p ez ρvx κ i+1/2
(39)
In 1D we obtain a preliminary update U˜ of the variables U by conservatively differencing the fluxes: ¢ 1t ¡ ∗ ∗ Fi+1/2 − Fi−1/2 . U˜ in+1 = Uin − 1xi
(40)
The final time-(n + 1) value of the variables U is obtained from the preliminary values U˜ by addition of the source terms S: Uin+1 = U˜ in+1 + 1t Si .
(41)
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
141
We discretize this in the general 3D (Cartesian) case as
0
n+1/2 n+1/2 ρi fi ¢ n+1/2 ¡ n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 n+1/2 ρi + vx,i f x,i + v y,i f y,i + vz,i f z,i 8 ∗ ∗ ∗ ∗ ¡ ¢ © gi, j+1/2,k ex − gi, j−1/2,k ex n+1/2 gi+1/2, j,k e y − gi−1/2, j,k e y − v y,i jk 1xi 1y j ∗ ∗ ∗ ∗ ¡ ¢ª g e − g e g e − g e n+1/2 i, j,k+1/2 x i, j,k−1/2 x i+1/2, j,k z i−1/2, j,k z − − vz,i jk 1z k 1xi ∗ ∗ ¢ © n+1/2 ¡ gi,∗ j+1/2,k ez − gi,∗ j−1/2,k ez gi, j,k+1/2 e y − gi, j,k−1/2 e y − v z,i jk 1y j 1z k ∗ ∗ ∗ ∗ ¡ ¢ª gi, j+1/2,k ex − gi, j−1/2,k ex n+1/2 gi+1/2, j,k e y − gi−1/2, j,k e y n+1 n+1 . ˜ − − v Ui = Ui + 1t x,i jk 1xi 1y j ∗ ∗ ¢ © n+1/2 ¡ gi,∗ j,k+1/2 ex − gi,∗ j,k−1/2 ex gi+1/2, j,k ez − gi−1/2, j,k ez − vx,i jk 1z k 1xi ∗ ∗ gi,∗ j,k+1/2 e y − gi,∗ j,k−1/2 e y ¢ª n+1/2 ¡ gi, j+1/2,k ez − gi, j−1/2,k ez − − v y,i jk 1y j 1z k n+1/2 n+1/2 ρi jk h i jk ex n+1/2 n+1/2 ρi jk h i jk e y n+1/2 n+1/2 ρi jk h i jk ez n+1/2
ρi jk
n+1/2
K i jk
(42) In 1D we use the 3D discretization above but retain only terms in ∂/∂ x and ∂ 2 /∂ x 2 and omit derivatives in all transverse directions. n+1/2 In the above expression, time-centered terms (e.g., ρi jk ) are estimated with n+1/2
qi jk
≈
¢ 1¡ n q + q˜ in+1 jk , 2 i jk
(43)
except for the g ∗ s appearing in the (v × G) terms. These are obtained at the half time step and cell edges as components of the Riemann problem solutions.
4. NUMERICAL METHOD: 2D AND 3D
To extend the 1D method described above to multiple spatial dimensions, we use a spatially unsplit fully corner-coupled second-order-accurate scheme after [6, 22]. In 2D, this predictor–corrector approach begins by estimating the 1D x and y fluxes at each cell edge, using the higher order 1D approach described in the previous section. These predictor fluxes, F˜ x and F˜ y , are given schematically as solutions to the Riemann problem R as ¡ ¡ n+1/2 ¢¢ x n+1/2 F˜ i+1/2, j = Fx R qx L ,i+1/2, j , qx R,i+1/2, j ¡ ¡ n+1/2 ¢¢ y n+1/2 F˜ i, j+1/2 = Fy R q y L ,i, j+1/2 , q y R,i, j+1/2 .
(44a) (44b)
142
MILLER AND COLELLA
The predictor fluxes are used to pose a corrector problem, wherein the edge values are augmented by transverse predictor fluxes. Schematically, ¡ ¡ n+1/2 ¢¢ n+1/2 Fx,i+1/2, j = Fx R q 0x L ,i+1/2, j , q 0x R,i+1/2, j ¡ ¡ n+1/2 ¢¢ n+1/2 Fy,i, j+1/2 = Fy R q 0y L ,i, j+1/2 , q 0y R,i, j+1/2
(45a) (45b)
with, for example, ¢ 1t ¡ ˜ y y F i, j+1/2 − F˜ i, j−1/2 21y j X ¡ y ¡ n ¢T ¢T 1t y + 0(γ , δ). v × e y × g˜ i, j+1/2 − g˜ i, j−1/2 γδ 21y j γ δ i, j
q 0 x L ,i+1/2, j = qx L ,i+1/2, j − n+1/2
n+1/2
(46)
The vector 0 is introduced to align the elements of the matrix (v × ∇ × g T )T with the appropriate elements of the vector q: 0(r, s) = (0, 0, 0, 0, 0, δγ 1 δδ1 , δγ 2 δδ1 , δγ 3 δδ1 , δγ 1 δδ2 , δγ 2 δδ2 , δγ 3 δδ2 , δγ 1 δδ3 , δγ 2 δδ3 , δγ 3 δδ3 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)T .
(47)
In 2D there are therefore four Riemann problems solved per cell: two in the predictor and two in the corrector steps. In setting up the corrector step, the components ρ, v, E, g, F p , and κ of the vectors q 0 are updated as indicated in Eq. (46). Our 1D Riemann solver also requires time-centered edge values of the stresses, (σ eα ) L/R in direction eα , and these components of q 0 are calculated by updating the (σ eα ) L/R components of q L/R with the change in stress accompanying the 0 − q L/R in E, g, F p , and κ using cell-centered thermodynamic derivatives. changes q L/R For example, (σ 0 ex ) L ,i+1/2, j = (σ ex ) L ,i+1/2, j +
µ
¯ ¶n ¡ 0 ¢ ∂σ ex ¯¯ E L ,i+1/2, j − E L ,i+1/2, j ¯ ∂E g,F p ,κ i j
¯ ¶n ¡ 0 ¢ ∂σ ex ¯¯ κ L ,i+1/2, j − κ L ,i+1,2, j + ¯ ∂κ g,F p ,E i j ¯ ¶n X µ ∂σ ex ¯ ¡ ¢ ¯ (gγ δ )0L ,i+1/2, j − (gγ δ ) L ,i+1/2, j + ¯ ∂gγ δ g6=gγ δ ,F p ,κ,E i j γδ µ ¶n X ∂σ ex ¯¯ ¡ p¢ ¡¡ p ¢0 ¢ ¯ + Fγ δ L ,i+1/2, j − Fγ δ L ,i+1/2, j . (48) p ¯ p ∂F γ δ g,F p 6=Fγ δ ,κ,E i j γδ µ
By employing this approximation we require only one equation-of-state evaluation per time step per cell for problems involving only elasticity. In problems that also include plasticity, additional equation-of-state evaluations are required for the computation of plastic source terms. In 3D there are two corrector steps: first, ¡ ¡ n+1/2 ¢¢ x n+1/2 F˜ i+1/2, j,k = Fx R qx L ,i+1/2, j,k , qx R,i+1/2, j,k
(49a)
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
¡ ¡ n+1/2 ¢¢ y n+1/2 F˜ i, j+1/2,k = Fy R q y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ n+1/2 ¢¢ z n+1/2 F˜ i, j,k+1/2 = Fz R qz L ,i, j,k+1/2 , qz R,i, j,k+1/2 ,
143
(49b) (49c)
and then ¡ ¡ ¢¢ n+1/2 x|y 0(y) n+1/2 F˜ i+1/2, j,k = Fx R q 0(y) x L ,i+1/2, j,k , q x R,i+1/2, j,k ¡ ¡ ¢¢ n+1/2 x|z 0(z) n+1/2 F˜ i+1/2, j,k = Fx R q x0(z) L ,i+1/2, j,k , q x R,i+1/2, j,k ¡ ¡ ¢¢ n+1/2 y|x 0(x) n+1/2 F˜ i, j+1/2,k = Fy R q 0(x) y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ ¢¢ n+1/2 y|z 0(z) n+1/2 F˜ i, j+1/2,k = Fy R q 0(z) y L ,i, j+1/2,k , q y R,i, j+1/2,k ¡ ¡ ¢¢ n+1/2 z|x 0(x) n+1/2 F˜ i, j,k+1/2 = Fz R q 0(x) z L ,i, j,k+1/2 , q z R,i, j,k+1/2 ¡ ¡ ¢¢ n+1/2 z|y 0(y) n+1/2 F˜ i, j,k+1/2 = Fz R q 0(y) z L ,i, j,k+1/2 , q z R,i, j,k+1/2 ,
(50a) (50b) (50c) (50d) (50e) (50f)
with, e.g., ¢ 1t ¡ ˜ y y F i, j+1/2,k − F˜ i, j−1/2,k 31y j X ¡ ¢T ¢T ¡ y 1t y + vi,n j,k × e y × g˜ i, j+1/2,k − g˜ i, j−1/2,k 0(γ , δ). (51) γδ 31y j γ δ
n+1/2
q 0(y) x L ,i+1/2, j,k = q x L ,i+1/2, j,k − n+1/2
The final fluxes, which enter the conservative differencing step of the integration, are then computed as ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 (52a) Fx,i+1/2, j,k = Fx R qx L ,i+1/2, j,k , qx R,i+1/2, j,k ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 Fy,i, j+1/2,k = Fy R q y L ,i, j+1/2,k , q y R,i, j+1/2,k (52b) ¡ ¡ 00 n+1/2 ¢¢ 00 n+1/2 Fz,i, j,k+1/2 = Fz R qz L ,i, j,k+1/2 , qz R,i, j,k+1/2 (52c) with, e.g., 00 n+1/2
¢ ¢ 1t ¡ ˜ y|z 1t ¡ ˜ z|y y|z z|y F i, j+1/2,k − F˜ i, j−1/2,k − F i, j,k+1/2 − F˜ i, j,k−1/2 21y j 21z k ¡ y|z ¢T ¢T 1t X¡ n y|z + 0(γ , δ) vi, j,k × e y × g˜ i, j+1/2,k − g˜ i, j−1/2,k γδ 21y j γ δ n+1/2
qx L ,i+1/2, j,k = qx L ,i+1/2, j,k −
+
¡ z|y ¢T ¢T 1t X¡ n z|y × ez × g˜ i, j,k+1/2 − g˜ i, j,k−1/2 0(γ , δ). v γδ 21z k γ δ i, j,k
(53)
There are a total of 12 Riemann solves per unit cell in 3D: 9 in the predictor steps and 3 in the corrector step. The σ components of the vectors q 0 and q 00 are computed as in the 2D case. 5. PLASTIC SOURCE TERMS
We present here an associated plasticity evolution equation for the rate of change of the plastic deformation tensor F p with time. The more common approach (e.g., [19, 23])
144
MILLER AND COLELLA
is to consider evolution equations for the plastic strain η p = 12 (F pT F p − I ), the plastic Green tensor C p = F pT F p , or the plastic Finger tensor b p = F p F pT . We choose instead to evolve the full nine-component plastic deformation tensor F p . This choice is necessary to be capable of modeling arbitrary crystal systems (see, e.g., [24]). For example, the elastic response of the lowest symmetry crystal system (triclinic) depends upon all six components of the elastic Green tensor. If one were to specify the total inverse deformation g, and either η p , C p , or b p , then all six components of C e could not be determined. Although our examples will make use of isotropic equation of state models (whose elastic invariants may be determined using g and C p ), our goal is to construct a framework of more general applicability. To motivate our choice of evolution equations for F p we begin by postulating the existence of a hyperelastic equation of state, ˆ F p , κ, S), E = E(g,
(54)
where S is the specific entropy. The material derivative of E is ∂E ∂E ˙ 1 ∂E ˙ p ∂E E˙ = κ˙ + S = − σβγ Fγ α g˙ αβ + 8, g˙ αβ + p Fαβ + ∂gαβ ∂κ ∂S ρ ∂Fαβ
(55)
where the second equality equates energy change with the sum of work and heat. Solving for entropy production (dissipation) we have 0 ≤ S˙ = − =
1 T
µ
¶ 1 ∂E ˙ p 1 ∂E σβγ Fγ α ∂E ∂S − + κ˙ + 8 g˙ αβ − p F ρ ∂gαβ T ∂Fαβ αβ T ∂κ ∂E
1 (9plast + 9therm ). ρT
(56)
Here ∂E/∂S = T is the temperature, and we have introduced the specific power of thermal dissipation, 9therm = ρ8,
(57)
and the specific power of plastic dissipation,
9plast = −ρ
∂E ˙ p ∂E − ρ κ˙ p F ∂κ ∂Fαβ αβ
e ˙ Fαβ − ϑ κ˙ = gβγ σγ δ Fδα p
p ˙p Fαβ − ϑ κ˙ = gβγ σγ δ Fδν gνα
= 6 : L p − ϑ κ˙
(58)
p ˙ ˙ Fαβ = (F p )−1 with ϑ = ρ∂E/∂κ being the work-hardening modulus. L νβ = gνα να Fαβ is the p
p
p
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
145
plastic distortion rate [14], and 6βν = gβγ σγ δ Fδν is the thermodynamic force conjugate to L p . The dependence of S˙ on g˙ αβ vanishes because σαβ = −ρ
∂E gγ α . ∂gγβ
(59)
In evaluating (58) we have assumed that E depends on g and F p only through the elastic deformation F e = F g p = (F p g)−1 , e.g., ˜ e , S, κ), E = E(F
(60)
¯ ∂E ¯¯ 1 e = − gβγ σγ δ Fδα . p ρ ∂Fαβ ¯S,κ
(61)
whence
Thermodynamics requires that the internal energy depends upon the volume, and we assume by (60) that this energy dependence is carried by the tensor F e , e.g., V = V0 det F e . For this to be true, it is necessary that det F p = 1 at all times (i.e., V = V0 det F = V0 det F e det F p ; V = V0 det F e iff det F p = 1.) Therefore, (60) assumes that plastic flow is volume-preserving. We postulate a plastic yield surface f = 0, which we represent for illustrative purposes with a Mises–Huber constitutive model written in terms of the Cauchy stress σ , a constant yield stress parameter σY , and the work-hardening modulus ϑ: r f (σ, ϑ) = kdev σ k −
2 (σY + ϑ). 3
(62)
Here, dev σ = σ − 13 (tr σ )I is the stress deviator, and kAk is the Schur norm of A, kAk2 = Aαβ Aαβ = tr(A T A). The flow model we adopt is derived from (62) by the postulate of maximum plastic dissipation [10, 13]. The plastic dissipation (58) is considered as a function of the variables 6 ˙ 9plast = 9plast (6, ϑ; L p , κ). ˙ The plastic dissipation and ϑ, with fixed parameters L p and κ; is then maximized with respect to 6 and ϑ, subject to the constraint that f = 0 during plastic flow. The resulting flow laws are dev(σ ) F˙ p = ζ F p g F kdev(σ )k r 2 κ˙ = ζ 3
(63) (64)
with ζ a parameter chosen to satisfy the Kuhn–Tucker complementarity conditions and the “consistency condition” [23] f =0
(65)
ζ ≥0
(66)
ζf = 0
(67)
ζ f˙ = 0 (if f = 0) .
(68)
146
MILLER AND COLELLA
The flow model (63) is consistent with the assumption that plastic flow is volumepreserving, p p (det˙F p ) = (det F p )gνδ F˙ δν µ ¶ (dev σ )βγ p p = ζ (det F p )gνδ Fδα gαβ Fγ ν kdev σ k
(dev σ )αα kdev σ k = 0 because tr(dev σ ) = 0,
= ζ (det F p )
(69)
and is therefore compatible with the assumption made in evaluating 9plast (58). As an example, we use a modified Mooney–Rivlin equation of state, ρ0 E(C e , S) =
λ(S) √ µ(S) µ(S) (ln det C e )2 + tr C e − log det C e 2 2 2 µ ¶ 1 ρ0 ϑ0 κ + e−ϑ1 κ + ρ ϑ1
(70)
where C e is the elastic Green tensor, C e = F eT F e .
(71)
This equation of state gives a work-hardening modulus, ϑ(κ) = ρ
∂E = ϑ0 (1 − e−ϑ1 κ ), ∂κ
(72)
in terms of two parameters: ϑ0 is the ultimate, asymptotic value of the work-hardening modulus, and ϑ1 dictates the rate of approach of the asymptotic limit. The combined elastic–plastic evolution problem is solved with a predictor–corrector strategy. The inverse total deformation g is advanced in accordance with the equations of motion, with the plastic deformation F p being conservatively advected. This step may predict a coordinate in state space that lies outside the convex manifold of permissible states f (σ, ϑ) ≤ 0, in which case a plastic corrector step is used to bring state back to the yield surface. The algorithmic approach is a return mapping algorithm [23], modified to require only one equation-of-state evaluation. Begin the iteration sequence with iteration index m = 0, F p(0) = F p,n+1 κ (0) = κ n+1 σ (0) = σ (g n+1 , F p,n+1 , κ n+1 )
(73)
ϑ (0) = ϑ(κ n+1 ) There is one equation-of-state evaluation at the beginning of the iteration in which the Cauchy stress σ , work-hardening modulus ϑ, and the derivatives ∂σ/∂F p |E,g,κ , ∂σ/ ∂κ|E,g,F p , and ∂ϑ/∂κ, are calculated. Next, evaluate the yield criterion f (m) = f (σ (m) , ϑ (m) ).
(74)
147
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
If m = 0 and f (m) ≤ ², then the state point is interior to the yield surface, and no plastic flow occurs. If f (0) > ² and | f (m) | ≤ ², then F p,n+1 ← F p(m)
(75)
κ n+1 ← κ (m) and stop. Otherwise, calculate 1ζ (m) = ζ (m+1) − ζ (m) using Newton’s method 1ζ (m) = − f (m)
µ
d f (m) dζ
¶−1 (76)
with d f /dζ estimated from d f (m) ≈ dζ
¯ ¶ ·µ ¶µ ˙ p ¶(m) µ ¯ ¶µ ¶(m) ¸ ∂ f (m) ∂σ ¯¯ ∂σ ¯¯ dF d κ˙ + ¯ ¯ p ∂σ ∂F κ,E,g dζ ∂κ F p ,E,g dζ µ ¶(m) µ ¶µ ¶(m) d κ˙ ∂ϑ ∂f . + ∂ϑ ∂κ dζ µ
(77)
Next, calculate revised estimates µ F˜ p = F p(m) +
∂ F˙ p ∂ζ
¶(m)
1ζ (m)
F p(m+1) = (det F˜ p )−1/3 F˜ p µ ¶(m) ∂ κ˙ (m+1) (m) κ =κ + 1ζ (m) ∂ζ ¯ µ ¶ µ ¯ ¡ p(m+1) ¢ ∂σ ¯¯ ∂σ ¯¯ (m+1) (0) p(0) =σ + −F F σ + ¯ p ∂F ∂κ ¯ E,g,κ
ϑ
(m+1)
= ϑ(κ
(m+1)
(78) ¶
¡
E,g,F p
κ (m+1) − κ (0)
¢
),
set m ← m + 1, and retest the stopping criterion. In this procedure we evaluate the equation of state once to determine σ and the thermodynamic derivatives ∂σ/∂F p |g,E,κ and ∂σ/∂κ|g,E,F p . The stress σ (m) , for m > 0, is approximated by first-order Taylor expansion about the initial m = 0 value. The method converges in 1 or 2 iterations, with ² = 10−6 , in each of the test problems involving plasticity described below. The framework described by Eq. (9) calls for rates of plastic deformation h and rates of work hardening K . In the example above, which is rate-independent, we use τ h = 1F p τ K = 1κ,
(79a) (79b)
where τ = 1t/2 in the predictor step of the method (Eqs. (30a) and (30b), and τ = 1t in the corrector (Eq. (41)). A generalization of this approach to rate-dependent plasticity is described in [19].
148
MILLER AND COLELLA
6. DISSIPATION
In certain problems in hydrodynamics it has been found that the higher order Godunov strategy we adapted here will give rise to spurious post-shock oscillations (e.g., [5]). A solution that rectifies this problem is the addition of a small amount of dissipation at strong shocks. This dissipation is added by introducing an additional slope limiter via a “flattening” parameter χ (see Eq. (27)). A variety of flattening strategies have been proposed. Perhaps the simplest variant, employed by Miller and Puckett [15], uses the divergence of the velocity to detect potential shocks, and uses a simple measure of shock strength, the ratio of pressure jump across a cell to the isentropic bulk modulus, |1P|/K S where K S = ∂ P/∂ log ρ|S , to compute a flattening measure. This introduces additional dissipation in regions where the pressure change is large compared to the bulk modulus—where linearization of the equation of state is expected to become error-prone. This strategy may introduce extra dissipation in regions that do not require it, however, as when a shock is spread over a large (>5 or 6) number of grid cells. It is therefore desirable to also include measures of the shock structure to minimize application of this dissipation mechanism. Elaborate strategies for computing χ are described by Colella and Woodward [7]. One of their strategies is to restrict the use of this dissipative mechanism to regions where the detected shock is steep. In our solid mechanics computations we found this strategy to be useful, and in conjunction with a measure of shock strength it provides judicious, adequate additional dissipation. We detect a strong shock by measuring in 1D the divergence of the velocity field, and calculating a normalized jump in stress. We define zi =
k(σ eα )i+1 − (σ eα )i−1 k∞ (det Aαα,i )1/3
(80)
as a measure of shock strength in the neighborhood of cell i in direction eα . The numerator is the maximum of the absolute value of the jump in those stress components that may change in direction eα 1D purely elastic flow, and the denominator is a mean modulus of the acoustic propagation tensor in direction eα . Following Colella and Woodward, we discriminate between steep and broad shocks by the ratio βi =
k(σ eα )i+1 − (σ eα )i−1 k∞ . k(σ eα )i+2 − (σ eα )i−2 k∞
(81)
In the limit βi = 12 , stress is approximately linear across five grid cells, and so a shock discontinuity is not being captured. When βi ≈ 1 the discontinuity is captured in three cells: the shock may be overly steep and postshock oscillations are expected. Accordingly, the minimum value that our flattening parameter χ should have, based upon shock steepness, is µ µ ¶¶ a1 − βi , , (82) χmin i = max 0, min 1, a1 − a0 where a0 and a1 are numerical constants. We use the values a0 = 0.75 and a1 = 0.85 in the computations presented here.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
149
A local shock-strength-sensitive flattening parameter χ˜ i , χmin i ≤ χ˜ i ≤ 1, is thus ¡ ¡ − zi ¢¢ min 1, max zz11 − , χmin i z0
(v · e j )i+1 < (v · e j )i−1
1
otherwise.
( χ˜ i =
(83)
In our example calculations we use the numerical values z 0 = 0.25 and z 1 = 0.75. In 1D we limit the slopes by the minimum over nearest neighbor cells of the local flattening parameter, χi = min(χ˜ i−1 , χ˜ i , χ˜ i+1 ).
(84)
In higher dimensions, we employ the same 1D local flattening parameters—measured separately in each direction. All slopes (∂q/∂ x, ∂q/∂ y, and ∂q/∂z) are limited by the same cell-valued flattening parameter, which is given by the minimum of the directional local measures. In 2D, χi j = min(χ˜ x,i−1, j , χ˜ x,i, j , χ˜ x,i+1, j , χ˜ y,i, j−1 , χ˜ y,i, j , χ˜ y,i, j+1 ),
(85)
and in 3D, χi jk = min(χ˜ x,i−1, j,k , χ˜ x,i, j,k , χ˜ x,i+1, j,k , χ˜ y,i, j−1,k , χ˜ y,i, j,k , χ˜ y,i, j+1,k , χ˜ z,i, j,k−1 , χ˜ z,i, j,k , χ˜ z,i, j,k+1 ).
(86)
7. ACCURACY
The term (v × ∇ × g T )T was introduced to the evolution equations of the inverse deformation gradient g to make the system of equations stable and well-posed when the gauge constraint ∇ × g T = 0 fails to be satisfied. Although the partial differential equations show that when satisfied initially, it will be satisfied for all times, numerical errors cause the constraint to be violated to some degree. We propose a modification of (9) to control inaccuracy that may arise from violation of the gauge constraint. The conservation law (11) indicates that G will be created by numerical errors as dipoles. Thus, a numerical strategy that will control this truncation error is to diffuse G, ∂G + ∇ · (vG − Gv) = D(∇ 2 G) ∂t
(87)
∂geα ∂ (gv) = (v × G)T eα − D(∇ × G)T eα . + ∂t ∂ xα
(88)
or, equivalently,
g is also related to the density via ρ = ρ0 det g,
(89)
150
MILLER AND COLELLA
where ρ0 is the mass density in the reference state F = g = I . Multiplying the g equations by ρ0 det(g)g −T , and summing over the nine components of g gives a conservation law for ρˆ ≡ ρ0 det(g): ∂ ρˆ + ∇ · (ρv) ˆ = 0. ∂t
(90)
Thus the continuity equation is embodied in the g equations as well. However, because of discretization errors the equivalence of ρˆ and the mass density ρ cannot be assured. To make the method strictly conservative, we keep ρ as a redundant variable, and we invoke a relaxation mechanism on g to enforce the condition ρˆ = ρ. This relaxation alone (not including the diffusion modification) is accomplished by writing ∂ ∂geα (gv) = (v × G)T eα + η + ∂t ∂ xα
µ
¶ ρ − 1 geα . ρˆ
(91)
The “continuity” equation for ρˆ is then Dgαβ D ρˆ ˆ = ρF ˆ βα = −ρ∇ ˆ · v + ρF ˆ βα (v × G)βα + 3η (ρ − ρ) Dt Dt ∂ ρˆ when G = 0. ˆ + ∇ · (ρv) ˆ = 3η (ρ − ρ) ∂t
(92) (93)
Including the diffusion and relaxation terms, the system of equations we will solve is ρvα ρ ρvv − σ e α α ρv ρ Evα − vβ σβα ρE gvδxα gex gvδ yα ∂ ge y ∂ gez + ∂ xα ∂t gvδ zα ρF p ex ρF p ex vα ρF p e y ρF p e v y α p ρF ez ρF p ez vα ρκ ρκvα 0 ρf ρ(8 + v · f ) T −D(∇ × G)T ex e (v × G) 1 (v × G)T e −D(∇ × G)T e 2 y = + (v × G)T e3 T −D(∇ × G) ez ρhe x ρhe y ρhez ρK
0 0 0 +η +η +η 0 0 0 0
¡ ¡ ¡
ρ ρ0 det g ρ ρ0 det g ρ ρ0 det g
¢ − 1 gex ¢ − 1 ge y . ¢ − 1 gez
(94)
151
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
Our discretization of the diffusion and relaxation terms takes the form ¡ ρ˜ n+1 ¢ n+1 − 1 g˜ ex −D(∇ × G)T,n ex + η ρ0 det n+1 g˜ n+1 n+1 ˜ x ge gex ¡ ¢ ρ˜ n+1 T,n n+1 ˜ −D(∇ × G) e + η − 1 g e ge y = ge y y ˜ y + 1t ρ0 det g˜ n+1 ¡ ρ˜ n+1 ¢ n+1 ˜ z i jk ge gez i jk −D(∇ × G)T,n ez + η ρ0 det − 1 g˜ ez n+1 g˜
,
(95)
i jk
where g˜ n+1 denotes g after flux differencing and evaluation of source terms in (10) (cf. (40)). The second derivatives of g appearing in the G diffusion term, £ 2 £ ∂ 2 gx x £ ∂ 2 gx y 2 2 ∂2g ∂ gx z − ∂∂zg2x x − ∂ x 2x y − ∂∂ yg2x z ∂ x∂z ∂ y∂ x ∂z∂ y ¤ ¤ 2 ∂ 2 gx y ¤ ∂ 2 gx z ∂ 2 gx x ∂ 2 gx z + ∂ gx y − ∂ 2 g2x x + ∂ y∂z − ∂z 2 + ∂z∂ x − ∂ x 2 ∂ x∂ y ∂y £ £ £ 2 ∂ 2 g yz ∂ 2 g yx ∂ 2 g yy ∂ 2 g yy ∂ 2 g yz ∂ x∂z − ∂∂zg2yx − − ∂ x∂ y ∂x2 ∂ y∂z ∂ y2 T (∇ × G) = , (96) ∂ 2 g yz ∂ 2 g yy ¤ ∂ 2 g yx ∂ 2 g yz ¤ + ∂ 2 g yy − ∂ 2 g yx ¤ − − + + 2 2 2 ∂ x∂ y ∂y ∂ y∂z ∂z ∂ x∂z ∂x £ ∂ 2 g £ £ 2 2 ∂ gzy ∂ gzy ∂ 2 gzx ∂ 2 gzx ∂ 2 gzz zz − − − ∂ x∂z ∂z 2 ∂ x∂ y ∂x2 ∂ y∂z ∂ y2 ¤ ¤ ¤ 2 2 2 2 2 2 ∂ g ∂ g gzz gzx + ∂ x∂zyy − ∂∂ yg2zx − ∂z 2zy − ∂∂ xg2zz + ∂∂ y∂z + ∂∂ x∂z are computed using time-n cell-centered values of g, with a standard three-point stencil for homogeneous second derivatives, e.g., ¢ ¢¶ ¡ µ ¡ n µ 2 ¶ n 2 gi+1, jk − ginjk 2 ginjk − gi−1, 1 ∂ g jk , (97) = − ∂ x 2 i jk 1xi 1xi+1 + 1xi 1xi + 1xi−1 and heterogeneous derivatives are computed with a four-point stencil, e.g., µ 2 ¶ n n n n gi+1, ∂ g j+1,k − gi+1, j−1,k − gi−1, j+1,k + gi−1, j−1,k . =4 ∂ x∂ y i jk (1xi−1 + 21xi + 1xi+1 )(1y j−1 + 21y j + 1y j+1 )
(98)
A von Neumann stability analysis of the diffusion update in (94), considered independently of other source terms or the basic solid mechanics equations, gives a bound on the diffusion coefficient: 2 h in 1D, 21t D≤ (99) 2 h in 2D or 3D. 41t This suggests an approximate overall Courant–Friedrichs–Lewy stability criterion of ( 1t (|v| + cmax ) CFL = CFL < 1,
1x
+
2D1t (1x)2
in 1D,
1t (|v| + cmax ) 1x
+
4D1t (1x)2
in 2D or 3D,
(100) (101)
where here it is assumed that 1x = 1y = 1z, a constant. The more rigorous CFL condition 4D1t CFL = max( 1t (|v|1x+ cmax ) , (1x) 2 ) (in 2D or 3D) would hold if the mechanics equations and G diffusion steps were performed sequentially.
152
MILLER AND COLELLA
The optimal damping conditions for the diffusion of G are obtained by choosing the empirical diffusion constant D to satisfy D=
(1x)2 , 4d1t
(102)
where d = 1, 2, 3 is the dimensionality of the problem. Similarly, optimal relaxation is obtained by choosing the empirical relaxation parameter η to satisfy η=
1 61t
(103)
for all dimensions. According to the approximate CFL condition (101), the optimal value of D will contribute 1/2 to the CFL value in 1D and 2D, and 1/3 in 3D, limiting the overall step size 1t by factors of 1/2 and 2/3 (respectively) relative to the D = 0 value. Thus, some of the examples presented below use smaller values of D than indicated by (102). In some cases, however, we find that values of CFL > 1 provide stable solutions [consistent with (101) being only an approximation]. Assumptions underlying our plastic yield model require that det F p be constant. The differential equations describing our plastic flow model F˙ p preserves det F p , but again numerical errors will lead to some violation of this constraint. To remedy this problem we renormalize the plastic deformation tensor at the end of each time step, F p ← (det F p )−1/3 F p .
(104)
8. EXAMPLES
8.1. Convergence: Elasticity To demonstrate the convergence properties of the algorithm we model in 1D the smooth flow resulting from an initial Gaussian-shaped disturbance. For these computations we use a hyperthermoelastic model of the Mooney–Rivlin variety, ρ0 E(C e , S) =
√ λ(S) µ(S) µ(S) (log det C e )2 + tr C e − log det C e 2 2 2 µ ¶ 1 −ϑ1 κ ρ0 ϑ0 κ+ e , + ρ ϑ1
(105)
where S is the entropy. Entropy dependence is introduced by supposing λ(S) = λ0 + λS f (S)
(106a)
µ(S) = µ0 + µS f (S)
(106b)
where f (S) is an unspecified function of the entropy. From this equation we evaluate σ (E, g, F p , κ), and other derivatives including the acoustic propagation tensors, by first solving this equation of state for f (S) and then differentiating E with respect to the elements of C e while holding f (S) constant. We use values ρ0 = 1, µ0 = λ0 = 0.6, and µS = λS =
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
153
p
0.01, with initial values g0 = 1.1I , F0 = I , κ = 0, and v = 0. The initial disturbance is generated by distributing internal energy from E 0 to E 00 , ¡ ¢ p E 0 = E g0 , F0 , f (S) = 0 00
E = 10E0 .
(107a) (107b)
These limiting values are used to construct a Gaussian initial profile, via Ei = E 0 ωi + (1 − ωi )E 00
(108)
¸ · 1 r2 ωi = √ exp − i 2 , 2a a 2π
(109)
with
where a 2 = 100 is the variance of the distribution, and where ri is the coordinate of the center of cell i, in the domain [0, 40]. Boundary conditions are reflecting at r = 0 and r = 40. We pick the time step 1t to satisfy the Courant–Friedrichs–Lewy constraint (101), with CFL = 0.8. This problem was chosen to give a nontrivial shockless flow, with initial conditions that strictly obey G = 0. Plasticity is not incorporated into this test problem, and no flattening is required. Figure 1 shows the initial and final conditions of this test problem in Cartesian geometry. At this scale, the difference between results at 40, 80, and 160 Cartesian points is not resolvable. A comparison of results using 40, 80, and 160 grid points is used to estimate the L 1 , L 2 , and L ∞ (max) norm rates of convergence using the volume-weighted variables (Table I). In Cartesian geometry the method exhibits approximately third-order convergence: as the number of grid cells is doubled, the error diminishes by a factor of 23 . Slightly lower rates of convergence are seen in cylindrical and spherical geometries, but in all cases the order exceeds 2.
FIG. 1. Initial conditions and computed results using 160 cells in Cartesian geometry for (a) density, ρ; and (b) stress, σx x .
154
MILLER AND COLELLA
TABLE I Convergence Test: Pure Elasticity Geometry
Field
L1
L2
L∞
Cartesian
ρ vx σx x σ yy , σzz
3.33 3.02 3.30 2.95
3.24 2.97 3.47 2.91
3.06 2.89 3.80 2.69
Cylindrical
ρ vr σrr σθθ σzz
2.51 2.84 2.82 3.12 3.22
2.68 2.70 2.78 3.24 3.17
2.80 2.53 2.68 3.30 2.97
Spherical
ρ vr σrr σθθ , σφφ
2.42 2.77 2.85 3.31
2.53 2.64 2.78 3.34
2.66 2.55 2.79 3.37
8.2. Convergence: Plasticity To assess the rate of convergence in a plasticity-dominated flow we pose a model problem similar to the purely elastic problem presented in (Fig. 2). A Gaussian distribution with width
FIG. 2. Initial conditions and computed results using 160 cells in Cartesian geometry for (a) density, ρ; (b) stress, σx x ; (c) plastic deformation Fxpx ; and (d) work hardening parameter κ.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
155
TABLE II Convergence Test: Elastic–Plastic Flow Geometry
Field
Cartesian
ρ vx σx x σ yy σzz Fxpx F yyp Fzzp κ
L1
L2
L∞
2.14 2.07 1.99 1.62 1.89 2.01 2.18 2.36 2.05
2.03 2.08 1.93 1.55 1.76 1.80 2.08 2.32 1.98
1.84 2.36 2.12 1.21 1.50 1.60 2.10 2.71 1.91
5 is used to vary g yy and gzz as functions of coordinate x according to gx x,i = 1.1 g yy,i = (1 + 9ωi )1.1
(110)
gzz,i = 1.1/(1 + 9ωi ) with homogeneous initial density, internal energy, and zero velocity. We use the equation of state (105) with yield model (62) and flow rates (63, 64). The equation-of-state parameters are as used in the purely elastic convergence test, and the plastic constitutive parameters are σY = 0.1, ϑ0 = 0.1, and ϑ1 = 10.0. The flow field in this problem is C 0 , which lowers the overall order of convergence. Density converges at greater than second order (Table II), but the tangential stress components converge only at first order. 8.3. Blake’s Problem Blake [3] presented an analytical solution to the problem of an unbounded solid medium characterized by an isotropic linear elastic equation of state, ρ0 E =
1 1 λ[tr(C e − I )]2 + µ tr(C eT C e − 2C e − I ), 8 4
(111)
loaded by a prescribed pressure boundary condition on the interior of a spherical cavity of initial radius a. We present a numerical solution to this problem in 1D spherical coordinates (see Appendix), with slight modification of the code to accommodate the moving boundary with prescribed flux (Neumann) boundary conditions. This problem is selected to verify the behavior of the elastic algorithm in the weak shock limit. The cavity wall represents a material interface across which the mass flux will be zero. Accordingly, the flux at this boundary is given by F(U B ) − sU B , where U B is the vector of conserved quantities at the boundary, s is the velocity of the boundary, and F(U B ) is the radial flux vector evaluated at the boundary. Blake’s solution provides u(r, t), the displacement of a mass element in the radial direction. In spherical coordinates, this gives
156
MILLER AND COLELLA
rise to an inverse deformation tensor (1 + ∂u/∂r )−1 0 g(r, t) = 0 (1 + u/r )−1 0
0
0
.
0
(112)
−1
(1 + u/r )
The velocity of the material interface is s = ∂u/∂t|r =a . Cell 1, whose left boundary is r = a at t = 0, and whose right boundary is fixed at a + 1r , has a volume which varies with time. Applying Gauss’s divergence theorem to this cell gives ¡ n+1/2 n+1/2 V1n+1 U1n+1 = V n U1n + 1t A1/2 F1/2 − A3/2 F3/2 ¡ n+1/2 n+1/2 ¢¢ n+1/2 + A¯ 1 H1/2 − H3/2 , (113) + 1t V¯ 1 G 1 where F denotes the radial flux component that enters as (1/r 2 )∂(r 2 F)/∂r , H denotes the radial flux component that enters as ∂ H/∂r (see Appendix), A¯ 1 is the average area (r 2 ) n+1/2 over r in [a − u(a, t), a + 1r ], V¯ 1 is the time-averaged cell volume, and G 1 is the cell-centered vector of (geometric) source terms, which we time-center with a predictor– corrector strategy. In general (see Wilkins’s problem below), an algebraic solution of this discretization is unstable. In the particular case of our discretization of Blake’s problem, however, |u(a, t)| ¿ 1r and so V1 does not vary appreciably with time and in particular is of order a 2 1r . Our solution to Blake’s problem therefore uses (113) as written. It is also necessary to modify the algorithm to account for the absence of cell values at i − 1 and i − 2. The gradient ∂q/∂r at i = 1 is obtained by first order forward finite difference with a van Leer limiter. The flattening parameter χ operates on a stencil that requires cell values at 0 and −1. However, for this weak problem additional flattening is never required, so the algorithm is modified by omission of the flattening computation (χ = 1). Following Trangenstein and Colella [25] we use parameters a = 0.1 m, ρ0 = 3000 kg/m3 , λ = 2.36 × 1010 Pa, and µ = 2.78 × 1010 Pa. The pressure inside the spherical cavity is 106 Pa, and the solution is plotted at time 1.6 × 104 s. We compare in Figs. 3–6 our computed results for radial stress, σrr = (λ + 2µ)(∂u/∂r ) + 2λ(u/r ),
(114)
σθ θ = σφφ = λ(∂u/∂r ) + 2(λ + µ)(u/r ),
(115)
¶ µ 1 2 P = − σii = − λ + µ [(∂u/∂r ) + 2(u/r )] , 3 3
(116)
hoop stress
pressure
and radial velocity vr against Blake’s analytical results. These results verify the method in the case of weak (linear) waves. The leading shock is captured in approximately five grid cells. A single stress undershoot precedes the shock, and a corresponding overshoot follows it, but the wave speed and amplitude are correctly modeled.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
157
FIG. 3. σrr calculated from Blake’s analytical solution at time 1.6 × 104 , compared with values computed using 200 grid cells.
FIG. 4. Hoop stress, σθθ = σφφ .
FIG. 5. Pressure, (σrr + σθ θ + σφφ )/3.
158
MILLER AND COLELLA
FIG. 6. Material velocity vr for Blake’s problem.
8.4. Wilkins’s Problem Wilkins’s flying plate problem [27] involves a 5-mm-thick aluminum plate impacting an initially stationary aluminum halfspace. The rear (left) surface of the flying plate is a free surface (vacuum). Initially, left- and right-traveling shocks propagate outward from the point of contact of the plate with the halfspace. When the left-traveling shock reaches the free surface, a right-traveling rarefaction is created, which ultimately overtakes the right-traveling shock. This problem incorporates plasticity. To model this problem, we modify our 1D algorithm to allow for the moving free-surface boundary. This is an example of volume-of-fluid front reconstruction applied to multi-fluid modeling, and details will be described in a future correspondence. Briefly, we modify the approach adopted for Blake’s problem using the flux redistribution ideas of Chern and Colella [4]. Application of this approach to stationary incompressible boundaries is described in [16], and to reaction front tracking in [1, 18]. Our implementation is similar, but the free-surface boundary moves at a velocity determined by the solid-vacuum Riemann problem. This problem is solved as described above for the solid–solid case but uses only the 3 × 3 stress component of the eigenvectors. This interface velocity, and the surrounding material velocities, are used with a volume-pushing algorithm (after [2]) to update the fractional occupancy of the interface cells. We construct a hyperelastic model of aluminum in close correspondence to Wilkins’s (rate model) description, with ¶ µ µZ µ ¶2/3 ¶ µ0 P(ρ 0 ) 0 ρ0 e + tr C , (117) dρ − 3 E(g, F p ) = 02 ρ 2ρ0 ρ where P(ρ) is the hydrostatic pressure (in GPa) P(ρ) = 72(ρ/ρ0 − 1) + 172(ρ/ρ0 − 1)2 + 40(ρ/ρ0 − 1)3 ,
(118)
with ρ0 = 2.7 kg/m3 . The shear modulus is µ0 = 24.8 GPa. The problem is perfectly plastic (no work hardening), and uses the von Mises yield surface function r 2 (119) σY f (σ ) = kdev σ k − 3 with constant flow stress σY = 0.2976 GPa.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
159
FIG. 7. Longitudinal stress σx x for Wilkins’s problem with impact velocity 0.8 km/s. Time in µs.
Computations with impact velocities of 0.8 km/s (Figs. 7 and 8) and 2.0 km/s (Figs. 9 and 10) were obtained with CFL = 0.80 and 500 Cartesian grid points. At 0.8 km/s, a plastic shock trails a leading elastic shock precursor. When the left-facing shocks reach the free surface, right-traveling elastic and trailing plastic rarefaction waves begin to overtake the initial right-facing shocks. The shock stress at 2.0 km/s is above the elastic limit, so only plastic shocks are formed. On rarefaction from the left free surface, a leading rightfacing elastic rarefaction is formed, followed by the plastic wave. These results are in good quantitative agreement with those of Wilkins. 8.5. A Test in 2D This test problem compares a 1D cylindrical coordinate computation against a 2D Cartesian result, for a problem with cylindrical symmetry. We use the modified Mooney– Rivlin model presented in Eq. (70) with initial conditions ρ0 = 1, g = 1.1I , F p = I , and κ = 0. The plasticity parameters are σY = 0.1, ϑ0 = 0.1, and ϑ1 = 10. All boundary conditions are reflecting. The material is initially at rest, except for a cylindrical shell r ∈ [5, 15], which moves toward the axis with a velocity of −1. This generates a diverging rarefaction, and a convergent shock, which reflects off the axis of symmetry.
160
MILLER AND COLELLA
FIG. 8. Mass density ρ for Wilkins’s problem with impact velocity 0.8 km/s. Time in µs.
In Figs. 11–15 we compare results from a 1D cylindrical calculation (500 cells, CFL = 0.8), and an equivalent 2D Cartesian calculation using 250 × 250 cells, also at CFL = 0.8. The 2D results are presented as 1D scatter plots in order to demonstrate the accurate preservation of cylindrical symmetry obtained with the spatially unsplit 2D method. The high-resolution 1D results and lower resolution 2D results are in good agreement, although there is some discrepancy in κ and σrr near the axis. Using this same 2D test we demonstrate the errors associated with the gauge constraints ρ − ρ0 det(g) = 0
(120)
G = ∇ × g T = 0.
(121)
and
These conditions are enforced in the computation by way of a relaxation term to satisfy (120) and a diffusion-like term to satisfy (121). In Fig. 16 we plot the left-hand side of
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
161
FIG. 9. Longitudinal stress σx x for Wilkins’s problem with impact velocity 2.0 km/s. Time in µs.
(120) comparing results from the computation presented in Figs. 11–15, and results from a similar computation in which, however, neither a relaxation nor a diffusion correction was applied. In Fig. 17 we plot the L 2 norm of the tensor ∇ × g T , comparing results from the computation with relaxation and diffusion to results from a computation using neither correction. These figures demonstrate over an order of magnitude reduction in density error is achieved by the relaxation mechanism. Approximately a factor of 2 reduction of ||G||2 is achieved by the diffusion mechanism. 8.6. A Test in 3D This test problem compares a 1D spherical coordinate computation against a 3D Cartesian result, for a problem with spherical symmetry. The equation of state is identical to the 2D test above, and the initial conditions are similar: a spherical shell r ∈ [5, 15] is given an initial velocity of −1. This computation, with 100 × 100 × 100 cells at CFL = 0.8, is underresolved. Nevertheless, there is good agreement between the 3D Cartesian results and the 1D spherical calculation and excellent preservation of spherical symmetry (Figs. 18–20).
162
MILLER AND COLELLA
FIG. 10. Mass density ρ for Wilkins’s problem with impact velocity 2.0 km/s. Time in µs.
FIG. 11. Density at time 4.
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
FIG. 12. Radial velocity at time 4.
FIG. 13. σrr at time 4.
FIG. 14. F p rr component of plastic deformation tensor.
163
164
MILLER AND COLELLA
FIG. 15. Work-hardening parameter κ.
FIG. 16. Error in density constraint, ρ − ρ0 det(g).
FIG. 17. Error in curl constraint,
p
Gαβ Gαβ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
FIG. 18. Density at time 4.
FIG. 19. Radial velocity at time 4.
FIG. 20. σrr at time 4.
165
166
MILLER AND COLELLA
9. CONCLUSIONS
We have presented a new method for the solution to equations of solid mechanics in one, two, and three spatial dimensions on Eulerian grids. Our method addresses the problem of gauge constraints (∇ × g T = 0) by adopting a nonconservation approach first proposed by [8] for the equations of magnetohydrodynamics. We write the partial differential equations of solid mechanics in such a way that the constraint, if applicable initially, holds true for all time. The constraint is violated by the truncation error of the method, and reinforced with an explicit diffusion term which annihilates the dipolar field of ∇ × g T . Another constraint of the system, a correspondence between density variation and the deformation field (ρ = ρ0 det g) is also satisfied for all times by the PDEs, if satisfied in the initial conditions. Truncation errors of the method are compensated with an explicit relaxation term. The method presented here does not incorporate artificial viscosity, but its solutions are sensitive to six adjustable parameters: D and η control accuracy of the gauge constraints, and a0 , a1 , z 0 , and z 1 in Eqs. (82) and (83) govern the introduction of dissipation near strong shocks to prevent overshoot and ringing by locally reducing the high-order Godunov method to first order. Our strategy for damping modes violating the curl gauge constraint, g T : = g T − λ∇ × ∇ × g T , = g T + λ(∇ 2 g T − ∇(∇ · g T )),
(122)
= g + λ(∇ g − ∇ Q(g )), T
2 T
2
T
(λ = 1tD), uses a single central difference operator acting on cell-centered variables. Here, Q(g T ) = ∇ −2 ∇(∇ · g T ) is the projection onto the curl-free part of g T . Defining P(x) = 1 − Q(x) as the projection onto the divergence-free part of x, and noting P Q = Q P = 0, we have P(g T ) := P(g T ) + λ∇ 2 P(g T );
(123)
thus we are diffusing the divergence-free part of g T without modifying the curl-free part. A similar scheme may be used to modify a vector field B subject to a divergence-free constraint, B := B + λ∇(∇ · B), Q(B) := Q(B) + λ∇ 2 Q(B)
(124)
with a single matrix-valued central difference operator for the projection ∇(∇ · B). This will directly target odd–even and checkerboard short-wavelength modes of ∇ · B by diffusing the curl-free part of B. The application of this extension to magnetohydrodynamics, where B is the magnetic field subject to gauge constraint divB = 0, is currently being investigated (R. Crockett, personal communication). APPENDIX: CYLINDRICAL AND SPHERICAL COORDINATES
The equations of solid mechanics in cylindrical and spherical coordinates (like those of gas dynamics) differ from the Cartesian equations by the existence of both spatial and
167
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
volumetric spatial derivatives, and by the introduction of “geometric source terms.” The coordinate transformation is accomplished by rotating the Cartesian basis vectors into the curved coordinate frame via the rotation matrices
Rcyl
0 0 1
cos θ −sin θ = sin θ cos θ 0 0
sin θ cos φ
cos θ cos φ
cos θ
cos θ sin φ −sin θ
Rsph = sin θ sin φ
−sin φ
(A.1)
cos φ , 0
(A.2)
where we adopt the standard curved coordinate notation x = r cos θ y = r sin θ
(A.3)
z=z in cylindrical coordinates and x = r sin θ cos φ y = r sin θ sin φ
(A.4)
z = r cos θ in spherical coordinates. R is the matrix of inner products of unit vectors in the curved coordinate system, eα0 , and the Cartesian system eβ ; Rαβ = eα0 · eβ . These rotation matrices transform the Cartesian tensors F, g, and σ transform as σcyl = R T σCart R, etc. and transform the velocity vector v as vcyl = R T vCart . In cylindrical coordinates, the system of transformed equations may be written (cf. Eq. 9) as
0 ! à ! ρvθ σrr ρ ρvv − σ0 rθ ρvvθ − σ eθ r ρv − 0 σr z ρE ρ Evθ − v T σ eθ 0 ge 0 ρ Evr − v T σ er 0 r ∂ 1 ∂ 1 ∂ ∂ gv gv ge 0 θ + + + r r ∂r 0 0 ge ∂t ∂r r ∂θ z 0 p p ρF er 0 ρF er vθ 0 p ρF p e v p 0 ρF eθ ρF e v θ θ r r 0 ρF p e ρF p e v p z θ r ρF e v z θ 0 ρF p ez vr ρκ ρκv ρvà r
ρκvr
0
θ
168
MILLER AND COLELLA
0 ¡ ¢ − ρv 2 + σrr − σθ θ (ρvθ v − σ ) r θ rθ 0 0 0 −g v − g v − g v θr r θθ θ θz z ρvz grr vr + gr θ vθ + gr z vz ρvvz − σ ez 0 ρ Ev − v T σ e z z 0 0 p p −Fr θ − Fθr ∂ 1 0 + = − p gv ∂z r ρvθ Frrp − Fθ θ p ρF er vz p −F zθ p ρF eθ vz ρF p e v p p z z Frr − Fθ θ ρκvz p p F + F ρv θ rθ θr p Fzr p −Fθ z Frpz ρv θ 0
0 0 ρf ρ(8 + v · f ) rr − vz ∂g∂zrr − vθ ∂g + r ∂θ
v ∂gr θ + v ∂gr z vθ (gr θ + gθr ) θ ∂r z ∂r r ∂gθθ ∂gθ z ∂gθr ∂gθr vθ (gθθ − grr ) vθ ∂r + vz ∂r − vθ r ∂θ − vz ∂z + r zr vθ ∂g∂rzθ + vz ∂g∂rzz − vθ ∂g − vz ∂g∂zzr + vθ rgzθ r ∂θ v ∂grr + v ∂gr z − v ∂gr θ − v ∂gr θ − vr (gr θ + gθr ) + vz gθ z r r ∂θ z r ∂θ r ∂r z ∂z r v ∂gθr + v ∂gθ z − v ∂gθθ − v ∂gθ θ + vr (grr − gθ θ ) + vz gr z r r ∂θ z r ∂θ r ∂r z ∂z r . + vr gzθ ∂gzr ∂gzz ∂gzθ ∂gzθ vr r ∂θ + vz r ∂θ − vr ∂r − vz z − r ∂grr ∂gr θ ∂gr z ∂gr z vθ gθ z vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ + r vθ gr z ∂gθr ∂gθθ ∂gθ z ∂gθ z vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ − r ∂gzr ∂gzθ ∂gzz ∂gzz vr ∂z + vθ ∂z − vr ∂r − vθ r ∂θ ρher ρheθ ρhez
(A.5)
ρK
This equation does not include the det g relaxation term, whose representation is unaffected by the change in variables, nor does it include the G diffusion correction, which will be described separately below. There is some latitude in the partitioning of terms between the LHS and the RHS geometric source vector. This is particularly evident in the stress terms appearing in the momentum equations. The choice of representations described here was chosen in order that some cancellation between σrr and σθ θ occur in the r -momentum source term.
169
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
The linearized equations of solid mechanics [cf. (28)], used in the construction of L and R edge states, also has a geometric source vector. Expressed in terms of the primitive variables q, but omitting stress components which are described later, we have −ρvr vθ2 + σrr /ρ − σθθ /ρ −v v + 2σ /ρ r θ rθ σr z /ρ (vr σθθ − vθ σr θ )/ρ Ã ! 0 0 0 Ã v g +v g +v g ! r θr θ θθ z θz −vr grr − vθ gr θ − vz gr z 0 Ã ! 0 1 . s= (A.6) 0 r 0 p Fr θ + Fθrp p p ρvθ Fθ θ − Frr p Fzθ p Fθ θ − Frrp ρvθ −Frpθ − Fθrp −Fzrp p Fθ z ρvθ −Frpz 0 0
The stress evolution equations, used in the predictor steps of the method, are (plastic source terms omitted) σ er ∂ σ eθ ∂t σ ez
+ vr
Arr = Aθr Azr
σ er ∂ σ eθ ∂r σ ez
vr · ∂ vθ ∂r vz
+ vθ
σ er ∂ σ eθ r ∂θ σ ez
Ar θ + Aθ θ Azθ
+ vz
vr · ∂ vθ r ∂θ vz
2vθ σr θ + vr (Ar θ )r θ − vθ (Ar θ )rr
σ er ∂ σ eθ ∂z σ ez
Ar z
where the tensors A are defined by (14).
vr + A θ z · ∂ vθ ∂z vz Azz
− vθ (σrr − σθθ ) + vr (Ar θ )θ θ − vθ (Ar θ )θr vθ σθ z + vr (Ar θ )r θ − vθ (Ar θ )zr −v (σ − σ ) + v (A ) − v (A ) θ rr θθ r θθ rθ θ θ θ rr 1 , + −2vθ σr θ + vr (Aθ θ )θ θ − vθ (Aθ θ )θr r −vθ σr z + vr (Aθθ )zθ − vθ (Aθθ )zr vθ σθ z + vr (Azθ )r θ − vθ (Azθ )rr −vθ σr z + vr (Azθ )θθ − vθ (Azθ )θr +vr (Azθ )zθ − vθ (Azθ )zr
(A.7)
170
MILLER AND COLELLA
The g relaxation term,−D(∇ ×∇ ×g T )T , transforms in cylindrical coordinates as −D × £ ∂ 2 g rz ∂r ∂z
∂ 2 gr θ r ∂r ∂θ ∂ 2 grr + grrr2 − grθθ2 r 2 ∂θ 2 ∂gr θ θr θθ + 2∂g − ∂g r 2 ∂θ r 2 ∂θ r ∂r
−
∂ 2 grr ∂z 2
+
£ ∂ 2 grr
2 2 rz − ∂∂rg2r θ + r∂∂θg∂z r ∂r ∂θ 2 rr − ∂∂zg2r θ − r∂g2 ∂θ + grr2θ + grθr2 ¤ rθ θr θz − ∂g − ∂g − ∂g r ∂r r ∂r r ∂z
£ ∂ 2 gr θ
2 gr z ∂ 2 grr − r∂2 ∂θ 2 + ∂r ∂z r ∂θ ∂z 2 θz − ∂∂rg2r z + grr2z + 2∂g r 2 ∂θ ¤ rr rz θθ − ∂g − ∂g + ∂g r ∂z r ∂r r ∂z
− ¤ + £ ∂ 2 g £ £ 2 2 2 2 2 2 2 2 ∂ gθr ∂ gθθ ∂ gθr ∂ gθ θ ∂ gθ z ∂ gθ θ ∂ gθ z ∂ gθr θz − + − + − + ∂r ∂z ∂z 2 r ∂r ∂θ r ∂r ∂θ ∂r 2 r ∂θ ∂z r ∂θ ∂z r 2 ∂θ 2 ∂r ∂z 2 2 2 ∂gθθ ∂ gθ θ ∂gθr ∂grr ∂ gθ z 2∂gr z ∂gr θ rr − ∂2 gθr2 − 2∂g + r 2 ∂θ − ∂z 2 − r 2 ∂θ + r ∂r − ∂r 2 − r 2 ∂θ + r ∂z r ∂θ r 2 ∂θ ¤ ¤ ¤ ∂g g g ∂gθ θ grr gθθ ∂gθr ∂gθ z gθ z rz − − + − + + ∂g + + r ∂rr θ + rr2θ + rθr2 2 2 2 r ∂z r ∂r r r r ∂z r ∂r r £ ∂ 2 g £ £ ∂ 2 gzr ∂ 2 gzθ ∂ 2 gzr ∂ 2 gzθ ∂ 2 gzz ∂ 2 gzθ ∂ 2 gzz ∂ 2 gzr zz − + − + − + ∂r ∂z ∂z 2 r ∂r ∂θ r ∂r ∂θ ∂r 2 r ∂θ ∂z r ∂θ ∂z r 2 ∂θ 2 ∂r ∂z ¤ ¤ ¤ 2 2 2 ∂ gzr ∂gzθ ∂ gzθ ∂gzr ∂gzθ gzθ ∂ gzz ∂gzr ∂gzz − ∂z 2 − r 2 ∂θ − r ∂r + r 2 − ∂r 2 + r ∂z − r ∂r − r 2 ∂θ 2 + r 2 ∂θ (A.8)
The transformed system of equations in spherical coordinates may be written as ρ ρvr ρv 0 ρvvr −σ er ρE ρ Evr − v T σ er ger 0 0 ∂ geθ 1 ∂ 2 ∂ gv 0 + r + ∂r 0 ∂t geφ r 2 ∂r 0 ρF p er ρF p er vr 0 ρF p e ρF p e v 0 θ θ r p p 0 ρF eφ ρF eφ vr ρκ ρκvr
ρvθ ρvvθ T ρ Evθ − v σ eθ 0 1 ∂ 0 ∂ 1 sin θ + + r ∂θ 0 r sin θ ∂θ ρF p er vθ ρF p e v θ θ ρF p eφ vθ ρκvθ
ρvφ ρvv − σ e 0 φ φ T −σ eθ ρ Ev − v σ eφ φ 0 0 0 0 ∂ + 1 gv r sin θ ∂φ gv ρF p er vφ 0 ρF p e v 0 θ φ p 0 ρF eφ vφ ρκvφ
171
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
0 −ρ vθ2 + vφ − 2σrr + σθ θ + σφφ − cot θ σr θ ¢ ¡ ρ vr vθ − cot θvφ2 − 3σr θ + cot θ (σφφ − σθ θ )
¡
¢ 2
ρ(vr vφ + vθ vφ cot θ ) − 3σr φ − 2 cot θ σθ φ 0 0 −gθr vr − gθ θ vθ − gθ φ vφ grr vr + gr θ vθ + gr φ vφ 0 −gφr vr − gφθ vθ − gφφ vφ −cot θ(gφr vr + gφθ vθ + gφφ vφ ) (grr vr + gr θ vθ + gr φ vφ ) + cot θ (gθr vr + gθ θ vθ + gθ φ vφ ) 1 =− ¡ p ¡ p p¢ p ¢ r −vθ Fr θ + Fθr − vφ Fr φ + Fφr ¡ p p¢ p p ρ vθ Frr − Fθ θ − vφ Fθ φ − cot θ vφ Fφr ¡ ¢ p p p p −vθ Fφθ + vφ Frr − Fφφ + cot θ vφ Fθr ¡ p p¢ p p vθ Frr − Fθ θ − vφ Fφθ − cot θ vφ Fr φ ¡ ¢ ¡ ¢ p p p p ρ vθ Fr θ + Fθr − cot θ vφ Fφθ + Fθ φ ¡ ¢ p p p p F + v F + cot θ v − F v F θ φ φ φr rθ θθ φφ ¡ p ¢ p p p −vθ Fθ φ + vφ Frr − Fφφ + cot θ vφ Fr θ ¡ ¢ p ρ vθ Frpφ + vφ Fθrp + cot θ vφ Fθpθ − Fφφ ¡ p ¡ p p ¢ p ¢ vφ Fφr + Fr φ + cot θ vφ Fφθ + Fθ φ 0 0 ρf ρ(8 + v · f ) ¡ ∂gr θ ∂grr ¢ ¡ ∂grr ¢ ∂g v (g + g ) + v (g + g ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂rr φ + θ r θ θr r φ r φ φr ¡ ∂gθr ¡ ∂gθθ ∂gθr ¢ ∂gθφ ¢ vθ (gθθ − grr ) + vφ (gθ φ + cot θgφr ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂r + r ¡ ¢ ¡ ¢ ∂gφθ ∂gφr ∂gφr ∂gφφ vθ gφθ + vφ (gφφ − grr − cot θgθr ) vθ ∂r − r ∂θ − vφ r sin θ ∂φ − ∂r + r ¡ ∂g ¢ ¡ ¢ −vr (gr θ + gθr ) + vφ (gφθ − gθ φ + cot θgr φ ) ∂gr θ ∂gr θ ∂grr rφ vφ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r ¡ ∂g ¢ ¡ ∂gθ θ ∂gθr ¢ vr (grr − gθ θ ) + vφ (gr φ + cot θ (gθφ + gφθ )) ∂gθθ θφ v φ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r + . ¡ ∂g ¢ ¡ ¢ ∂gφθ ∂gφθ ∂gφr −vr gφθ + vφ (−gr θ + cot θ (gφφ − gθθ )) φφ vφ r ∂θ − r sin θ ∂φ − vr ∂r − r ∂θ + r ¡ ∂grr ¢ ¡ ¢ ∂gr φ ∂gr φ −vr (gφr + gr φ ) + vθ (gθ φ − gφθ − cot θgr φ ) ∂gr θ vr r sin θ ∂φ − ∂r − vθ r ∂θ − r sin θ ∂φ + r ¡ ∂gθr ¢ ¡ ¢ ∂gθφ ∂gθφ −vr (gθφ + cot θgφr ) − vθ (gr φ + cot θ (gφθ + gθφ )) ∂gθ θ vr − − − v + θ r ∂θ r sin θ ∂φ ∂r r sin θ ∂φ r ¢ ¡ ¢ ¡ ∂gφr ∂gφφ ∂gφφ ∂gφθ vr (grr − gφφ + cot θgθr ) + vθ (gr θ + cot θ (gθ θ − gφφ )) vr − ∂r − vθ r ∂θ − r sin θ ∂φ + r sin θ ∂φ r ρhe r ρhe θ ρheφ ρK
[A.9]
172
MILLER AND COLELLA
Again, the det g source term, not included above, is unaffected by the transformation of variables. The G diffusion term is described separately below. The geometric source terms in the vector s [cf. (28)] corresponding to the primitive variables q, but omitting the direction-dependent stress terms, are
−ρ(2vr + cot θ vθ ) 2 2 vθ + vφ + (2σrr − σθ θ − σφφ + cot θ σr θ )/ρ
−vr vθ + cot θ v 2 + (3σr θ + cot θ (σθ θ − σφφ ))/ρ φ −vr vφ − cot θ vθ vφ + (3σr φ + 2 cot θ σθ φ )/ρ (vr (σθ θ + σφφ ) + vθ (−σθr + cot θ σφφ ) − vφ (σφr + cot θ σθ φ ))/ρ 0 0 0 vr gθr + vθ gθ θ + vφ gθ φ −vr grr − vθ gr θ − vφ gr φ 0 vr gφr + vθ gφθ + vφ gφφ cot θ (vr gφr + vθ gφθ + vφ gφφ ) 1 . s= −vr grr − vθ gr θ − vφ gr φ − cot θ (vr gθr + vθ gθ θ + vφ gθ φ ) r ¡ p ¢ ¡ ¢ p p p vθ Fr θ + Fθr + vφ Fr φ + Fφr ¡ p ¢ vθ F − F p + vφ F p + cot θ vφ F p ρ θθ θφ φr rr ¡ ¢ p p p p vθ Fφθ + vφ Fφφ − Frr − cot θ vφ Fθr ¡ p ¢ p p p vθ Fθ θ − Frr + vφ Fφθ + cot θ vφ Fr φ ¡ p ¢ ¡ ¢ −vθ F + F p + cot θ vφ F p + F p ρ rθ θr φθ θφ ¡ p ¢ p p p −vθ Fφr − vφ Fr θ + cot θ vφ Fφφ − Fθ θ ¡ p ¢ p p p vθ Fθ φ + vφ Fφφ − Frr − cot θ vφ Fr θ ¡ p ¢ p p p −vθ F − vφ F + cot θ vφ F − F ρ rφ θr φφ θθ ¡ p ¡ p p ¢ p ¢ −vφ Fφr + Fr φ − cot θ vφ Fφθ + Fθ φ 0 (A.10) The stress evolution equations, used in the predictor steps of the method, are (nongeometric source terms omitted)
σ er σ er σ er σ er ∂ ∂ ∂ ∂ σ eθ σ eθ + vr σ eθ + vθ σ eθ + vφ ∂t ∂r r ∂θ r sin θ ∂φ σ eφ σ eφ σ eφ σ eφ vr vr vr Ar θ Ar φ Arr ∂ ∂ ∂ vθ = Aθr · vθ + Aθ θ · vθ + A θ φ · ∂r r ∂θ r sin θ ∂φ Aφr vφ Aφθ vφ Aφφ vφ
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
173
[2(vθ σr θ + vφ σr φ ) + vr (Ar θ )r θ − vθ (Ar θ )rr + (vr + vθ cot θ )(Ar φ )r φ − vφ (Ar φ )rr − cot θ vφ (Ar φ )r θ ] [−vθ (σrr − σθ θ ) + vφ (σθ φ + cot θ σr φ ) + vr (Ar θ )θ θ − vθ (Ar θ )θr + (vr + vθ cot θ)(Ar φ )θ φ − vφ (Ar φ )θr − cot θ vφ (Ar φ )θ θ ] σ − v (σ − σ + cot θ σ ) + v (A ) − v (A ) [v θ θφ φ rr φφ rθ r rθ rθ θ r θ φr + (vr + vθ cot θ )(Ar φ )φφ − vφ (Ar φ )rr − cot θ vφ (Ar φ )φθ ] [−vθ (σrr − σθ θ ) + vφ (σθ φ + cot θ σr φ ) + vr (Aθ θ )r θ − vθ (Aθ θ )rr + (vr + vθ cot θ)(Aθ φ )r φ − vφ (Aθ φ )rr − cot θ vφ (Aθ φ )r θ ] σ − v cot θσ ) + v (A ) − v (A ) + (v + v cot θ ) [−2(v 1 θ rθ φ θφ r θθ θθ θ θ θ θr r θ + r × (Aθ φ )θ φ − vφ (Aθ φ )θr − cot θ vφ (Aθ φ )θ θ ] [−vθ σr φ − vφ (σr θ + cot θ (σθ θ − σφφ )) + vr (Aθ θ )φθ − vθ (Aθ θ )φr + (vr + vθ cot θ)(Aθ φ )φφ − vφ (Aθ φ )φr − cot θ vφ (Aθ φ )φθ ] [vθ σθ φ − vφ (σrr − σφφ + cot θ σr θ ) + vr (Aφθ )r θ − vθ (Aφθ )rr + (vr + vθ cot θ )(Aφφ )r φ − vφ (Aφφ )rr − cot θ vφ (Aφφ )r θ ] [−vθ σr φ − vφ (σr θ + cot θ (σθ θ − σφφ )) + vr (Aφθ )θ θ − vθ (Aφθ )θr + (vr + vθ cot θ)(Aφφ )θ φ − vφ (Aφφ )θr − cot θ vφ (Aφφ )θ θ ] [−2vφ (σr φ + cot θσθφ ) + vr (Aφθ )φθ − vθ (Aφθ )φr + (vr + vθ cot θ ) × (Aφφ )φφ − vφ (Aφφ )φr − cot θ vφ (Aφφ )φθ ]
The g relaxation term −D(∇ ×∇ ×g T )T transforms in spherical coordinates as ((∇ × ∇ × g T )T )rr =
((∇ × ∇ × g T )T )r θ =
((∇ × ∇ × g T )T )r φ =
∂ 2 gr φ ∂ 2 gr θ ∂gr θ ∂ 2 grr ∂ 2 grr 2∂gθr + − 2 2 − 2 2+ 2 + 2 2 r sin θ∂r ∂φ r sin θ ∂φ r ∂r ∂θ r ∂θ r ∂θ r ∂θ ∂gθ θ 2∂gφr ∂gr φ cot θ ∂grr ∂gφφ − + 2 + − − r ∂r r ∂r r sin θ ∂φ r 2 sin θ ∂φ r 2 ∂θ 2grr gθ θ cot θ gr θ 2cot θ gθr gφφ cot θ∂gr θ + − 2 + 2 − 2 + + r ∂r r r r r2 r2 (A.12a) ∂gφφ ∂ 2 gr θ ∂ 2 gr θ 2∂gr θ ∂ 2 grr ∂ 2 gr φ − − 2 2 2 − 2 − + 2 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂ φ r ∂θ r ∂r cot θ ∂gr φ 2∂gφθ ∂gθ φ gr θ ∂gθr + 2 + − + 2 − r ∂r r sin θ ∂φ r 2 sin θ ∂φ r 2 sin θ ∂φ r cot θ gφφ cot θ gθ θ − + (A.12b) 2 r r2 ∂ 2 gr φ ∂ 2 gr φ 2∂gθ φ ∂ 2 gr θ ∂ 2 grr ∂gφθ − − + − 2 + 2 r 2 sin θ∂θ∂φ r 2 ∂θ 2 r sin θ ∂r ∂φ ∂r 2 r ∂θ r ∂θ 2∂gr φ cot θ ∂gr θ cot θ ∂gr φ ∂gθ θ gr φ ∂gφr − − 2 − − 2 + 2 − r ∂r r ∂r r sin θ ∂φ r 2 ∂θ r sin θ ∂φ r gr φ cot θ gθ φ cot θ gφθ + 2 2 + + (A.12c) 2 r r2 r sin θ
174
MILLER AND COLELLA
((∇ × ∇ × g T )T )θr =
((∇ × ∇ × g T )T )θ θ =
((∇ × ∇ × g T )T )θ φ =
((∇ × ∇ × g T )T )φr =
((∇ × ∇ × g T )T )φθ =
((∇ × ∇ × g T )T )φφ =
∂ 2 gθr ∂ 2 gθr ∂gθ θ ∂ 2 gθ φ ∂ 2 gθ θ 2∂grr − 2 2 − + 2 + − 2 r sin θ∂r ∂φ r sin θ ∂φ 2 r ∂r ∂θ r 2 ∂θ 2 r ∂θ r ∂θ ∂gr θ 2cot θ ∂gφr ∂gθ φ cot θ ∂gθr cot θ ∂gφφ + + 2 + 2 − − 2 r ∂r r sin θ ∂φ r sin θ ∂φ r ∂θ r ∂r cot θ gφφ gr θ gθr cot θ gθ θ cot θ∂gθ θ + 2 + 2 2 − + (A.12d) + r ∂r r r2 r2 r sin θ ∂ 2 gθr ∂ 2 gθ φ ∂grr ∂ 2 gθ θ ∂ 2 gθ θ 2∂gθ θ + 2 + − − 2 2 − 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂φ 2 r ∂r r ∂r cot θ ∂gθ φ cot θ ∂gφφ ∂gr φ 2cot θ∂gφθ + 2 − + 2 + 2 r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r sin θ ∂φ gθ θ cot θ gr θ gθ θ gφφ gφφ (A.12e) + 2 − 2 − 2 2 + 2 2 + r r r2 r sin θ r sin θ ∂ 2 gθ φ ∂ 2 gθ φ 2∂gθ φ ∂ 2 gθ θ ∂ 2 gθr 2∂gr φ − − − + − 2 r 2 sin θ∂θ∂φ r 2 ∂θ 2 r sin θ ∂r ∂φ ∂r 2 r ∂θ r ∂r cot θ ∂gφθ ∂gr θ cot θ ∂gθ φ cot θ ∂gφr cot θ∂gθ θ − + 2 − − − 2 2 2 r sin θ∂φ r ∂θ r sin θ ∂φ r ∂θ r ∂r gθ φ gφθ gθ φ cot θ gr φ gφθ (A.12f) − 2 + 2 + 2 2 + 2 2 − r r r2 r sin θ r sin θ ∂ 2 gφφ ∂ 2 gφθ ∂gφθ ∂ 2 gφr ∂ 2 gφr ∂gr φ + + 2 + − 2 2 − 2 2 2 r sin θ∂r ∂φ r sin θ ∂φ r ∂r ∂θ r ∂θ r ∂θ r ∂r ∂gφφ cot θ ∂gφr 2∂grr cot θ ∂gφθ 2cot θ∂gθr + 2 − − 2 + − 2 r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r sin θ ∂φ r ∂r gφr cot θ gθ φ gr φ cot θ∂gθ φ cot θ gφθ + (A.12g) + 2 + 2 2 + + r ∂r r r2 r2 r sin θ ∂ 2 gφr ∂ 2 gφφ ∂gr φ ∂ 2 gφθ ∂ 2 gφθ 2∂gφθ + + 2 − − − 2 2 2 2 2 r ∂r ∂θ ∂r r sin θ ∂θ ∂φ r sin θ ∂φ r ∂θ r ∂r 2cot θ ∂gθ θ 2∂gr θ cot θ ∂gθ φ 2gθ φ cot θ∂gφφ − 2 − 2 + − 2 + 2 r sin θ∂φ r sin θ ∂φ r sin θ ∂φ r 2 ∂θ r gφθ gθ φ 2cot θ gr φ + 2 2 + 2 2 + (A.12h) r2 r sin θ r sin θ ∂ 2 gφθ 2 r sin θ∂θ∂φ
∂ 2 gφr ∂gr θ ∂ 2 gφφ ∂ 2 gφφ + + 2 − 2 2 r ∂θ r sin θ ∂r ∂φ ∂r 2 r ∂θ ∂grr cot θ ∂gφθ cot θ ∂gφφ cot θ ∂gθ θ 2∂gφφ + − 2 − + − r ∂r r ∂r r sin θ ∂φ r 2 ∂θ r 2 ∂θ gφφ cot θ∂gθr gθ θ + 2 2 − 2 2 + r ∂r r sin θ r sin θ −
(A.12i)
We have implemented these equations in 1D, direction r , with only slight modifications to the strategy described for Cartesian geometry. Schematically, we represent the overall system of equations in the form ∂ AF(U ) ∂ H (U ) ∂U + + = G(q, r ) + S(q). ∂t ∂V ∂r
(A.13)
175
EULERIAN GODUNOV METHOD FOR SOLID MECHANICS
Here we distinguish between the area-weighted volumetric flux terms, AF, and the spatial flux terms H . The geometric source terms are represented by G(q, r ), and the plastic source terms are S(q). Note that in strict 1D-r flow, there is no angular dependence to any flow variable, and therefore terms proportional to cot θ (for example) vanish identically. As in the Cartesian case, we solve the time-centered edge Riemann problems to deduce ∗,n+1/2 single-valued time-centered edge states Ui+1/2 . These edge states are then used to construct the flux terms F and H , which are used to compute a preliminary update U˜ n+1 via the difference scheme ¢ ¢ 1t ¡ 1t ¡ ∗ ∗ ∗ ∗ − Ai−1/2 Fi−1/2 Ai+1/2 Fi+1/2 − Hi+1/2 − Hi−1/2 . U˜ in+1 = Uin − Vi 1ri
(A.14)
Next, we modify the preliminary update by inclusion of the geometric source terms. This is made second-order using a predictor-corrector strategy, ¡ ¢ 0 U˜ i = U˜ in+1 + 1t G qin , σin µ ¶n ∂σ (q˜ 0 − q n ) σ0 = σn + ∂q ¢ 1t ¡ ¡ n n ¢ 00 U˜ i = U˜ in+1 + G qi , σi + G(q˜ i0 , σi0 ) . 2
(A.15)
The plastic source terms are then evaluated at the half time step, giving the final result ¢ 1¡ n qi + q˜ i00 2 = U˜ i0 + 1t S(q¯ i ).
q¯ i = Uin+1
(A.16)
ACKNOWLEDGMENTS G.M. benefitted from correspondences with B. Plohr (SUNY at Stony Brook), V. Barcilon (University of Chicago), and M. Ortiz (Caltech).
REFERENCES 1. J. B. Bell, P. Colella, and M. L. Welcome, Conservative front-tracking for inviscid compressible flow, in AIAA 10th Comput. Fluid Dynamics Conf., Honolulu, Hawaii, 1991. 2. J. B. Bell, C. N. Dawson, and G. R. Shubin, An unsplit, higher order Godunov method for scalar conservation laws in multiple dimensions, J. Comput. Phys. 74, 1 (1988). 3. F. G. Blake, Spherical wave propagation in solid media, J. Acoust. Sci. Am. 24, 211 (1952). 4. I.-L. Chern and P. Colella, A Conservative Front Tracking Method for Hyperbolic Conservation Laws, Technical Report (LLNL, 1987), UCRL-97200. 5. P. Colella, A direct Eulerian MUSCL scheme for gas dynamics, SIAM J. Sci. Stat. Comput. 6, 104 (1985). 6. P. Colella, Multidimensional upwind methods for hyperbolic conservation laws, J. Comput. Phys. 87, 171 (1990). 7. P. Colella and P. R. Woodward, The piecewise parabolic method (PPM) for gas dynamical simulations, J. Comput. Phys. 54, 174 (1984). 8. S. K. Godunov, Symmetric form of magnetohydrodynamics equations, Chislennye Metody Mekhaniki Sploshnoi Sredy, Novosibirsk 3, 26 (1972), in Russian.
176
MILLER AND COLELLA
9. S. K. Godunov and E. I. Romensky, Thermodynamics, conservation laws, and symmetric forms of differential equations in mechanics of continuous media, in Computational Fluid Dynamics Review 1995 (Wiley, New York, 1995), pp. 19–31. 10. R. Hill, The Mathematical Theory of Plasticity (Oxford University Press, Oxford, U.K., 1950). 11. Mindy Fruchtman Lai, A Projection Method for Reacting Flow in the Zero Mach Number Limit, Ph.D. thesis (Dept. of Mechanical Engineering, Univ. of California, Berkeley, 1993). 12. E. H. Lee and D. T. Liu, Finite-strain elastic-plastic theory with application to plane-wave analysis, J. Appl. Phys. 38, 19 (1967). 13. J. Lubliner, Plasticity Theory (McMillan, New York, 1990). 14. J. Mandel, Plasticit´e Classique et Viscoplasticit´e (Springer-Verlag, New York, 1972). 15. G. H. Miller and E. G. Puckett, A high-order Godunov method for multiple condensed phases, J. Comput. phys. 128, 134 (1996). 16. R. B. Pember, J. B. Bell, P. Colella, W. Y. Crutchfield, and M. L. Welcome, An adaptive Cartesian grid method for unsteady compressible flow in irregular regions, J. Comput. Phys. 120, 278 (1995). 17. L. Petzold, Differential/algebraic equations are not ODEs, SIAM J. Sci. Stat. Comput. 3, 367 (1982). 18. J. E. Pilliod, A Second-Order Unsplit Method for Modeling Flames in Two-Dimensional Compressible Flow, Ph.D. thesis (University of California, Davis, 1996). 19. B. J. Plohr and D. H. Sharp, A conservative formulation for plasticity, Adv. Appl. Math. 13, 462 (1992). 20. K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. DeZeeuw, A solution-adaptive upwind scheme for ideal magnetohydrodynamics, J. Comput. Phys. 154, 284 (1999). 21. W. J. Rider, Filtering non-solenoidal modes in numerical simulations of incompressible flow, Int. J. Numer. Meth. Fluids 28, 789 (1998). 22. J. Saltzman, An unsplit 3D upwind method for hyperbolic conservation laws, J. Comput. Phys. 115, 53 (1994). 23. J. C. Simo and T. J. R. Hughes, Computational Inelasticity (Springer-Verlag, New York, 1997). 24. G. F. Smith and R. S. Rivlin, The strain-energy function for anisotropic elastic materials. Trans. Am. Math. Soc. 88, 175 (1958). 25. J. A. Trangenstein and P. Colella, A higher-order Godunov method for modeling finite deformation in elasticplastic solids, Comm. Pure Appl. Math. 44, 41 (1991). 26. B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys. 32, 101 (1979). 27. M. L. Wilkins, Calculation of elastic-plastic flow, in Methods in Computational Physics, edited by B. Alder, S. Fernbach, and M. Rotenberg (Academic, New York, 1964), Vol. 3, p. 211.