NUMERICAL STUDY OF PHASE TRANSITION IN VAN DER WAALS FLUID

. In this article, we use a relaxation scheme for conservation laws to study liquid-vapor phase transition modeled by the van der Waals equation, which introduces a small parameter (cid:15) and a new variable. We solve the relaxation system in Lagrangian coordinates for one dimension and solve the system in Eulerian coordinates for two dimension. A second order TVD Runge-Kutta splitting scheme is used in time discretization and upwind or MUSCL scheme is used in space discretization. The long time behavior of the ﬂuid is numerically investigated. If the initial data belongs to elliptic region, the solution converges to two Maxwell states. When the initial data lies in metastable region, the solution either remains in the same phase or converges to the Maxwell states depending to the initial perturbation. If the initial state is in the stable region, the solution remains in that region for all time.


1.
Introduction. Van der Waals equation in dealing with the behavior of phase transition is very important. The equation of state is following where p is pressure, V is volume and T is temperature. It is easy to see that p is a monotonically decreasing function of V for T > 1. For T < 1, there is an interval of V where p increases as V increases, see [8] and [18]. Phase transition takes place around the region where p(V ) is increasing. In this article, we are concerned with the equations of isothermal motion for a compressible fluid in Lagrangian coordinates V t − u x = 0, where u is the velocity. The eigenvalues of (2) are ± −p (V ), which become complex in the elliptic region, i.e., p (V ) > 0. For a fixed temperature T < 1, the system is unstable when α < V < β which corresponding to elliptic region. When 2. The relaxation systems. Consider the 1-D relaxation system where u ∈ R n , v ∈ R n , x ∈ R 1 , t > 0 and A = diag{a 1 , a 2 , ..., a n }, a i > 0 (1 ≤ i ≤ n). For sufficiently small, it is expected that by solving (3) properly, one can obtain good approximations to the original conservation laws In (3), is a small positive parameter called the relaxation time and A satisfies Liu's sub characteristic condition where F (u) is the Jacobian matrix of F . Letting u = (V, u) T and F (u) = (−u, p(V )) T , and introducing a new variable v ∈ R 2 , the relaxation system of van der Waals equation in Lagrangian coordinates is (3). One simple choice of the matrix A is i.e., they are chosen to be the maximal(real) eigenvalue of (2).
2.1. The relaxation schemes. A spatial discretization to (3) in conservation form can be written as where mesh size h j = x j+1/2 − x j−1/2 , x j+1/2 is the spatial grid points. We denote by u j+1/2 the approximate point value of u, and by u j the approximate cell average of u in the cell We shall consider the upwind scheme and the MUSCL scheme, which are described in [12].
The upwind scheme. The first order semi-discrete upwind approximation to the relaxation system (3) is following The MUSCL scheme. Instead of using the piecewise constant interpolation which yields a first order approximation, the MUSCL uses the piecewise linear interpolation, which is a second order scheme. The MUSCL for the relaxation system (3) is componentwise as where u and v are the p-th components of u and v respectively, F (p) is the p-th component of F , σ ± j is the slope of v ± √ a p u on the j-th cell, defined by Sweby's notation [19], The minmod slope φ(θ) = max (0, min(1, θ)) is used in this article. The Time Discretization. A second order TVD Runge-Kutta splitting scheme to (3) gives where for the upwind scheme for the MUSCL scheme, the p-th components of u j+1/2 and v j+1/2 are calculated by From (10), the algorithm will break down if ∆t = O( ). Therefore, we choose ∆t = O( ).
3.1. One dimensional case. The equations (3) are solved numerically with the van der Waals equation of state at The system is elliptic for α < V < β, where α = 0.7186 and β = 1.5285. The Maxwell line intersects the curve of p(V ) at V 1 = 0.6034 and V 2 = 2.3488. We use the upwind or MUSCL scheme in space and TVD Runge-Kutta algorithm (9)-(18) in time. The CFL number is defined as where a i is introduced in Section 2. The periodic boundary condition is employed in the computational domain (−π, π). Some cases the results of upwind scheme are more stable than that of MUSCL scheme. But the MUSCL scheme is more accurate than upwind scheme. Our numerical results in the following figures show that (α, β) is a unstable region, (V 1 , α) and (β, V 2 ) is metastable region, V > V 2 and V < V 1 corresponds to stable region. All the computations are performed in uniform mesh and after refinement, the results are convergent. In the following figures, dashed lines represent the Maxwell states and the elliptic region is between the dotted lines. We now present our numerical results.
(1) The initial data that lie entirely in the elliptic region with The is taken to be 0.0025. The initial data V 0 correspond to a uniform liquid-vapor mixing state, u 0 is a small perturbation of velocity from equilibrium. From the numerical results, we clearly see the oscillations are generated inside the elliptic region due to the instability. The solution gradually evolves into equilibrium states (Maxwell states) in the stable part of the liquid phase and vapor phase, respectively, and those two states are connected by phase boundaries. We plot V at various times in Figure 1. We can see that V stabilized after time T = 50. We observe a stable configuration of phase separation.
(2) The initial data involve a larger perturbation of velocity, which is We report different from 10 −8 to 0.05. The solution is stable after time T = 50(see Figure 2). Figure 3 shows that when is increased to 0.005, the oscillations are much weaker. After T = 50, the oscillations merge and the system reach the equilibrium states with two separate phases, i.e. Maxwell states. We increase to 0.025 and 0.05, the two equilibrium states are different the previous two cases, the phase interface is much smoother and less sharp, and we also see fewer oscillations at the early stage of the phase transition in Figure 4 and 5 respectively, because larger produces greater dissipation. After long time simulation, both upwind scheme and MUSCL scheme can obtain the same results for the above two cases. The results (1) and (2) verify that if the initial data are a constant state in elliptic region (α < 1 0.9 < β), the instability leads to phase transition and eventually phase separation.
(3) Next we consider the initial data is in metastable region, i.e. (V 1 , α) or (β, V 2 ). Our numerical results in Figure 6 shows there is no phase transition and the solution will remain in the hyperbolic region for (V 0 , u 0 ) T = ( 1 1.6 , 0.1e −4x 2 ) T , where V 1 < 1 1.6 < α. But if we take initial data as (V 0 , u 0 ) T = ( 1 1.6 , 0.5e −4x 2 ) T , with higher amplitude, some part of the solution goes into the elliptic region and consequently phase transition can be seen ,which is shown in Figure 7. In this case, the MUSCL scheme can resolve better result than upwind scheme since MUSCL is second order accurate in space.
(4) For non constant initial states, if average volume which defined by 23, we can see the phase transition and the solution tends to the Maxwell states in Figure 8.
(5) For a constant state below the Maxwell state , the solution is stable which shown in Figure 9. For a non constant state initially, long time behavior is entirely hyperbolic, which can be observed in Figure 10. For case (4) and case (5), the upwind scheme and MUSCL scheme obtain much similar numerical results.
All one dimensional results are shown with 1025 grids. Actually, we have checked the convergence of the relaxation scheme by using a sequence of refined meshes, i.e. grid points are 257, 513, 1025, 2049. With various initial data, we confirm that average mass V 0 of the system determines the stability behavior and the final equilibrium states, which match the results obtained in [11]. We can use much less grid points than [11] to resolve the long time behavior for different initial data.

