Source-linked AI summary
Linear, Second order and Unconditionally Energy stable schemes for The Viscous Cahn-Hilliard Equation with hyperbolic relaxation using the invariant energy quadratization method
Xiaofeng Yang, Jia Zhao
TL;DR
The challenge is to discretize the hyperbolic and nonlinear terms at high temporal order while preserving energy stability; the authors report that no known VCH-HR schemes combine easy implementation, unconditional energy stability, and second-order accuracy. The paper uses invariant energy quadratization to transform the free energy and semi-explicitly discretize nonlinear terms in two second-order time-marching schemes. The two schemes are second order accurate, easy-to-implement linear systems, and unconditionally energy stable; 2D and 3D simulations demonstrate stability and accuracy.
Problem
The challenge is to discretize the hyperbolic and nonlinear terms at high temporal order while preserving energy stability; the authors report that no known VCH-HR schemes combine easy implementation, unconditional energy stability, and second-order accuracy.
Source: §1 · p. 2Method
The paper uses invariant energy quadratization to transform the free energy and semi-explicitly discretize nonlinear terms in two second-order time-marching schemes.
Source: §1 · pp. 2–3Results
The two schemes are second order accurate, easy-to-implement linear systems, and unconditionally energy stable; 2D and 3D simulations demonstrate stability and accuracy.
Source: §1 · p. 3Takeaways & Limitations
Each time step requires solving a symmetric positive definite linear system, which can be handled with fast matrix solvers such as CG or other Krylov methods.
Source: §1 · p. 3The model uses either periodic boundaries or no-flux boundary conditions for φ and μ.
Source: §2 · p. 4
Abstract
from arXiv · showhide
In this paper, we consider numerical approximations for the viscous Cahn-Hilliard equation with hyperbolic relaxation. This type of equations processes energy-dissipative structure. The main challenge in solving such a diffusive system numerically is how to develop high order temporal discretization for the hyperbolic and nonlinear terms, allowing large time-marching step, while preserving the energy stability, i.e. the energy dissipative structure at the time-discrete level. We resolve this issue by developing two second-order time-marching schemes using the recently developed "Invariant Energy Quadratization" approach where all nonlinear terms are discretized semi-explicitly. In each time step, one only needs to solve a symmetric positive definite (SPD) linear system. All the proposed schemes are rigorously proven to be unconditionally energy stable, and the second-order convergence in time has been verified by time step refinement tests numerically. Various 2D and 3D numerical simulations are presented to demonstrate the stability, accuracy and efficiency of the proposed schemes.
1. Introduction
The paper targets a gap in practical second-order, unconditionally energy-stable schemes for the viscous Cahn–Hilliard equation with hyperbolic relaxation. It uses IEQ to develop two linear schemes and reports stability and accuracy in 2D and 3D simulations.
1. Introduction
- The VCH-HR model adds viscosity and hyperbolic relaxation to the classical Cahn–Hilliard equation, with the latter associated with strongly non-equilibrium decomposition during rapid solidification.The introduction also notes theoretical study of VCH and VCH-HR, including well-posedness, sharp-interface limits, and global attractors.Source: §1 · p. 1
- Unconditional energy stability supports capturing long-time dynamics and addressing stiffness, but larger time steps still increase computational error.The introduction motivates second-order schemes as a way to seek larger time steps while maintaining desirable accuracy.Source: §1 · p. 2
- To address the lack of easy-to-implement, unconditionally energy-stable, second-order VCH-HR schemes, the paper develops an IEQ-based approach.The authors identify discretizing viscosity, hyperbolic inertia, and nonlinear double-well stiffness as numerical difficulties.Source: §1 · p. 2
- IEQ transforms the free-energy potential into quadratic form, enabling semi-explicit treatment of nonlinear terms and a linear system.The two schemes use Crank–Nicolson and Adam–Bashforth (BDF2) time discretizations.Source: §1 · pp. 2–3
- The authors report two second-order schemes that are unconditionally energy stable, require solving linear systems, and show stability and accuracy in 2D and 3D simulations.Their linear operators are symmetric positive definite, allowing use of fast matrix solvers such as CG or other Krylov methods.Source: §1 · p. 3
2. Model Equations
The model describes a binary alloy through a volume fraction and diffusion flux, with the viscous and hyperbolic-relaxation terms incorporated into its energy. Its continuous energy dissipation law states that total free energy decays over time.
2. Model Equations
- The model considers a binary alloy in a bounded domain in two or three dimensions, and defines the volume fraction φ and diffusion flux J through a balance law.The domain boundary is Lipschitz continuous.Source: §2 · p. 3
- The total free energy includes interfacial and nonlinear bulk-potential contributions, plus a term weighted by the relaxation parameter α.The text identifies ε as the interfacial-width parameter and allows double-well or Flory–Huggins bulk potentials.Source: §2 · p. 3
- The model permits periodic boundaries or no-flux conditions for φ and the chemical potential μ.Source: §2 · p. 4
- When α = β = 0, the system reduces to the standard Cahn–Hilliard system, which conserves local mass density.Source: §2 · p. 4
- The model’s total free energy decays in time according to its energy dissipation law.The stated dissipation rate includes terms involving the gradient of an auxiliary variable and the time derivative of the volume fraction.Source: §2 · pp. 4–5
3. Numerical Schemes
This section develops two second-order IEQ schemes for the transformed nonlinear system and examines their energy stability and convergence. Numerical tests cover relaxation effects, coalescence dynamics, and 3D phase separation under different conditions.
3. Numerical Schemes
- IEQ rewrites the bounded-below nonlinear potential in quadratic form, yielding an equivalent transformed system with the same continuous energy law as the original system.For logarithmic potentials, the section regularizes the bulk free energy with a C2 piecewise function and takes B = 1.Source: §3 · pp. 5–6
3.1. Crank-Nicolson Scheme.
- The Crank-Nicolson scheme is unconditionally energy stable, and its linear system has a unique solution with a symmetric positive definite operator.The theorem establishes existence and uniqueness in H1(Ω).Source: §3.1 · pp. 8–9
3.2. Adam-Bashforth Scheme.
- The Adam-Bashforth BDF2 scheme provides an alternative second-order method and is unconditionally energy stable.The paper also states that its linear system is well-posed with a symmetric positive definite operator.Source: §3.2 · pp. 10–11
4.1. Convergence test.
- Both CN2 and BDF2 show second-order convergence in time in the refinement test at t = 4.Errors were measured as differences between solutions at adjacent time-step refinements.Source: §4.1 · p. 13, Figure 1 · p. 14
4.2. The viscous and hyperbolic relaxation effects for the coalescence of two kissing bubbles.
- With α = β = 0 or α = 1, β = 0, the bubbles reach the same steady state, while β = 1, α = 0 slows coalescence.For α = 1, β = 0, the interface shows a sawtooth profile at t = 1 and 2 before becoming circular.Source: §4.2 · p. 14
- Changing β can dramatically affect the coalescence speed, while Scheme 1 and Scheme 2 predict the same dynamics across the tested α and β values.The tests use α = 0, 0.5, 1 and β = 0, 0.5, 1.Source: §4.2 · p. 15, Figure 5 · p. 15
4.3. Energy.
- The modified energy (3.11) is a second-order approximation of the original energy (2.14), and its energy decays with time steps as large as δt = 1.The comparison uses a fully implicit δt = 2^-12 solution as the benchmark; the time-step tests use δt = 0.00625, 0.01235, 0.25, 0.5, and 1.Source: §4.3 · p. 16, Figure 6 · p. 16, Figure 7 · p. 16
4.4. Spinodal Decomposition in 3D.
- With α = 0 and β = 0.9, φ0 = 0.5 forms a uniform two-layer structure by t = 1800, while φ0 = 0.3 yields cylindrical accumulation of MMs.The 3D simulations use random initial data; the figure captions show snapshots through t = 1800 and t = 2000, respectively.Source: §4.4 · p. 17, Figure 8 · p. 17, Figure 9 · p. 18
- The spherical segment under no-flux boundaries has lower energy than the sphere under periodic boundaries; both cases use φ0 = 0.3 and α = β = 0.The no-flux result is a 1/8 spherical segment at a corner, while the periodic result is a sphere with parts at each cube corner.Source: §4.4 · pp. 17–18
- With the Flory-Huggins potential, different initial values of φ0 produce dramatically different dynamics and distinct steady states.The paper leaves detailed correlations between initial values and final states out of scope.Source: §4.4 · p. 18
5. Concluding Remarks
The paper develops two second-order-in-time IEQ schemes for the viscous Cahn–Hilliard equation with hyperbolic relaxation. The schemes are energy stable and linear to implement, with conclusions also addressing logarithmic potentials and possible extensions.
5. Concluding Remarks
- The two IEQ schemes are second-order in time, unconditionally energy stable, and require only linear equations at each time step.Their induced linear system is symmetric positive definite, allowing Krylov methods with mass-lumping preconditioners.Source: §5 · p. 18
- For logarithmic potentials, the authors describe their schemes as the first linear and accurate schemes with provable energy stability.Source: §5 · p. 18
- The authors state that the approach can extend to gradient flows with complex free-energy nonlinearities and to consistent finite-dimensional Galerkin approximations.The proofs use a variational formulation with trial and test functions in the same space.Source: §5 · pp. 18–19
- The 3D spinodal-decomposition experiments include no-flux and periodic boundary conditions with φ0 = 0.3 and α = β = 0.Snapshots are shown at t = 1, 2, 5, and 50, alongside free-energy evolution plots.Source: Figure 10 · p. 19, Figure 11 · p. 19
- The Flory–Huggins 3D spinodal-decomposition simulations use α = β = 0.5 with initial φ0 values of 0.5 and 0.3.The figures show φ at t = 100, 200, 500, and 4000.Source: Figure 12 · p. 20