Numerical Investigation of a Thermoelastic Model for a Propagating Phase Boundary

A brief discussion of a numerical investigation of admissibility criteria for one-dimensional propagating phase boundaries in thermoelastic solids with nonconvex energies via viscosity-capillarity regularization, undertaken under the supervision of Prof. Fosdick.

By the late 1990s the case of a single conservation law in one spatial dimension with a flux function f which was neither strictly convex nor strictly concave was becoming well-understood; see, for example, [21,22].In this case, the sound speed is given by f (u), was such that the chord connecting (u − , f (u − )) and (u + , f (u + )) cut across the graph of f (u), as in Fig. 1, then the admissible solution consists of, in this instance, a shock from u − to the tangent (sonic) state u t , followed by a smooth spreading wave (similarity solution) from u t to u + ; this is only possible if u − and u + are separated by an inflection point of f .Professor Fosdick was curious if such a sonic connection also occurred in the system arising from a thermoelastic solid with a nonconvex energy, and if this would lead to a selection criterion for admissible shocks by augmenting the thermoelastic constitutive equations with viscosity and capillarity which were then allowed to go to zero.In order to introduce capillarity into the model in a thermodynamically consistent way, the nonclassical thermodynamic framework of Dunn and Serrin [23,24] was utilized, which includes interstitial working; thus the total internal energy ε, the total stress T , the interstitial working , and the heat flux q depend on the absolute temperature θ > 0, the deformation gradient F > 0, the thermoelastic internal energy e(θ, F ), the thermoelastic stress σ (θ, F ), the reference mass density ρ 0 , and three constants μ, δ, and κ, the coefficients of viscosity, capillarity, and thermal conductivity, respectively, in the following way: Here X ∈ R is the material coordinate, "X" subscripts denote partial derivatives, and a " ˙" indicates a material time derivative.A model thermoelastic Helmholtz free energy ψ = e − θη, where η is the entropy, given by η = −∂ψ(θ, F )/∂θ , was used for extensional motions (F ≥ 1) of a thermoelastic solid: where a > 0, b < 0, c > 0, and d > 0 are constants; the total Helmholtz free energy is ψ + The constant c has the physical interpretation of being the specific heat at constant deformation gradient c v .
From (1), the thermoelastic stress is The spinodal boundary consists of those (θ, F ) for which ∂σ/∂F = 0; from (2): Below the critical temperature θ cr = 9b 2 /(20ad), it is straightforward to solve (3) for two real roots Deformation gradients F ∈ [1, α(θ)] are said to be in the (low strain) α-phase, while those in [β(θ), ∞) are said to be in the (high strain) β-phase.(Deformation gradients between α(θ) and β(θ) correspond to statically unstable states.)It is also possible to solve (3) for θ in terms of F , and then substitute this value into (2) to obtain an explicit representation of the spinodal boundary: The metastable, or global convexity of the free energy, boundary for this material model must be found numerically by determining two distinct deformation gradients for each temperature which have equal stresses (2) and Gibbs functions The isentropic sound speed ς is given by ρ 0 ς 2 = ∂σ (η, F )/∂F evaluated at η = η(θ, F ), or, equivalently, and ς = 0 determines the hyperbolic boundary (This is not plotted in Fig. 2, but approaches the spinodal boundary from below as the specific heat c → ∞.)With a = 1, b = −8, c = 10, and d = 5, the constitutive relation (1) leads to a specific model which is depicted in Fig. 2. For these choices, the critical temperature θ cr = 5.76.
Normalizing the (constant) reference mass density so that ρ 0 = 1, the governing equations for mass, momentum and energy in one-dimensional material coordinates are where v is the velocity.With viscosity and capillarity it is expected that the solution will remain smooth, so the jump conditions are not needed.(All references to discontinuities in the numerical solutions are to small regions of large gradients.)For investigation of the phase portraits referenced below, note that if (6a) is differentiated with respect to t , and (6b) is differentiated with respect to X, this system may be transformed into two equations in which v does not appear: The first task was to look for shock profiles, i.e., solutions of (6a)-(6c) depending only on the variable ξ = (X − Vt)/ √ δ which may be used to approximate discontinuous solutions in the limit as δ → 0. (The velocity V, the shock velocity, is determined below.)Seeking solutions that are functions of ξ and integrating, assuming that the limits lim ξ →±∞ (F, v, θ )(ξ ) = (F ± , v ± , θ ± ) exist due to the presence of viscosity, leads to the following system of ordinary differential equations. where The Hugoniot function based on the right-hand state is and H + = 0 represents those states possible to reach with a discontinuity.The phase portraits of (7b), converted to two first-order equations, and (7c) in R 3 connecting the end states (F ± = 0, F ± , θ ± ) give the shock profiles.One observation is that as the coefficients of viscosity, capillarity, and thermal conductivity tend to zero, the ratios m and k should be kept constant.This aspect of the work however is not reported on further in this article.
From the assumptions made it follows that so the slope of the Rayleigh line connecting (F + , σ + ) to (F − , σ − ) determines the possible shock velocities.In the classical case, V must satisfy admissibility criteria, i.e., that the shock is supersonic with respect to the material it is moving into, and subsonic with respect to the material behind it.This determines whether a shock or rarefaction is the stable solution.For our material it is possible to choose initial data on the Hugoniot curve H + = 0 such that the shock velocity computed by the Rayleigh line is subsonic with respect to the material on both sides.This is the case we wished to investigate.

