有限元求解Cahn Hilliard 相场方程

积分酱

已于 2024-07-01 04:14:03 修改

阅读量510

点赞数 7

分类专栏：有限元文章标签：算法笔记数学建模

于 2024-06-23 06:37:43 首次发布

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

版权

有限元专栏收录该内容

7 篇文章 1 订阅

订阅专栏

Cahn Hilliard 方程, 相场

因为最近在学习用有限元求解Cahn Hilliard 方程内容，所以我把我的汇报内容放在这里，希望可以和大家一起学习，交流。因为我实在是懒得翻译，就放的英文版，但是英语都不多，我的英语水平也有限。如果有什么问题的话，也希望大家可以提供给我建议。

目标方程

$\frac{\partial \varphi }{\partial t} - M_c\Delta w = f_1\\ w + r\epsilon \Delta \varphi = f_2$
这里 $\varphi(x,y,t),w(x,y,t)$ 是需要求解的未知函数， $f_1(x,y,t),f_2(x,y,t)$ 是已知函数， $M_c,r,\epsilon$ 是参数，为了简单，我全部设置为1. 求解区域为 $\Omega = [0,1] \times [0,1]$ ,时间区域为 $T_d[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}$

弱格式

第一个方程

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

Multiply a function a test function
multiply 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.
Green’s formula
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)$

第二个方程

$r\epsilon \Delta \varphi = f_2$

multiply a test function $v$
multiply a function a test function $v (x, y)$ on both side of the the above equation.
$r\epsilon \Delta \varphi = f_2$
$\int_\Omega wvdxdy + \int_\Omega r\epsilon \Delta \varphi vdxdy = \int_\Omega f_2vdxdy$
$\varphi(x,y,t)$ and $w (x, y, t)$ are called trail function and $v (x, y)$ is called a test function.
Green’s formula
using Green’s formula (divergence theory, integration by parts in multi-dimension).
$\int_\Omega r\epsilon \Delta \varphi vdxdy =\int_\Omega( r\epsilon \nabla \cdot \nabla \varphi) vdxdy= \int_{\partial\Omega} ( r\epsilon \nabla \varphi \cdot \vec n)v - \int_\Omega r\epsilon \nabla \varphi \cdot \nabla vdxdy$

Assume it is Dirichlet boundary condtion. Then $v = 0$ on $\partial \Omega$
$\int_\Omega wvdxdy- \int_\Omega r\epsilon \nabla \varphi \cdot \nabla vdxdy = \int_\Omega f_2vdxdy$
$\nabla \varphi \cdot \nabla v = \varphi_xv_x+\varphi_yv_y$

$\int_\Omega w vdxdy -\int_\Omega r\epsilon \varphi_xv_x dxdy - \int_\Omega r\epsilon\varphi_yv_ydxdy = \int_\Omega f_2vdxdy$
Let
$\int_\Omega w vdxdy$
$a3(\varphi,v) = - \int_\Omega r\epsilon \nabla \varphi \cdot \nabla vdxdy$

$\int_\Omega f_2vdxdy$ \

Weak formulation: find $w\in H^1(0, T; H^1(\Omega))$ and $\varphi \in H^1(0, T; H^1(\Omega))$ , such that \
$a3(\varphi,v) = (f2,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$

第二个方程

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
$a3(\varphi_h,v_h) = (f2,v_h)$
for any $v_h \in U_h$

It also can be express as:
$\int_\Omega w_hv_hdxdy- \int_\Omega r\epsilon \nabla \varphi_h \cdot \nabla v_hdxdy = \int_\Omega f_2v_hdxdy$
or
$\int_\Omega w_h v_hdxdy -\int_\Omega r\epsilon (\varphi_h)_x (vh)_x dxdy - \int_\Omega r\epsilon(\varphi_h)_y (vh)_ydxdy = \int_\Omega f_2v_hdxdy$

Since
$\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)$ .\

$\int_\Omega \left( \sum_{j=1}^{N_b} w_{j}(t)\phi_j\right)\phi_idxdy- \int_\Omega r\epsilon \nabla \left(\sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j\right) \cdot \nabla \phi_idxdy = \int_\Omega f_2\phi_idxdy$

$\sum_{j=1}^{N_b} w_{j}(t)\int_\Omega \phi_j \phi_idxdy + \sum_{j=1}^{N_b} \varphi_{j}(t)\int_\Omega - r\epsilon \nabla \phi_j\cdot \nabla \phi_idxdy = \int_\Omega f_2\phi_idxdy$

