Visco-elastic relaxation with a van der waals type stress

Int. 1. Engng Sci. Vol. 32, No. 2. pp. 327-338, Printed in Great Britain. All rights reserved

0020-7225194

194

$6.00 + 0.00

Copyright@ 1994 PergamonPressLtd

VISCO-ELASTIC RELAXATION WITH A VAN DER WAALS TYPE STRESS E. BRUCE

PITMAN?

and YIGONG

NI

Department of Mathematics, State University of New York, Buffalo, NY 14214-3093, U.S.A. Abstract-Contemporary science models certain phase transitions as conservation laws with a non-monotone stress-strain relation. The system of conservation laws are often hyperbolic except in a region in phase space, where the system is elliptic or of composite type. It is possible to prove existence and uniqueness of solutions for specific models which exhibit such change of type, but notions of admissibility and the connections between the mathematical theory and engineering experiments are not fully understood. In this paper we model phase transition by a visco-elastic relaxation law, and describe numerical experiments on this system. In particular we examine special initial data which, in the context of the mixed system, gives rise to multiple solutions.

INTRODUCTION

Certain dynamic phase transitions are modeled by conservation laws with non-monotone stress-strain relations. The best known non-monotone system is the van der Waals gas at sufficiently low temperature. James [8] proposed a non-monotone dynamic model for crystal transitions in steel. Since then, many others have examined different aspects of dynamics involving non-monotone stress-strain relations. Important to this study is the work of Shearer who solved the Riemann problem for a van der Waals type model [ll] and showed that there exist two solutions to the Riemann problem for data in a certain range [12]. At the same time, Slemrod [13] was examining a regularized phase transition model which included viscositycapillarity terms. Affouf and Caflisch conducted a numerical study of shock-like solutions for the visco-capillarity system [3]. More recently Abeyaratne and Knowles [l, 21 have proposed a theory requiring additional constitutive information in the form of a kinetic relation and a nucleation criterion, and have shown how this information leads to a well-posed initial value problem and unique solution. A different kind of regularization has been proposed by Suliciu and co-workers, who introduce a visco-elastic regularization which allows for a pressure relaxation over time (see [6,14] and references cited there). Among other results, Suliciu [14] proposed solving the visco-elastic system numerically by using the method of characteristics, and Faciu [6] proved an L* bound on the difference Ip -prefl, where p is the pressure and pref is the reference (non-monotone) pressure relation. In this paper we employ a second order Godunov scheme to study the visco-elastic equations with a van der Waals type reference pressure. The scheme yields high resolution without numerical oscillation, and allows for a careful examination of Ip -prefl; our results suggest an L” bound on this difference. We then compare our visco-elastic calculations with calculations on the mixed hyperbolic-elliptic system using a simple Lax-Friedrichs scheme, and comment on a formal expansion result which indicates when solutions of the visco-elastic system may be close to solutions of the mixed system. The next section describes the governing systems of equations for the visco-elastic and mixed systems, and presents some of the mathematical facts we will use in this study. In Section 3, we describe our version of the Godunov scheme and in Section 4 show computational results. GOVERNING

EQUATIONS

The basic formulation of the equations describing the deformation of a uniform elastic bar with constant density in one dimension is to consider the motion Q, which maps the particles in tThis research supported in part by the Air Force Office of Scientific Research under Grant 900076. 327

328

E. B. PITMAN

and Y. NI

the bar from their original position (say the interval [0, l] at time t= 0)to new positions @(x, t) at time t. The strain in the bar is given by u = a,@ and the velocity of the particles in the bar is TV= a,@,. The equality of mixed partials is then equivalent to the conservation of mass, and mass and momentum balance are written

a,u - a,v = 0

(1)

a,v- a,p = 0.

(2)

Here p is the stress in the bar, usually considered as a function of strain, i.e. p =p(u). non-monotone van der Waals type stress-strain relation that we use here is given by p = p=‘(u)

The system (l-3)

= u3/3 - u.

A

(3)

is hyperbolic with eigenvalues

except for u E (-1, 1) when the system is elliptic. The Riemann problem for this system consists in solving (1,2,3) for x E (- m, ~4) with discontinuous initial data (u, V) = (u,, u,) for x 0. In order to study hyperbolic systems of PDEs, the mathematical theory requires imposition of an admissibility criterion, which determines which discontinuous solutions are to be allowed. Shearer [ll] solved the Riemann problem, imposing a chord condition for admissibility. Subsequent work by Shearer [12] and Slemrod [13] examined viscous and visco-capillarity admissibility criteria for the system (l-3). These criteria allow discontinuous solutions that are the limit, as E * 0, of the regularized system

a,u- a,v = BE a$ a,v- axp= E a$ - ~2

a$.

For the viscous case, B = 1 and A = 0; for the visco-capillarity case, B = 0 and 0
a,(p - E%) = -

(P

- Pref) t .

Here E* is (essentially) the dynamic Young’s modulus, and r is a relaxation time. The stress p is allowed to deviate from the reference stress pref, but this deviation is quenched over a timescale is always hyperbolic, with determined by the relaxation time r. The system (1,2,4) characteristic speeds Izo= 0 and A+ = f E. In fact the system is semi-linear, and, for fixed t > 0, standard existence and uniqueness results apply; see, for example, [5]. (REMARK:The original model proposed by Suliciu et al. involves a more complicated time derivative term; equation (4) is a mathematically tractable simplification they proposed.) Formally (1,2,4) approaches (1,2,3) in the limit t+ 0, but the precise nature of this convergence is not known. Let us describe an expansion procedure due to Liu (see [lo]) for systems with relaxation; this expansion is similar to the Chapman-Enskog expansion for the Boltzman equations. The principle idea here is to expand p about the pressure pref and examine the resulting asymptotics. At first order, p =prcf; substituting into (1,2), we recover the mixed system (1,2,3). Now writing p = pref + tp, , equation (4) gives -p, = (a#-

~2 a,+ = ((prep a+ - E* a,v)

Visco-elastic

relaxation

329

From equation (l), d,u = &v, so -p,

= ((p”?

3,~ - E* 3,~).

Substituting this expression into (2) yields the system d,u - a,v = 0 3,~ -

axPIe'= ta,((~*-(Pref)') a,v).

(5)

(6)

The system (5,6) is linearly stable whenever E* > (pref)‘, irrespective of the non-monotonicity of p”‘. Of course, this calculation is purely formal and obtaining error estimates for the expansion procedure is crucial. For example, when (pref)’ < 0, the first order system may be linearly ill posed and exhibit exponential blow up of Fourier modes; continuing the expansion to second order may not be justified. We are currently examining conditions under which an expansion procedure like the above can be justified. We mention that the magnitude of (prcf)’ in equations (4) and (6) depends on the solution of the particular initial value problem being considered, and one must exercise care when selecting a value of E. In the computations reported here, we set E = 3 which is large enough to satisfy the stability constraint for the problems we study. The stability constraint E* > (p”‘)’ can be derived for strictly hyperbolic systems following Whitham [17]. In [6], Faciu proved L* estimates for the difference Ip -prefl of (1,2,4), obtaining a bound which approaches zero as t+ 0. This bound was motivated in part by numerical results of [14] for the system (1,2,4) using a method of characteristics. In the next section we describe a second order Godunov method which suggests a pointwise bound on Ip -prefj.

NUMERICAL

METHOD

We give a brief description of the numerical techniques we use to solve (1,2,4), highlighting in particular the second order Godunov method for hyperbolic conservation laws. For further details about the Godunov method, see [l-5,16]. It is convenient to write (1,2,4) as a matrix system

a,u+AaJJ=R(U) where U = (u, v, p)‘,

R(U) is the vector

= (o,o2_(p-Pref) > T and

The superscript T denotes transpose. Divide the computational domain into I equal length subintervals of length AX, and write a timestep as At. For the dependent variables write, e.g. ul for u(i Ax, n At), where we interpret ul as the average of u in the cell ((i - 1) Ax, (i + i) Ax) at time n At. Assume we know U at time 12At and all grid points i = 1, . . . , I + 1; we describe how to advance U at all grid points to time (n + 1) At. The first step is to split the relaxation effects from the propagation effects. That is, we solve in turn

a,u=R(U) and

a,u+ a,(Au)=o.

E. B. PITMAN

330

and Y. NI

In fact, we use a Strang splitting to maintain second order accuracy in time. Thus we solve a,U = R(U) for a half-timestep, then solve d,U + d,(AU) = 0 for a full timestep, then solve d,U = R(U) for another half-step. To solve the relaxation ODE, we use a 4th order Runge-Kutta scheme. The propagation equation is solved by a Godunov method. The second order Godunov method consists of four steps: slope determination; characteristic tracing; local Riemann solver; conservative differencing. We discuss these steps in turn. In order to attain second order accuracy in space, we construct a piecewise linear interpolant function from the grid values r/l. First we must project U into the space of characteristic variables. Define P to be the matrix whose columns are the right eigenvectors r_, c,, and r+ of A. The characteristic variable is W = P-‘U. We would like to use the centered difference (Wj+, - Wi_,)/2 AX as the slope of the interpolant. However in order to prevent nonphysical oscillations, we limit this slope so that the interpolating function in ((i - 1) AX, (i + 4) AX) is always bounded by its neighboring cell-average values, WY__, and WY+,,. Thus, define structure coefficients PI = P;’ AYU, where A? denotes either the forward, backward, or centered differences operators (for q =f, b, c, respectively). Set

Here, sgn(z) gives the sign of the argument.

Now define the limited slopes to be

The second step in the Godunov method is to define cell-edge time-centered via characteristic tracing. To this end, Taylor expand U as

Now replace the time differentiation d,U by space differentiation Factoring and replacing d, by the limited slope, we can write

-A

predictor steps

d,U using the equations.

In this expansion, all waves are propagated to the cell edge at (i + f) AX. However only waves associated with the two non-negative characteristic values ought to be used for this tracing. We throw away waves associated with negative characteristic speeds by defining projection matrices Z* as E* = P diag[ i (1 f sgn(l))]P-i where the diagonal matrix diag[$(l f sgn(n)] has entries 3~10, depending on the sign of the corresponding eigenvalue. This corrects the cell-edge predictor, which we now write as

Similarly define

r/y-+,:::U~-~~,(I+A;~) R

=

Ail/.

Using equations (7-8) we have defined two time-centered values next step in the method is to obtain a unique half-time predictor applications to hyperbolic conservation laws, this step corresponds problem at (i + 1) Ax, with data fJ~~,‘,& and Us!,:&. For linear

of I!/ at every cell edge. The value at the cell edges. In to solving the local Riemann equations, waves superpose

Visco-elastic

and we can simply project compute

the difference

K =

relaxation

r/y,?&

P-‘( r/;:&

331

- Uyz&?, onto the eigenvectors.

Thus,

- r/;;&).

The stationary state is then either U’ = Ui + K’r_ + K2cB or u’ = Vi+, - K2ro- K3r+. Here, Kk is the kth component of the vector K. In fact, the zero speed wave splits the stationary state, so we average, defining (9) The final step in the method is a conservative

differencing to update the equations:

The Godunov scheme is stable under the CFL restriction

However a more restrictive timestep is required for the Runge-Kutta and to produce correct results. In practice, we use

integration

to be stable

r At = const. min @ E ‘sup IPI I with a const. -0.9. The method is formally second order accurate in space and time. We mention here the work of [4,9] which clearly shows that numerical schemes for conservation laws with stiff source terms may give incorrect results if all length and timescales are not adequately resolved. Simple calculations of the system (1,2,4) on course grids illustrate this point.

NUMERICAL

RESULTS

Here we present numerical results using the basic scheme, demonstrating convergence and the resolution of the Godunov method. We then report on calculations which illustrate interesting aspects of the relaxation system, and we compare these results with calculations for the system (1,2,3). In the first series of plots, we consider the initial data

The initial data for p is p =p”‘(u). For this data, the solution of (1,2,3) remains in the half-plane u > 1 and is therefore always hyperbolic. Figure 1 shows the solution at time t = 0.3. Unless otherwise noted, all calculations are performed on a grid of 400 points, and the relaxation time t = l/100. To illustrate convergence of the method, Fig. 2 shows the results of calculations for the same data and the same r, but with different grids. We computed on grids of 200,400 and 800 points, and plot the difference U800 - U200 and U800 - WOO, where, for example, U800 is the computed solution u on the 800 point grid. We have multiplied these differences by a factor of 10 to better illustrate the convergence. We remark that, in refining the space grid, the time step is automatically refined because of the CFL condition.

332

E. B. PITMAN

and Y. NI

U, V, P, Pref k0.3 2.0 -u --” --_ P _____.‘pre,

Fig. 1. Solution of the Riemann problem at time t = 0.3. Data is (u, u) = (1.1, -0.5) for x ~0 and (u, u) = (1.1,0.5) for x >O. For all our computations, data for p is p =p”‘(u). Unless otherwise noted, calculations are performed on a grid of 408 points, and the relaxation time r = l/100.

CONVERGENCE U8OOU200 and U8OOU400 U8OOU200 U8OOU400

-0.10

-0.20 L -1 .o

-0.5

0.0 X

0.5

1.0

Fig. 2. The difference U8CJO- U208 and UgoO - 1/400, times 10, where U800, U400 and U200 are the solutions u computed on a grid of 880, 480, and 280 points, respectively. Here the time is t = 0.3.

Visco-elastic relaxation

333

P-Pref r=1/5011100 l/200

X

Fig. 3. The difference p -pnr(u)

for t = l/50, l/100, l/200 at t = 0.3.

Figure

3 shows the effect of decreasing the relaxation time. In Fig. 3 we plot the difference for z = l/50, 1/100, l/200. Based on these and other similar computations, we note (i) the difference appears to decrease monotonically at each gridpoint, and (ii) the decrease seems to be linear in r. We now study the relaxation system for two other sets of initial data, and compare the results with a Lax-Friedrichs calculation for the system (1,2,3). The Lax-Friedrichs scheme for this system can be written

p -p”‘(u)

u;+’

=

u;+,+ u;_,

+z

2 v;+’ =

At (by+1 - tG-1)

v;+, + vi”_, At +z ((P”‘F+: 2

- (P’“‘Kd

where (pref)l denotes p”‘(u;). Consider the Riemann data

The analytic solution to the system (1,2,3) is always hyperbolic, but u jumps across the elliptic region. Any numerical solution computed with a method which includes viscosity will include points inside the elliptic region. Figure 4(a) shows the numerical solution of the mixed system and Fig. 4(b) shows the numerical solution of the visco-elastic system. The linear waves of the visco-elastic system are obviously smeared during the calculation, but the computed results are quite similar (modulo the spikes in the Lax-Friedrichs solution near x = 0). The results of Fig. 5 are rather different. In this example the data is (-1.1,

@, This

U)=I(-l.I,

-0.5)

XC0

0.5)

x >o.

is an example of the Riemann data for which Shearer (121 found two solutions to the system (1,2,3); see also [7]. One of his solutions consists of four waves, a rarefaction followed

E. B. PITMAN

334

and Y. NI

U, V, Pref (cl) 2.0

k0.3

r

1.0

5 a’ >

0.0

1

-1 .o

-2.0

c

--_-_---.-_-___

0,

1

-1 .o

________. ~,_______________

I_________ \

-0.5

I

-----_

I:‘__.___________

1

___. _ ____ _.._

\ I

0.0 X

0.5

1.0

U, V, P, Pref (b)

1=0.3

Fig. 4. Solution of the Riemann problem at time I = 0.3. Data is (u, u) = (-1.1, -0.0) for x < 0 and (u, u) = (1.1,O.O) for x > 0. Figure 4(a) shows the results of a Lax-Friedrichs computation on a 400 point grid; Fig. 4(b) shows the results using the relaxed system.

by a phase jump followed by a phase jump back followed by another rarefaction. This is the solution computed by the Lax-Friedrichs scheme [Fig. 5(a)]. The second solution Shearer found consists of a pair of rarefaction waves, without any change in phase. The visco-elastic system [Fig. 5(b)] appears to start toward this rarefaction solution, generating a pair of weak expansion waves, but the relaxed system develops phase jumps as u enters the elliptic region. (REMARK: We have performed a series of refinements for this data in order to verify the behavior of the solution. The results shown in Fig. 5(b) used 1000 grid points and the

Visco-elastic relaxation

335

U, V, Pref (a)

MI.3

2.0

I--+ -U --v _____.‘pref

I

I

1.0 -

L._> \ Lo_“______________________.__ _____________________________\ : ,___ I ,%_--_--_--__-_--_--_--_--. : b :‘, ’: ‘5.____I !:.____, I I 0.0 :--I h’ :

% Ii > 1

-_---_-_---_-____

1 i

\ L_-1.0

-2.0 -1.0

i !

,i

I 1 I

0.0 X

-0.5

1.0

0.5

U, V Pref lb)

1=0.3 1000 grid points

3.0 -U ---

v

__---.

pre,

-

2.0

1.0 .____________________________________~

--

_-------______________________________ ________-___________

?=1/1000 0.0 -_-_-_____________-_-~ -1 .o

-2.0 1 -1 .o

-0.5

0.0 X

0.5

1.0

Fig. 5. Solution of the Riemann problem at time 1= 0.3. Data is (u, v) = (-1.1, -0.5) for x < 0 and (u, u) = (-1.1,O.S) for x > 0. Figure 5(a) shows the results of a Lax-Friedrichs computation on a 400 point grid; Fig. 5(b) shows the results using the relaxed system. For this plot, the calculations for the relaxed system used 1000 grid points and a relaxation time r = l/1000. Note scale in Fig. 5(b).

l/1000.) For clarity, Fig. 5(b) shows only u, V, and pref; the solution for p is close to the displayed p”‘. In Fig. 6 we illustrate the slow propagation of this change-of-phase solution. For the 400 point grid with t = l/100, we show u at times t = 0.1, 0.3, 0.5, 0.7, 0.9. By t = 0.9,the leading weak elastic wave has propagated past the end of the grid. However the phase jump moves slowly through the domain. relaxation

time

z =

E. B. PITMAN and Y. NI

336

U AT SUCCESSIVE TIMES t =0.1,0.3,0.5,0.7, 0.9 I

-1.o

I

I -2.0

-1.0

-0.5

0.0

0.5

1.o

X Fig. 6. The computed u at successive times t = 0.1, 0.3, 0.5, 0.7, 0.9, with data as in Fig. 5. For this plot, calculations were done using a 400 point grid and t = l/100.

Finally we present the solution to a Cauchy problem for the visco-elastic system. The initial data is U(X, 0) = sin(nx)

V(X, 0) = 0

P(% 0) = pref(u(x, 0)).

In Figs 7(a) and 7(b) we show the solution at times 0.15 and 0.3. Notice how, at t = 0.3, u is about to dive back through the elliptic region in a pair of new phase jumps. Calculations suggest that successively more phase jumps develop in time, leading to finer scale structure in the solution.

DISCUSSION

Following Suliciu, we study a mixed type system of conservation laws by examining a related visco-elastic system. We have used a Godunov method to numerically solve the governing PDEs. We report here on convergence of the method and present some interesting computational results. In particular, we find For most data, the solution of the Riemann problem for the relaxed system is near the solution of the mixed system. For such data, it appears that p + prcf pointwise almost everywhere. Using special data for which the Riemann problem of the mixed system has (at least) two solutions, the solution of the relaxed system converges to a pair of weak rarefaction waves separated by a pair of phase jumps, and this solution appears to be different from both the solutions displayed by Shearer. For the Cauchy problem, the solution of the relaxed system is capable of developing successively finer-scale structure as time evolves. Further study is needed to understand the relation of solutions of the relaxed system to solutions of the mixed-type system under various proposed admissibility criteria. Based on our computations, we are currently examining conditions under which solutions of the relaxed system converge to solutions of the mixed system.

Visco-elastic relaxation

U, V, P, Pref t=o.15

(a) 2.0

#’I’ ,x-L-., 1’ ,/ x: / ‘i / /

-2.0 1 -1.o

\

1

.L))

1

U

: ’ : ’’ : \ :9 ,8

--V ,’ --P /# _____.pref

-___AI, -0.5

0.0 X

.*’

0.5

1.0

0.5

1.o

U, V, P, Pref kO.3

lb)

-1.0

337

-0.5

0.0 X

Fig. 7. Solution of the Cauchy problem, with data u(x, 0) = sin(nx), u(x, 0) = 0, p(x, 0) = p”‘(u(x, 0)). Figure 7(a) shows the solution at t = 0.15; Fig. 7(b) shows the solution at t = 0.3. We have truncated the scale in Fig. 7(b) for sake of comparison with Fig. 7(a); the v component is off-scale near x = -1, 0, 1. For this calculation we return to a 400 point grid with t = l/100. Acknowledgenrenf.r-Our understanding of the mixed system has benefited from conversations with Michael Shearer, Marshall Slemrod, and Victor Roytburd. Computations were performed on the Cray Y-MP at the Pittsburgh Supercomputing Center under grant DMS!XXlO26P.

REFERENCES (11 R. ABEYARATNE [2] R. ABEYARATNE ES 32:2-J

and J. KNOWLES, Arch. Rd. Mech. Anal. 114, 119-154 (1991). and J. KNOWLES, J. Appl. Me&. 55,491-492 (1988).

E. B. PITMAN and Y. NI

338

[3] M. AFFOUF and R. CAFLISCH, SIAM J. Appl. Math. 51,605-634 (1991). [4] P. COLELLA, A. MAJDA and V. ROYTBURD, SIAM J. Sci. Stat. Compur. 7, 1059-1080 (1986). [S] R. COURANT and D. HILBERT, Methods of Mathematical Physics. Wiley, New York (1974). [6] C. FACIU, Int. J. Engng Sci. 29, 1103-119 (1991). 17) H. HATTORI, The entropy-rate admissibility criterion and the double phase boundary problem. In Contemporary Mathematics (Edited by B. KEYFITZ and H. KRANZER), Vol. 60. AMS, Providence, RI (1987). [8] R. JAMES, Arch. Ratl. Mech. Anal. 66, 59-125 (1983). [9] R. LEVEQUE and H. YEE, Conservation laws with stiff source terms. NASA Preprint (1988). [lo] G.-Q. CHEN and T. P. LIU Convergence of hyperbolic conservation laws with relaxation effect. MSRI Preprint (1992). [11] M. SHEARER, J. Difi Eqns 46,426-443 (1982). 1121 M. SHEARER. Arch. Ratl. Mech. Anal. 93. 45-59 (1986). i13j M. SLEMROD, Arch. Ratl. Mech. Anal. Si, 301-315 (1983). [14] 1. SULICIU, Int. J. Engng Sci. 28, 829-841 (1990). [ 151 J. TRANGENSTEIN, “A comparison of finite difference and finite element methods for elastic-plastic problems”. In Shock Waves and Viscous Profiles (Edited by M. SHEARER). SIAM, Philadelphia, PA (1991). [16] J. TRANGENSTEIN and P. COLELLA, Commun. Pure. Appl. Marh. 44,41-100 (1991). [17] G. B. WHITHAM, Linear and Nonlinear Waves. Wiley, New York (1978). (Received

7 December

1992; accepted 9 February

1993)