A bodner-partom visco-plastic dynamic sphere benchmark problem

V International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2013 S. Idelsohn, M. Papadrakakis and B. Schrefler (Eds) A BODNER-PARTOM VISCO-PLASTIC DYNAMIC SPHERE BENCHMARK PROBLEM Brandon M. Chabaud∗, Jerry S. Brock†and Todd O. Williams† ∗Los Alamos National Laboratory (LANL) Los Alamos, NM 87545, United States email: c[email protected]v †Los Alamos National Laboratory (LANL) Los Alamos, NM 87545, United States email: jsbroc[email protected]v, [email protected]v Key words: spherical shell, analytic solution, Bodner-Partom model Abstract. Developing benchmark analytic solutions for problems in solid and fluid mechanics is very important for the purpose of testing and verifying computational physics codes. Our primary objective in this research is to obtain a benchmark analytic solution to the equation of motion in radially symmetric spherical coordinates. An analytic solution for the dynamic response of a sphere composed of an isotropic visco-plastic material and subjected to spherically symmetric boundary conditions is developed and implemented. The radial displacement uis computed by solving the equation of motion, a linear secondorder hyperbolic PDE. The plastic strains εp rr and εp θθ are computed by solving two nonlinear first-order ODEs in time. We obtain a solution for uin terms of the plastic strain components and boundary conditions in the form of an infinite series. Computationally, at each time step, we set up an iteration scheme to solve the PDE-ODE system. The linear momentum equation is solved using the plastic strains from the previous iteration, then the plastic strain equations are solved numerically using the new displacement. We demonstrate the accuracy and convergence of our benchmark solution under spatial mesh, time step, and eigenmode refinement. 1 INTRODUCTION In this paper we derive an analytic solution for the geometrically linear dynamic response of a sphere for the purpose of comparing it to the results from computational physics codes. An example of a spherical solid mechanics problem is the Blake problem, which has been of considerable long-term [1, 2] and recent [3] interest to the scientific community. The Blake problem consists of a spherical inclusion within an infinite medium. 1 Bodner-Partom visco-plastic dynamic sphere Benchmark problem 671 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams In contrast, the dynamic sphere problem considered here develops an analytic solution for the dynamic response of a spherically symmetric shell composed of a visco-plastic, isotropic material and subjected to prescribed time-dependent boundary conditions imposed at both the inner and outer radii (riand ro, respectively). The boundary conditions are arbitrary and may be any one of the following: applied radial stresses, applied radial strains, or applied radial displacements. The physical problem under consideration is illustrated in the figure below. Figure 1: 2-D schematic of dynamic sphere problem. In spherical coordinates the stress-strain constitutive relation for a visco-plastic isotropic material is given by the linear system [4] σ=Cεe(1) =C(ε−εp) (2) =Cε −Cεp(3) =Cε +λp,(4) λp=−Cεp.(5) Here σis the total stress tensor, λpis the plastic stress tensor, εis the total linear strain tensor, εeis the elastic strain tensor, εpis the plastic strain tensor, and Cis the constant fourth-order stiffness tensor. For a spherically symmetric shell, the shear components of the strain (εrθ,εrϕ,εθr,εθϕ,εϕr,εϕθ) are all identically zero, so the only stresses on the shell are normal stresses (σrr,σθθ,σϕϕ). Also, by the isotropic symmetry of the material, the tangential stresses and strains satisfy σθθ =σϕϕ,λp θθ =λp ϕϕ,εθθ =εϕϕ, and εp θθ =εp ϕϕ. The governing equations of the dynamic sphere problem come from the balance laws of continuum mechanics. We formulate these laws in material coordinates, so the balance of mass equation is automatically satisfied. Additionally, we assume that the shell is kept at constant temperature, so the balance of energy equation can also be neglected. The initial-boundary value problem (IBVP) to be solved, which comes from the equation of motion, is c−2 rutt =urr +2 rur−2 r2u+c−2 rfλ;t>0,r i≤r≤ro(6) fλ=1 ρ0 [∂λp rr ∂r +2 r(λp rr −λp θθ)]; t>0,r i≤r≤ro(7) 2 672 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams u(r, 0)=0,u t(r, 0) = 0; ri≤r≤ro(8) αBur(rB,t)+βBu(rB,t) = BCB(t)−dλp rr(rB,t); B=i, o;t>0 (9) In this system u(r, t) is the radial displacement, and in (8) we have set initial displacement and velocity equal to zero. fλis the contribution of the plastic stress to the acceleration of the material. This term is given by (7), where fλis defined in terms of the radial plastic stress λp rr and the tangential plastic stress λp θθ. From (5), we see that the plastic stresses are given in terms of the plastic strains by λp rr =−Crrrrεp rr −2Crrθθεp θθ,λ p θθ =−Crrθθεp rr −(Crrrr +Crrθθ)εp θθ.(10) Boundary conditions are given by (9), where B=irefers to conditions imposed on the inner surface and B=orefers to conditions imposed on the outer surface. BCi(t) and BCo(t) are prescribed time-dependent boundary conditions imposed on the inner and outer surfaces, respectively. The speed of sound through the material is cr=Crrrr ρ0, and ρ0is the initial density. For later times, the density ρis given by the linearized Lagrangian mass equation ρ=ρ0/(1 + εrr +2εθθ).(11) The values of the material constants αi,αo,βi, and βodepend on whether displacement (Dirichlet), strain (Neumann), or stress (Robin) boundary conditions are imposed [5, 6, 7]. The parameter dthat appears in (9) is equal to zero when Dirichlet or Neumann conditions are imposed and is equal to unity when Robin conditions are imposed. The system of equations given by (6)-(10) is not complete. Evolution equations must be given for the plastic strain components. In this work the plastic strain tensor evolves according to the Bodner-Partom (BP) model for visco-plasticity [8]. For a radially symmetric spherical shell problem, the equations are ˙εp rr =Srrγ, ˙εp θθ =Sθθγ; (12) εp rr(r, 0)=0,ε p θθ(r, 0) = 0; (13) Srr =σrr −1 3(σrr +2σθθ),S θθ =σθθ −1 3(σrr +2σθθ); (14) In this system ˙εp rr =∂εp rr ∂t and ˙εp θθ =∂εp θθ ∂t are the plastic strain rates. (13) gives initial conditions for the ODEs (12). According to the BP model, these rates are proportional to the stress deviator components Srr and Sθθ, which are defined by (14). The factor γthat appears in (12) is the proportionality coefficient for the BP model. We refer the reader to [8] for the details of how to compute it. 2 ANALYTIC SOLUTION Our primary objective is to obtain an analytic expression for u(r, t) in terms of λp rr and λp θθ that solves the system (6)-(9). In order to obtain a solution, we introduce the variable 3 673 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams w(r, t)=r1/2u(r, t). The transformed equation of motion is rc−2 rwtt =rwrr +wr−µ2 rw+r3/2f, (15) where µ=3/2 and f=c−2 rfλ. We decompose winto two terms w=w+˜w. We take w(r, t)=γ0(t)+γ1(t)rand choose the coefficients γ0and γ1so that wsatisfies the transformed boundary condition equations. This amounts to solving the linear system β∗ oα∗ o+β∗ oro β∗ iα∗ i+β∗ iriγ0(t) γ1(t)=ki(t) ko(t)(16) where α∗ i=αi √ri ,α ∗ o=αo √ro ,β ∗ i=βi √ri−αi 2r3/2 i ,β ∗ o=βo √ro−αo 2r3/2 o (17) are transformed boundary condition coefficients and ki(t) = BCi(t)−dλp rr(ri,t),k o(t) = BCo(t)−dλp rr(ro,t) (18) are transformed boundary conditions. From this system and the form chosen for w, it is easy to show that wcan be expressed as w(r, t)=[γ0,o +γ1,or]ko(t)+[γ0,i +γ1,ir]ki(t),(19) where γ0,o,γ1,o,γ0,i, and γ1,i are constants. The transformed IVBP solution ˜w(r, t) satisfies rc−2 r˜wtt =L(˜w)+ ˜ f, (20) ˜w(r, 0)=w(r, 0) −w(r, 0),˜wt(r, 0)=wt(r, 0) −wt(r, 0),(21) α∗ B˜wr(rB,t)+β∗ B˜w(rB,t)=0,B=i, o (22) where L(w)=rwrr +wr−µ2 rw(23) is the Sturm-Liouville operator and ˜ f=L(w)−c−2 rrwtt +r3/2f. (24) Thus we have reduced the problem of finding the displacement uto the problem of solving the homogeneous IBVP (20)-(22). The solution to (20)-(22) has the form ˜w(r, t)= ∞  n=1 an(t)ψn(r),(25) 4 674 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams where ψn(r) are eigenfunctions which are computed from solving the autonomous problem (with ˜ f≡0) and an(t) are time-dependent coefficients that satisfy the second-order-intime ODE d2an dt2+c2 rλnan=Fn.(26) Here λnis the eigenvalue corresponding to eigenfunction ψnand Fn(t) is a source function given by Fn(t)=c2 rro ri ˜ f(r, t)ψn(r)dr ro rirψ2 n(r)dr .(27) The eigenfunctions are given by ψn(r)=c1,nJµ(rλn)+c2,nYµ(rλn),(28) where Jµand Yµare Bessel functions of order µof the first and second kind, respectively, and c1,n and c2,n are constants that depend on eigenvalue λn. The solution to the ODE (26) is given by an(t)=c3,n cos crλnt+c4,n sin crλnt+1 cr√λnt 0 Fn(τ) sin crλn(t−τ)dτ, (29) where the constants c3,n and c4,n are given by c3,n =ro rir˜w(r, 0)ψn(r)dr ro rirψ2 n(r)dr ,c 4,n =1 cr√λnro rir˜wt(r, 0)ψn(r)dr ro rirψ2 n(r)dr .(30) For the derivation of this solution, we refer the reader to [9, 10]. Thus we have an analytic expression for the displacement u=r−1/2(w+˜w) in terms of λp rr and λp θθ solving the original IBVP (6)-(9). 3 COMPUTATION OF ANALYTIC SOLUTION Our ultimate objective in developing the analytic solution derived in the previous section is to compare it with results of computational physics codes. The analytic solution of the dynamic sphere problem is cast as an infinite series solution. In practice, however, we truncate the Fourier-Bessel series solution (25) after a finite number of terms. We also notice that our benchmark solution requires us to evaluate several space and time integrals, so computationally we must apply spatial and time meshes. Additionally, the plastic stress and strain equations (10)-(14) must be solved numerically. There is no analytic expression for these quantities. In this section we describe the algorithm we use to compute all relevant quantities at each time step. For any integer n≥1, we assume that the solution has been computed up to time tn−1. At each time step, we iterate between the equation of motion (6) and the plasticity equations (10)-(14) to obtain a solution. For convenience we define the plastic strain increment to be ∆εp(tn)=εp(tn)−εp(tn−1). The procedure is described in the following algorithm. 5 675 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams Algorithm 1 (Decoupling Equations at Time Step n)Compute u,εp, and λpat time tn. Set εp,(0) rr (tn)=εp rr(tn−1),εp,(0) θθ (tn)=εp θθ(tn−1),λp,(0) rr (tn)=λp rr(tn−1),λp,(0) θθ (tn)=λp θθ(tn−1). for iteration k=1,2,··· do Compute u(k)(tn)with analytic solution using λp,(k−1) rr (tn),λp,(k−1) θθ (tn)in plastic force term. Substitute u(k)and λp,(k−1) into (12) and solve for εp,(k) rr (tn),εp,(k) θθ (tn)using explicit forward Euler finite difference method. Substitute εp,(k) rr (tn),εp,(k) θθ (tn)into (10) to obtain λp,(k) rr (tn),λp,(k) θθ (tn). If ||∆εp,(k)(tn)−∆εp,(k−1)(tn)||L∞≤tol||∆εp,(k−1)(tn)||L∞, Set u(tn)=u(k)(tn),εp(tn)=εp,(k)(tn),λp(tn)=λp,(k)(tn). Exit loop. end for In this algorithm, the L∞norms are computed over the spatial mesh and tol is a userdefined tolerance that determines when the iteration scheme converges. We point out here that the analytic solution for ugiven in the previous section contains spatial integrals of spatial derivatives of λp rr and temporal integrals of time derivatives of λp rr. It is possible to use finite difference approximations to evaluate these derivatives. However, the best approach, leading to the fastest convergence, is to use integration by parts to remove as many plastic stress derivatives as possible from the analytic solution. We also note that when evaluating the plastic stress integrals, sixth-order Newton-Cotes or a similar high-order numerical quadrature scheme must be used to compute the integrals accurately. A lower-order quadrature scheme, such as the trapezoid rule, will result in nonphysical oscillations. 4 SELF-CONVERGENCE ANALYSIS In this section we consider a test problem and examine how the computed solution converges with respect to space and time meshes as well as with respect to number of terms (eigenmodes) taken in the series solution. We impose stress boundary conditions on the inner and outer radii, so d= 1 in (9). We take the interior of the shell to be a void, so BCi(t) = 0, and on the outer radius we take a time-varying smooth jump stress condition. Figure 2 gives a plot of the outer radius stress boundary condition. The curve takes the form of a hyperbolic tangent function. The exact form and parameters used for this boundary condition can be found in Section 2 of [11]. In all of our simulations, we take inner radius ri= 1 cm, outer radius ro= 2 cm, and we run our simulations out to a final time of T=3.5 microseconds. The constants appearing in the boundary conditions (9) are given by αi=αo=Crrrr,β i=2 ri Crrθθ,β o=2 ro Crrθθ.(31) The stiffness tensor coefficients that appear in the problem are computed in terms of 6 676 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams Figure 2: Stress boundary condition on outer radius. engineering constants [12]. For an isotropic material, the coefficients are given by Crrrr =K+4 3G, Crrθθ =K−2 3G, (32) where the constants Kand Gare the bulk and shear moduli of the material, respectively. Table 1 gives the parameter values, which are taken from Case (d) of [11], that we use in our simulations. The parameter values given in the table do not correspond to any real material. Our only interest here is verification of physics codes. Table 1: Parameters used in convergence simulations. Parameter Value ri1 cm ro2 cm ρ01000 kg/m3 K100 GPa G60 GPa In order to assess the convergence of our benchmark analytic solution, we will take both a qualitative and quantitative approach. Qualitatively, we will simply plot relevant physical quantities for different size meshes and show that the plots converge. To demonstrate quantitative convergence, we will compute L∞norms of percent errors with respect to extremely refined reference solutions and show that these norms are small. For example, the L∞norm of the percent error between the displacement uand a reference 7 677 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams displacement uref is computed by the expression % error = ||u−uref ||L∞ ||uref ||L∞×100%,(33) where we take the L∞norm either in time over the interval (0,T) along shell boundaries or in space over the spherical shell (ri,r o). The quantities of interest in our convergence analysis are displacement, velocity, and radial and tangential total strain, plastic strain, and stress. We begin by considering spatial convergence. To assess convergence qualitatively, we ran simulations with nl = 500 eigenmodes, nt = 8000 equally spaced time intervals, and uniform spatial meshes with varying numbers of intervals. The top two plots in Figure 3 Figure 3: Top Left: Spatial convergence plot of utvs. rthrough radial shell at time T=3.5 microseconds. Top Right: Zoomed-in plot of ut. Bottom: Log-log plot of L∞(ri,r o) percent error of utthrough shell vs. nr. 8 678 Brandon M. Chabaud, Jerry S. Brock and Todd O. Williams show the velocity utvs. radial position rthrough the shell at final computation time T=3.5 microseconds for uniform spatial meshes of size nr = 2000, 2500, 3000, 3500, 4000, and 4500 intervals. The zoomed-in plot at the top right of the figure shows that the velocity profile is converging as the mesh is refined. The bottom plot in Figure 3 gives a log-log plot of the L∞(ri,r o) norm of the percent error of utvs. nr for nr varying from 500 to 5500 in increments of 500. A mesh of nr = 6000 intervals was used to compute a reference solution for (33). We notice that the percent error corresponding to nr = 3500 is approximately 10−4.5%. Of all the quantities of interest, the one with the highest error through the shell is radial plastic strain, with an L∞(ri,r o) norm of percent error of approximately 10−3.5% for nr = 3500. Computing L∞(0,T) norms of percent errors on the inner and outer radii gives even smaller errors. Therefore, because these errors are so small, we take nr = 3500 as an appropriate spatial resolution for a “converged” solution. To assess temporal convergence, we ran simulations with nr = 3500 uniform spatial intervals, nl = 500 eigenmodes, and uniform time meshes of varying sizes. The top two plots in Figure 4 show the radial plastic strain εp rr vs. time tat the inner radius rifor time meshes of size nt = 5500, 6000, 6500, 7000, and 7500 uniform time steps. The zoomed-in plot at the top right of the figure shows that the radial plastic strain profile is converging as smaller step sizes are used. The bottom plot in Figure 4 gives a log-log plot of the L∞(0,T) norm of the percent error of εp rr at rivs. nt for nt varying from 2500 to 7500 in increments of 500. A mesh of nt = 8000 intervals was used to compute a reference solution for (33). We notice that the percent error corresponding to nt = 6500 is less than 0.01%. Radial plastic strain has the largest L∞(0,T) norm of percent error at ri of all the quantities of interest. We obtain similar results when we compute L∞(0,T) norms of percent errors at the outer radius ro. Computing L∞(ri,r o) norms of percent errors at final time T, we find once again that εp rr exhibits the largest error, this time with error approximately 0.1%. Because these errors are so small, we take nt = 6500 as an appropriate temporal resolution for a “converged” solution. Finally, we assess eigenmode convergence. We ran simulations with nr = 3500 uniform spatial intervals, nt = 6500 uniform time steps, and varying numbers of eigenmodes. The top two plots in Figure 5 show the radial strain εrr vs. time tat the outer radius rofor nl = 100, 150, 200, 250, and 300 eigenmodes. The zoomed-in plot at the top right of the figure shows that the radial strain profile is converging as more eigenmodes are used. The bottom plot in Figure 5 gives a log-log plot of the L∞(0,T) norm of the percent error of εrr at rovs. nl for nl varying from 100 to 450 in increments of 50. A reference solution with nl = 500 eigenmodes was used to compute the percent error in (33). We notice that the percent errors for all of the eigenmodes shown in the log-log plot are less than 10−5%. For nl = 300 (fourth mark from the right in the log-log plot), the error is less than 10−6.5%. Of all the quantities of interest, velocity has the largest L∞(0,T) norm of percent error at ro, slightly less than 10−5% when nl = 300. In the L∞(0,T) norm at ri, velocity still has the largest error, approximately 10−3.6% when nl = 300. In the L∞(ri,r o) norm at final time T, the quantity with the largest error is radial stress, with 9 679