AllenCahn1D

Derivation of the exact solution

The Allen Cahn equation admits a 1D solution that can be used to test the implementation. The energy functional in 1D is

\[F[\alpha]=\int_\Omega\left(\lambda W(\alpha)+\frac{1}{2}\kappa\left(\frac{\partial\alpha}{\partial x}\right)^2\right)dx\]

The Allen-Cahn equation is

\[\frac{\partial\alpha}{\partial t}=-L\frac{\delta F}{\delta\alpha}\]

where the variational derivative can be derived from the Euler-Lagrange equation resulting in

\[\frac{\partial\alpha}{\partial t}=-L\left(\lambda\frac{dW}{d\alpha}-\kappa\frac{\partial^2\alpha}{\partial x^2}\right).\]

At steady-state, (\(\frac{\partial\alpha}{\partial t}=0\)) the Allen-Cahn equation simplifies to

\[\lambda\frac{dW}{d\alpha}=\kappa\frac{d^2\alpha}{dx^2}.\]

This is a nonlinear ODE, but it can be reduced by multiplying both sides by \(\frac{d\alpha}{dx}\):

\[\lambda\frac{dW}{d\alpha}\frac{d\alpha}{dx}=\kappa\frac{d^2\alpha}{dx^2}\frac{d\alpha}{dx}\]

Using the chain rule, we can further simplify:

\[\lambda\frac{dW}{dx}=\frac{\kappa}{2}\frac{d}{dx}\left(\frac{d\alpha}{dx}\right)^2\]

When integrated, we obtain

\[ \begin{align}\begin{aligned}\int\lambda\frac{dW}{dx}dx&=\int\frac{\kappa}{2}\frac{d}{dx}\left(\frac{d\alpha}{dx}\right)^2dx\\\lambda W(\alpha)+C_1&=\frac{\kappa}{2}\left(\frac{d\alpha}{dx}\right)^2+C_2\end{aligned}\end{align} \]

Combining the constants of integration into a single constant \(C_0=C_2-C_1\):

\[\lambda W(\alpha)=\frac{\kappa}{2}\left(\frac{d\alpha}{dx}\right)^2+C_0\]

Solving for \(\frac{d\alpha}{dx}\):

\[\frac{d\alpha}{dx}=\pm\sqrt{\frac{2}{\kappa}(\lambda W(\alpha)-C_0)}\]

For a monotonic interface profile connecting two phases, we require \(C_0=0\) to ensure that \(\frac{d\alpha}{dx}=0\) at the minima of \(W(\alpha)\):

\[\frac{d\alpha}{dx}=\pm\sqrt{\frac{2\lambda}{\kappa}W(\alpha)}\]

As is often the convention, we let \(W\) represent a double-well chemical potential

\[W(\alpha)=\alpha^2(1-\alpha)^2\]

which has a convenient square root. After substituting, we now have a nonlinear, but trivially solved differential equation:

\[\frac{d\alpha}{dx}=\gamma\alpha(1-\alpha)\]

where \(\gamma=\pm\sqrt{\frac{2\lambda}{\kappa}}\). This differential equation can be solved by separation of variables or by substitution as shown here:

\[u=\alpha^{-1}\implies\alpha=u^{-1},\frac{d\alpha}{dx}=-u^{-2}\frac{du}{dx}\]

Substituting \(u\) for \(\alpha\) we find

\[\frac{du}{dx}+\gamma u=\gamma,\]

which has the solution:

\[u(x)=Ce^{-\gamma x}+1\]

where \(C\) is a constant. Reversing the substitution, we get the following solution for \(\alpha\):

\[\alpha(x)=\frac{1}{Ce^{-\gamma x}+1}\]

We desire a solution centered at \(x=0\) such that \(\alpha(-x)=1-\alpha(x)\). This condition is only satisfied when \(C=1\); so, the interfacial solution is:

\[\boxed{\alpha(x)=\frac{1}{e^{-\gamma x}+1}}\]

Determination of interface width \(\varepsilon\)

The slope at the center of the interface (\(x=0\)) is