3.2.
Two dimensional case. We solve two dimensional relaxation system (20) of van der Waals equation in Eulerian coordinates with p(ρ) = −3ρ 2 + 8T ρ 3−ρ and T = 0.9. The computational domain is (−π, π) × (−π, π). The CFL number is defined as where a i and b i are introduced in Section 2. The periodic boundary conditions are used in both x direction and y direction.
All two dimensional numerical results have similar behavior to the one dimensional case. The convergence of the 2-D relaxation scheme has also been checked by a sequence of refined meshes, i.e. grid points are 64 × 64, 128 × 128, 256 × 256, 512 × 512. For space discretization, both upwind scheme and MUSCL scheme work well for these 2-D examples.
4. Conclusions. We use relaxation scheme to study the dynamics of liquid-vapor phase transition for a van der Waals fluids. The relaxation system has a linear hyperbolic structure which enables us to avoid the time consuming Riemann solvers, proper implicit time discretizations can be taken to overcome the stability constraints brought in by the stiffness, and we do not need to solve nonlinear systems of algebraic equations. We solve the conservation laws in Lagrangian coordinates for one dimension and solve the system in Eulerian coordinates for two dimension. The large time behavior of the fluid is numerically investigated. If the initial data Figure 14. (ρ 0 , u 0 , v 0 ) T = (1.7, 0.01e −4(x 2 +y 2 ) , 0.01e −4(x 2 +y 2 ) ) T , = 0.002, 256 × 256 grid points.
belongs to elliptic region, the solution converges to two Maxwell states. When the initial data lies in metastable region, the solution either remains in the same phase or converges to the Maxwell states depending to the initial perturbation. If the initial state is in the stable region, the solution remains in that region for all time. Our results are agreed with the results in [11], which is in Eulerian coordinates and we use much less grid points in this article. In the future, we are going to use the relaxation scheme to study phase transition through Navier-Stokes-Cahn-Hilliard system with van der Waals equation of state.