有限元求解相场模型，Cahn—Hilliard 方程（非线性）

积分酱

已于 2024-11-18 06:06:52 修改

阅读量1.7k

点赞数 24

分类专栏：有限元文章标签：算法

于 2024-07-15 03:47:27 首次发布

本文链接：https://blog.csdn.net/wv0uue68/article/details/140427099

版权

有限元专栏收录该内容

3 篇文章

订阅专栏

目标方程

The given system of equations is the Cahn-Hilliard equation：

$\frac{\partial \varphi }{\partial t} - M_c\Delta w = f_1 \\ w + r\epsilon \Delta \varphi - r f(\varphi) = f_2$

where $\varphi(x,y,t)$ and $w (x, y, t)$ are the unknown functions, and $f_1(x,y,t)$ , $f_2(x,y,t)$ , $M_c$ , $r$ , and $\epsilon$ are known parameters.

Where $\varphi(x,y,t)$ and $w (x, y, t)$ are unknown functions. $f_1(x,y,t)$ , $f_2(x,y,t)$ , $M_c$ , $r$ , $\epsilon$ are known.

Here $f(\varphi) =F'(\varphi)= \frac{1}{\epsilon}\varphi(\varphi^2-1)$ , $F(\varphi)$ is Double-well potential function.

The space domain is $\Omega = [0,1] \times [0,1]$ , and the time domain is $T = [0, 1]$

$\nabla w = (w_x,w_y) \\ \Delta w =\nabla \cdot \nabla w = w_{xx} + w_{yy}\\ \nabla \varphi = (\varphi_x,\varphi_y)\\ \Delta \varphi =\nabla \cdot \nabla \varphi = \varphi_{xx} + \varphi_{yy}$

凸分裂（Convex splitting）