\[\left.\frac{d\alpha}{dx}\right|_{x=0}=\left.\frac{\gamma e^{-\gamma x}}{(e^{-\gamma x} + 1)^2}\right|_{x=0}=\frac{\gamma}{4}\]

The function changes from approximately 0 to 1 across the interface. Using a linear approximation slope—i.e., slope \(\approx\) (change in \(\alpha\))/(change in \(x\))—the interface width \(\varepsilon\) is as follows:

\[\frac{\gamma}{4}\approx\frac{1}{\varepsilon}\implies\varepsilon=\frac{4}{\gamma}=\boxed{\sqrt{\frac{8\kappa}{\lambda}}}\]

Determination of interface energy density \(\sigma\)

The interface energy \(\sigma\) is given by substituting the solution for \(\alpha\) into the free energy functional:

\[ \begin{align}\begin{aligned}\sigma&=F\left[\frac{1}{e^{-\gamma x}+1}\right]\\&=\int_{-\infty}^\infty\left(\lambda W(\alpha)+\frac{1}{2}\kappa\left(\frac{d\alpha}{dx}\right)^2\right)dx\end{aligned}\end{align} \]

Using the chain rule, the derivative of \(\alpha\) is

\[\frac{d\alpha}{dx}=-\frac{1}{\left(e^{-\gamma x}+1\right)^2}\left(-\gamma e^{-\gamma x}\right)=\frac{\gamma e^{-\gamma x}}{\left(e^{-\gamma x}+1\right)^2}.\]

Conveniently, this can be written as

\[\frac{d\alpha}{dx}=\gamma\alpha(1-\alpha).\]

Squaring this derivative, we can make another convenient substitution:

\[\left(\frac{d\alpha}{dx}\right)^2=\gamma^2\alpha^2(1-\alpha)^2=\gamma^2W(\alpha)\]

Plugging this back into the integrand:

\[\sigma=\int_{-\infty}^\infty\left(\lambda W(\alpha)+\frac{1}{2}\kappa\gamma^2W(\alpha)\right)dx\]

Recalling that \(\gamma^2=\frac{2\lambda}{\kappa}\):

\[ \begin{align}\begin{aligned}\sigma&=2\lambda\int_{-\infty}^{\infty}W(\alpha)dx\\&=2\lambda\int_{-\infty}^\infty\alpha^2(1-\alpha)^2dx\end{aligned}\end{align} \]

We can trivially solve this integral with a simple substitution:

\[ \begin{align}\begin{aligned}u=\alpha&\implies\frac{du}{dx}=\frac{d\alpha}{dx}=\gamma\alpha(1-\alpha)\\&\implies dx=\frac{du}{\gamma\alpha(1-\alpha)}=\frac{du}{\gamma u(1-u)}\end{aligned}\end{align} \]

To adjust the bounds of integration, we note:

\[ \begin{align}\begin{aligned}x&:-\infty\rightarrow\infty\\u&:0\rightarrow1\end{aligned}\end{align} \]

Substituting these variables into the integrand:

\[ \begin{align}\begin{aligned}\sigma&=2\lambda\int_0^1u^2(1-u)^2\frac{du}{\gamma u(1-u)}\\&=\frac{2\lambda}{\gamma}\int_0^1u(1-u)du\\&=\frac{\lambda}{3\gamma}\end{aligned}\end{align} \]

Applying one final substitution, we find the following expression for the interface energy \(\sigma\):

\[\boxed{\sigma=\frac{1}{3}\sqrt{\frac{\lambda\kappa}{2}}}\]

Inverting the relationship

It may be more desirable to define behavior in terms of the interface energy \(\sigma\) and interface width \(\varepsilon\). From the previously derived equation for \(\sigma\), we can show:

\[ \begin{align}\begin{aligned}\lambda&=\frac{18\sigma^2}{\kappa}\\\kappa&=\frac{18\sigma^2}{\lambda}\end{aligned}\end{align} \]

From the definition of \(\gamma\), we can use the following substitutions:

\[ \begin{align}\begin{aligned}\kappa&=\frac{2\lambda}{\gamma^2}\\\lambda&=\frac{\kappa\gamma^2}{2}\end{aligned}\end{align} \]