$\sum_{j=1}^{N_b} w_{j}(t)\int_\Omega \phi_j \phi_idxdy + \sum_{j=1}^{N_b}\varphi_{j}(t) \left (\int_\Omega - r\epsilon \frac{\partial \phi_j}{\partial x} \frac{\partial \phi_i}{\partial x}dxdy + \int_\Omega - r\epsilon \frac{\partial \phi_j}{\partial y} \frac{\partial \phi_i}{\partial y}dxdy \right) = \int_\Omega f_2\phi_idxdy$

Define the stiffness matrix
$A_2= \left[\int_\Omega - r\epsilon \nabla \phi_j\cdot \nabla \phi_idxdy\right]_{i,j=1}^{N_b,N_b}$
$A_2= A_{21} + A_{22}$
$A_{21} = \left[\int_\Omega - r\epsilon \frac{\partial \phi_j}{\partial x} \frac{\partial \phi_i}{\partial x}dxdy\right]_{i,j=1}^{N_b,N_b}$
$A_{22} = \left[\int_\Omega - r\epsilon \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_2 = \left[\int_\Omega f_2\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
$A_2\vec X_1+ M\vec X_2=\vec b_2$

全离散

第一个方程

$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, 1]$ into $M_m$ elements with mesh size $\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)$
Assume $\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

Then we can have the ireration scheme:

$\frac{M}{\Delta t}\vec X_1^{m+1}+\theta A_1^{m+1}\vec X_2^{m+1} = \frac{M}{\Delta t}\vec X_1^{m}+\theta \vec b_1^{m+1} + (1-\theta)(\vec b_1^{m} - A_1^{m}\vec X_2^{m})$

第二个方程

$A_2^{m+1}\vec X_1^{m+1}+ M\vec X_2^{m+1}=\vec b_2^{m+1}$

So the systems becomes to

$\frac{M}{\Delta t}\vec X_1^{m+1}+\theta A_1^{m+1}\vec X_2^{m+1} = \frac{M}{\Delta t}\vec X_1^{m}+\theta \vec b_1^{m+1} + (1-\theta)(\vec b_1^{m} - A_1^{m}\vec X_2^{m}) \\ A_2^{m+1}\vec X_1^{m+1}+ M\vec X_2^{m+1}=\vec b_2^{m+1}$

$A^{m+1} = \begin{bmatrix} \displaystyle \frac{M}{\Delta t} & \displaystyle \theta A_1^{m+1} \\[0.4cm] A_2^{m+1} & M \end{bmatrix}$

$\vec{X}^{m+1} = \begin{bmatrix}\displaystyle \vec{X}_1^{m+1}\\ \vec{X}_2^{m+1}\end{bmatrix}$

$\vec{b}^{m+1} = \begin{bmatrix} \displaystyle \frac{M}{\Delta t}\vec{X}_1^{m} + \theta \vec{b}_1^{m+1} + (1-\theta)(\vec{b}_1^{m} - A_1^{m}\vec{X}_2^{m}) \\[0.4cm] \vec{b}_2^{m+1} \end{bmatrix}$

Then the system is
$A^{m+1} \vec X^{m+1} = \vec b^{m+1}$
We solve the equation. Then we can get the solution $\vec X^{m+1}$ .

结果

The exact solution: $\varphi = e^{x+y+t}$ , $w=xy(1-\frac{x}{2})(1-y)e^{x+y+t}$
Set $M c = 1$ , $r = 1$ , $\epsilon=1$ , $\theta = 1$ , then compute $f_1$ and $f_2$ using Matlab

%% Matlab
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 of w
lap_w = laplacian(w, [x, y]);

f1 = ut-Mc*lap_w;

% compute the Laplacian of u
Lap_u=laplacian(u,[x,y]);
f2 = w+r*eplson*Lap_u;


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

fun_f1 =
    @(x,y,t)exp(t+x+y)-x.*exp(t+x+y).*(x./2.0-1.0).*2.0-y.*exp(t+x+y).*(y-1.0)-x.*exp(t+x+y).*(x./2.0-1.0).*(y-1.0).*2.0-y.*exp(t+x+y).*(x./2.0-1.0).*(y-1.0).*2.0-x.*y.*exp(t+x+y).*(y-1.0)-x.*y.*exp(t+x+y).*(x./2.0-1.0).*2.0-x.*y.*exp(t+x+y).*(x./2.0-1.0).*(y-1.0).*2.0