For example, Cahn-Hilliard equation is :
$u_t = \epsilon ^2\Delta^2u + \Delta f(u)$
$f = F^{'}$ , and $F=\frac{1}{4} (u^2-1)^2$ is a Double-well potential function.\
It also can be writen as energy functional form like:
$\int_\Omega \left(\frac{\epsilon^2}{2}|\nabla u|^2 + F(u)\right)d\Omega$

The convex splitting method is one of the commonly used approaches to construct unconditionally uniquely solvable and energy-stable numerical schemes. It was initially proposed by Eyre [2] in the study of unconditionally energy-stable schemes for the Cahn-Hilliard equation. Subsequently, this method has been extended and successfully applied to various other phase field models. \

The basic idea of the convex splitting method is to decompose the energy functional into the sum of a convex part and a concave part, treat the terms corresponding to the convex part implicitly, and the terms corresponding to the concave part explicitly. The resulting semi-implicit scheme is referred to as the convex splitting scheme.

The energy functional corresponding to the Cahn-Hilliard equation has the following decomposition:\
$\int_\Omega \left(\frac{\epsilon^2}{2}|\nabla u|^2 + F_{+}(u)\right)d\Omega - \int_\Omega \left(\frac{\epsilon^2}{2} F_{-}(u)\right)d\Omega$

where $F_{+}(u) = \frac{1}{4}(u^4+1)$ and $F_{-}(u) = \frac{1}{2}u^2$ satisfy $F(u) = F_{+}(u) - F_{-}(u)$ .

The first-order convex splitting scheme for the Cahn-Hilliard equation can be written as follows:

$\frac{u^{n+1} - u^n}{\Delta t} = \Delta w^{n+1}$

$w^{n+1} = -\epsilon^2\Delta u^{n+1} + \left( f_{+}(u^{n+1}) - f_{-}(u^n)\right)$
where $f_{+}(u) = F'_{+}(u) = u^3$ , $f_{-}(u) = F'_{-}(u) = u$ , and $u^n$ is the numerical approximation of $u(t_n)$ .

弱格式

第一个方程

$\frac{\partial \varphi }{\partial t} - M_c\Delta w = f_1$

multiply a function a test function $v (x, y)$ on both side of the the above equation.
$\frac{\partial \varphi }{\partial t}v - (M_c\Delta w) v = f_1v$
$\int_\Omega\frac{\partial \varphi }{\partial t}vdxdy - \int_\Omega(M_c\Delta w) vdxdy = \int_\Omega f_1vdxdy$

$\varphi(x,y,t)$ and $w (x, y, t)$ are called trail function and $v (x, y)$ is called a test function.
using Green’s formula (divergence theory, integration by parts in multi-dimension).
$\int_\Omega(M_c\Delta w) vdxdy =\int_\Omega(M_c\nabla \cdot \nabla w) vdxdy= \int_{\partial\Omega} (M_c \nabla w \cdot \vec n)v - \int_\Omega M_c\nabla w \cdot \nabla vdxdy$
Assume it is Dirichlet boundary condtion. Then $v = 0$ on $\partial \Omega$
$\int_\Omega\frac{\partial \varphi }{\partial t}vdxdy +\int_\Omega M_c\nabla w \cdot \nabla vdxdy = \int_\Omega f_1vdxdy$
$\nabla w \cdot \nabla v = w_xv_x+w_yv_y$
$\int_\Omega\frac{\partial \varphi }{\partial t}vdxdy +\int_\Omega M_c w_xv_xdxdy + \int_\Omega M_c w_yv_ydxdy = \int_\Omega f_1vdxdy$

Let $(\varphi_t,v) = \int_\Omega\frac{\partial \varphi }{\partial t}vdxdy$
$=\int_\Omega M_c\nabla w \cdot \nabla vdxdy$
$\int_\Omega f_1vdxdy$ \

$H^1(0, T; H^1(\Omega)) = \left\{ v(t, \cdot) \, \middle| \, \frac{\partial v}{\partial t}(t, \cdot) \in H^1(\Omega), \, \forall t \in [0, T] \right\}$
Weak formulation: find $w\in H^1(0, T; H^1(\Omega))$ and $\varphi \in H^1(0, T; H^1(\Omega))$ , such that \
$(\varphi_t,v) + a1(w,v) = (f1,v)$
for any $v\in H^1_0(\Omega)$

第二个方程

在这里插入图片描述

半离散

第一个方程

Assume there is a finite dimensional subspace $U_h = \text{span}\{\phi_j\}_{j=1}^{N_b}$ ,
where $\{\phi_j\}_{j=1}^{N_b}$ are the global finite element basis functions.
Then the Galerkin formulation is to find $\varphi_{h} \in H^1(0,T;U_h)$ and $w_{h} \in H^1(0,T;U_h)$
such that
$(\varphi_{h_t},v_h) + a1(w_h,v_h) = (f1,v_h)$
for any $v_h \in U_h$
It also can be express as:
$\int_\Omega\varphi_{h_t}v_hdxdy +\int_\Omega M_c\nabla w_h \cdot \nabla v_hdxdy = \int_\Omega f_1v_hdxdy$
or
$\int_\Omega\varphi_{h_t}v_hdxdy+\int_\Omega M_c (w_h)_x(v_h)_xdxdy + \int_\Omega M_c(w_h)_y(v_h)_ydxdy = \int_\Omega f_1v_hdxdy$
Since $\varphi_h \in H^1(0,T;U_h)$ and $w_{h} \in H^1(0,T;U_h)$ , and $U_h = \text{span}\{\phi_j\}_{j=1}^{N_b}$ , then
$\varphi_{h}(x,y,t) = \sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j(x,y)$

$w_{h}(x,y,t) = \sum_{j=1}^{N_b} w_{j}(t)\phi_j(x,y)$
for some coefficients $\varphi_{j}(t)(j=1,...,N_b)$ and $w_{j}(t)(j=1,...,N_b)$ .\
We choose $V_h = \phi_i(x,y)(i=1,...,N_b)$ .\
Then we can have
$\int_\Omega\left( \sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j\right)_t\phi_idxdy +\int_\Omega M_c\nabla \left(\sum_{j=1}^{N_b} w_{j}(t)\phi_j\right) \cdot \nabla \phi_i dxdy = \int_\Omega f_1\phi_idxdy$

$\int_\Omega( \sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j)_t\phi_idxdy +\int_\Omega M_c \sum_{j=1}^{N_b} w_{j}(t)\nabla \phi_j \cdot \nabla \phi_i dxdy = \int_\Omega f_1\phi_idxdy$

$\sum_{j=1}^{N_b}\varphi'_{j}(t) \int_\Omega\phi_j\phi_idxdy +\sum_{j=1}^{N_b} w_{j}(t)\int_\Omega M_c \nabla \phi_j \cdot \nabla \phi_i dxdy = \int_\Omega f_1\phi_idxdy$

$\sum_{j=1}^{N_b}\varphi'_{j}(t) \int_\Omega\phi_j\phi_idxdy +\sum_{j=1}^{N_b} w_{j}(t)\left(\int_\Omega M_c \frac{\partial\phi_j}{\partial x}\frac{\partial\phi_i}{\partial x}dxdy +\int_\Omega M_c \frac{\partial\phi_j}{\partial y}\frac{\partial\phi_i}{\partial y}dxdy\right) = \int_\Omega f_1\phi_idxdy$

Define the stiffness matrix
$A_1= \left[\int_\Omega M_c \nabla \phi_j \cdot \nabla \phi_i dxdy\right]_{i,j=1}^{N_b,N_b}$
$A_1= A_{11} + A_{12}$
$A_{11} = \left[\int_\Omega M_c \frac{\partial\phi_j}{\partial x}\frac{\partial\phi_i}{\partial x}dxdy\right]_{i,j=1}^{N_b,N_b}$
$A_{12} = \left[\int_\Omega M_c \frac{\partial\phi_j}{\partial y}\frac{\partial\phi_i}{\partial y}dxdy\right]_{i,j=1}^{N_b,N_b}$
Define the mass matrix
$\left[\int_\Omega\phi_j\phi_idxdy\right]_{i,j=1}^{N_b,N_b}$

Define the load vector
$\vec b_1 = \left[\int_\Omega f_1\phi_idxdy \right]_{i=1}^{N_b}$
Define the unknown vector
$\vec X_1=\left[\varphi_j(t) \right]_{j=1}^{N_b}$
$\vec X_2=\left[w_j(t) \right]_{j=1}^{N_b}$

Then we obtain the system
$M\vec X_1'(t) + A_1\vec X_2 = \vec b_1$

第二个方程

在这里插入图片描述

全离散

第一个方程

$M\vec X_1'(t) + A_1\vec X_2 = \vec b_1$
$\vec X_1'(t) = \frac{b_1 - A_1\vec X_2}{M}$

Assume that we have a uniform partition of $[0, T]$ into $M_m$ elements with time step $\Delta t$ .
The mesh nodes are $t_m = m \Delta t, m = 0,1,...,M_m$ .
Assume $\vec X_1^{m}$ is the numerical solution of $\vec X_1(t_m)$ and $\vec X_2^{m}$ is the numerical solution of $\vec X_2(t_m)$ .

$\theta$ -scheme (Finite difference (FD) method):
$\frac{X_1^{m+1} - X_1^{m}}{\Delta t} = \theta\frac{b_1^{m+1} - A_1^{m+1}X_2^{m+1}}{M} + (1-\theta)\frac{b_1^{m} - A_1^{m}X_2^{m}}{M}, m = 0,1,...,M_m$
where:
$\theta = 0$ : forward Euler scheme;
$\theta = 1$ : backward Euler scheme;
$\theta =1/2$ : Crank-Nicolson scheme.

第二个方程

在这里插入图片描述

牛顿迭代

在这里插入图片描述

牛顿线性化

在这里插入图片描述

结果

在这里插入图片描述

syms x y t;

% parameters
Mc=1;
r=1;
eplson=1;


u= exp(x+y+t);
w=x*y*(1-x/2)*(1-y)*exp(x+y+t);



ut = diff(u,t);
% compute the Laplacian op
lap_w = laplacian(w, [x, y]);

f1 = ut-Mc*lap_w;

Lap_u=laplacian(u,[x,y]);
fu=1/eplson * u *(u^2-1);
f2 = w+r*eplson*Lap_u-r*fu;

% function f1 f2
fun_f1 = matlabFunction(f1, 'Vars', [x, y,t])

fun_f2 = matlabFunction(f2, 'Vars', [x, y,t])



wx=diff(w,x);
wy=diff(w,y);
% function of derivative 
fun_w_analytic_der_x = matlabFunction(wx, 'Vars', [x, y,t])
fun_w_analytic_der_y = matlabFunction(wy, 'Vars', [x, y,t])

ux=diff(u,x);
uy=diff(u,y);