These relationships then give us the following expressions:

\[ \begin{align}\begin{aligned}\lambda&=\frac{9\sigma^2\gamma^2}{\lambda}=3\sigma\gamma\\\kappa&=\frac{36\sigma^2}{\kappa\gamma^2}=\frac{6\sigma}{\gamma}\end{aligned}\end{align} \]

Recalling that \(\gamma=\frac{4}{\varepsilon}\), we can finalize the inversion:

\[ \begin{align}\begin{aligned}\lambda&=\frac{12\sigma}{\varepsilon}\\\kappa&=\frac{3\sigma\varepsilon}{2}\end{aligned}\end{align} \]

Writing the free energy in terms of these quantities gives:

\[F\left[\alpha\right]=\int_\Omega\left(\frac{12\sigma}{\varepsilon}W(\alpha)+\frac{3\sigma\varepsilon}{4}\left(\frac{d\alpha}{dx}\right)^2\right)dx\]

Dimensional analysis

Given that \(\sigma\) has units of \(\left[\mathrm{E}\ \mathrm{L}^{-2}\right]\) and \(\varepsilon\) has units of length \(\left[\mathrm{L}\right]\), it becomes clear from the inverted relationships that \(\lambda\) and \(\kappa\) must have the following units:

\[ \begin{align}\begin{aligned}\left[\lambda\right]&=\left[\mathrm{E}\ \mathrm{L}^{-3}\right]\\\left[\kappa\right]&=\left[\mathrm{E}\ \mathrm{L}^{-1}\right]\end{aligned}\end{align} \]

eps_0.50_sigma_0.05

	Two-dimensional
	Serial
	Validated using check script
	./bin/alamo-2d-g++ tests/AllenCahn1D/input ch.eps="0.50" ch.sigma="0.05"

eps_0.50_sigma_0.1

	Two-dimensional
	Serial
	Validated using check script
	./bin/alamo-2d-g++ tests/AllenCahn1D/input ch.eps="0.50" ch.sigma="0.1"

eps_0.05_sigma_0.05

	Two-dimensional
	Serial
	Validated using check script
	./bin/alamo-2d-g++ tests/AllenCahn1D/input ch.eps="0.05" ch.sigma="0.01"

Input file (../../tests/AllenCahn1D/input)

#@
#@ [eps_0.50_sigma_0.05]
#@ exe=alamo
#@ dim=2
#@ args = ch.eps=0.50
#@ args = ch.sigma=0.05
#@ check-file = eps_0.50_sigma_0.05
#@
#@ [eps_0.50_sigma_0.1]
#@ exe=alamo
#@ dim=2
#@ args = ch.eps=0.50
#@ args = ch.sigma=0.1
#@ check-file = eps_0.50_sigma_0.1
#@
#@ [eps_0.05_sigma_0.05]
#@ exe=alamo
#@ dim=2
#@ args = ch.eps=0.05
#@ args = ch.sigma=0.01
#@ check-file = eps_0.05_sigma_0.01
#@

alamo.program = allencahn

### OUTPUT ###

plot_file = tests/AllenCahn1D/output

### MESHING ###

# Specify system units
system.length = m
system.time   = s

amr.plot_int = 1000
amr.max_level = 3
amr.max_grid_size = 500000
amr.blocking_factor = 2
amr.grid_eff = 0.8
amr.n_cell = 64 4 0

geometry.prob_lo = -1.0_m  0.000_m 0.0_m # [ m ] 
geometry.prob_hi =  1.0_m  0.125_m 0.0_m # [ m ]

geometry.is_periodic = 0 0 0

timestep  = 1.0_ms
stop_time = 5.0_s

alpha.bc.constant.type.xlo = neumann
alpha.bc.constant.type.xhi = neumann
alpha.bc.constant.type.ylo = neumann
alpha.bc.constant.type.yhi = neumann

alpha.ic.type = expression
alpha.ic.expression.region0 = "x>0"

ch.L = 1.0_1/Pa/s
ch.eps = 0.100.0_m
ch.sigma = 0.01_J/m^2