Results
Equations (6a)-(6c) were used to devise a finite difference scheme, which was used to approximate the solutions given arbitrary initial data.Two examples of the results for (smoothed) Riemann initial data with the right-side state (F + , θ + ) in the globally stable portion of the α-phase and the left-side state (F − , θ − ) in the globally stable portion of the β-phase follow.Specifically, For the presented simulations, the following parameters were used: μ = 10 −3 , δ = 10 −8 , and κ = 10 −3 .Using (8), the slope of the Rayleigh line R corresponds to a shock velocity of |V| = 0.605, which is readily seen to be smaller than either of the sound speeds in (9a)-(9b).Simulations with velocity jump v + − v − = ±0.457are illustrated in Figs. 3 and 4. The upper plots in Figs. 3 and 4 show the computed solutions (F black, θ gray) as a function of X at a time large enough to discern the structure (t = 4).More illuminating are the plots for these solutions in the F -σ plane in the lower graphs, which include the isotherm θ = 4, the Hugoniot curve H + = 0, as well as the metastable, spinodal, and hyperbolic boundaries.
Figure 3 (for v + − v − = +0.457)represents the case where the material X > 0 is moving to the right, and the material X < 0 is moving to the left.In this case it is expected to see the creation of additional high strain β-phase material.The calculated intermediate states are estimated as:  4 (for v + − v − = −0.457)represents the case where the material X > 0 is moving to the left, and the material X < 0 is moving to the right.In this case it is expected to see the creation of additional low strain α-phase material.The calculated intermediate states estimated as: The slope of the Rayleigh line R in Fig. 4, computed using (8) with the calculated front and back states, is 0.1287 corresponding to a phase boundary velocity of ±0.3587; this is to be compared with the calculated velocity of −0.3644.Here the sound speed in state 1 is 0.3390, a much closer match to the phase boundary velocity.(The discontinuity connecting the + state to the state 3 is traveling at a speed V = 1.7343, which is again intermediate to the sound speeds on either side.)This investigation was more thoroughly documented in [25].
Fig. 4 The solution of (6a)-(6c) with initial Riemann data (9a)-(9b) and v + − v − = −0.457.In this case it is expected to see the creation of additional low strain α-phase material.The Rayleigh line R connecting the + and − states does not represent the discontinuity in the calculated solution, but rather there is a spreading wave from the − state to state 1, near the tangent point, the spinodal boundary and the hyperbolic boundary, followed by a slowly propagating phase boundary to state 2 in the α-phase and shown as the computed Rayleigh line R , a contact discontinuity at X = 0 to state 3, and finally a (smoothed) classical shock connecting back to the initial data in the + state

Discussion
A large variety of simulations, with viscosity, capillarity and thermal conductivity at increasing smaller values were tantalizing, but did not unambiguously lead to the conclusion that there was always a sonic connection for the stable numerical solution.The spinodal and hyperbolic boundaries could also have been implicated.We concluded that an analytical approach was warranted.

Fig. 1
Fig.1If the states u − and u + are such that the chord connecting them cuts across the graph of f (u), the admissible solution jumps from the state u − to the state u t , where the chord is tangent to the flux function, followed by a smooth similarity solution to u +

Fig. 2
Fig.2Several isotherms of the thermoelastic stress (2), along with the spinodal boundary (4) (dotted), and the metastable boundary (dashed) (8) slope of the Rayleigh line R in Fig.3, computed using(8)with the calculated front and back states, is 0.0583, which corresponds to a phase boundary velocity of ±0.2415; this is to be compared with the calculated velocity of +0.2349.The crucial comparison is with the sound speed of state 1: 0.581.Thus while the plots are intriguing, the quantitative agreement is lacking.(The discontinuity connecting the − state to the state 3 is traveling at a speed |V| = 2.279, which is intermediate to the sound speeds on either side.)