fun_f2 =:
    @(x,y,t)exp(t+x+y).*2.0+x.*y.*exp(t+x+y).*(x./2.0-1.0).*(y-1.0)

Table1：Results for $\varphi$ ，Using linear basis function， Backforward Euler scheme（ $\theta =1$ ）
$\begin{array}{ccccccccc} \hline h_1 & h_2 & dt&L_\infty \text{ error} & \text{Convergence } L_\infty & L_2 \text{ error} & \text{Convergence } L_2 & H^1 \text{ semi-error} & \text{Convergence } H^1 \\ \hline 1/4 & 1/4 & 1/8 & 1.3626e-01 & & 4.9320e-02 & & 8.8899e-01 & \\ 1/8 & 1/8 & 1/64 & 3.6311e-02 & 1.9079 & 1.2185e-02 & 2.0171 & 4.4348e-01 & 1.0033 \\ 1/16 & 1/16 & 1/512 & 9.3747e-03 & 1.9536 & 3.0301e-03 & 2.0077 & 2.2161e-01 & 1.0008 \\ 1/32 & 1/32 & 1/4096 & 2.3818e-03 & 1.9767 & 7.5579e-04 & 2.0033 & 1.1079e-01 & 1.0002 \\ \hline \end{array}$

Table2：Results for $w$ ，Using linear basis function， Backforward Euler scheme（ $\theta =1$ ）
$\begin{array}{ccccccccc} \hline h_1 & h_2 & dt&L_\infty \text{ error} & \text{Convergence } L_\infty & L_2 \text{ error} & \text{Convergence } L_2 & H^1 \text{ semi-error} & \text{Convergence } H^1 \\ \hline 1/4& 1/4& 1/8 & 2.2236e-01& & 4.4761e-02& &7.6775e-01& \\ 1/8& 1/8 & 1/64 & 6.4205e-02 &1.7922& 1.2208e-02&1.8771 & 3.8709e-01&9.8796e-01\\ 1/16&1/16 & 1/512 & 1.7240e-02 & 1.8969 &3.2441e-03 &2.0104 &1.9411e-01&9.9579e-01 \\ 1/32&1/32 & 1/4096& 4.4660e-03 &1.9487& 8.3737e-04 &2.1018 & 9.7133e-02&9.9886e-01\\ \hline \end{array}$

Table3：：Results for $\varphi$ , Using Quadratic basis function， Backforward Euler scheme（ $\theta =1$ ）
$\begin{array}{ccccccccc} \hline h_1 & h_2 & dt&L_\infty \text{ error} & \text{Convergence } L_\infty & L_2 \text{ error} & \text{Convergence } L_2 & H^1 \text{ semi-error} & \text{Convergence } H^1 \\ \hline 1/4 & 1/4 & 1/8 & 2.7228e-03 & & 1.1655e-03 & & 2.8968e-02 & \\ 1/8 & 1/8 & 1/64 & 3.6743e-04 & 2.8896& 1.5276e-04 & 2.9316& 7.1757e-03 & 2.0133 \\ 1/16 & 1/16 & 1/512 & 4.6325e-05 & 2.9876 & 1.9377e-05 & 2.9788 & 1.7893e-03 & 2.0037 \\ 1/32 & 1/32 & 1/4096 & 5.8420e-06 &2.9873 & 2.4375e-06 & 2.9909 & 4.4704e-04 & 2.0009 \\ \hline \end{array}$

Table4：Results for $w$ ，Using Quadratic basis function， Backforward Euler scheme（ $\theta =1$ ）
$\begin{array}{cccccccc} \hline h_1 & h_2 &dt & L_\infty \text{ error} & \text{Convergence } L_\infty & L_2 \text{ error} & \text{Convergence } L_2 & H^1 \text{ semi-error} & \text{Convergence } H^1 \\ \hline 1/4& 1/4& 1/8 & 3.6298e-02& &1.9126e-02& &1.1580e-01& \\ 1/8& 1/8 & 1/64 & 4.7576e-03 & 2.9316& 2.4981e-03&2.9366 & 2.2043e-02&2.3932\\ 1/16&1/16 & 1/512 & 6.0271e-04 & 2.9807 & 3.1466e-04 & 2.9889 & 4.9292e-03&2.1609 \\ 1/32&1/32 & 1/4096& 7.5641e-05 & 2.9942& 3.9415e-05 &2.9970 & 1.1932e-03& 2.0465\\ \hline \end{array}$