首页 > 其他分享 >有限元求解Cahn Hilliard 相场方程

有限元求解Cahn Hilliard 相场方程

时间：2024-06-23 10:59:17浏览次数：25  
标签：有限元 phi int varphi 相场 Hilliard dxdy Omega vec

Cahn Hilliard 方程, 相场

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

目标方程

∂ φ ∂ t − M c Δ w = f 1 w + r ϵ Δ φ = f 2 \frac{\partial \varphi }{\partial t} - M_c\Delta w = f_1\\ w + r\epsilon \Delta \varphi = f_2 ∂t∂φ​−Mc​Δw=f1​w+rϵΔφ=f2​
这里 φ ( x , y , t ) , w ( x , y , t ) \varphi(x,y,t),w(x,y,t) φ(x,y,t),w(x,y,t) 是需要求解的未知函数， f 1 ( x , y , t ) , f 2 ( x , y , t ) f_1(x,y,t),f_2(x,y,t) f1​(x,y,t),f2​(x,y,t)是已知函数， M c , r , ϵ M_c,r,\epsilon Mc​,r,ϵ是参数，为了简单，我全部设置为1. 求解区域为 Ω = [ 0 , 1 ] × [ 0 , 1 ] \Omega = [0,1] \times [0,1] Ω=[0,1]×[0,1],时间区域为 T d [ 0 , 1 ] T_d[0,1] Td​[0,1].

∇ w = ( w x , w y ) Δ w = ∇ ⋅ ∇ w = w x x + w y y \nabla w = (w_x,w_y)\\ \Delta w =\nabla \cdot \nabla w = w_{xx} + w_{yy} ∇w=(wx​,wy​)Δw=∇⋅∇w=wxx​+wyy​

∇ φ = ( φ x , φ y ) Δ φ = ∇ ⋅ ∇ φ = φ x x + φ y y \nabla \varphi = (\varphi_x,\varphi_y)\\ \Delta \varphi =\nabla \cdot \nabla \varphi = \varphi_{xx} + \varphi_{yy} ∇φ=(φx​,φy​)Δφ=∇⋅∇φ=φxx​+φyy​

弱格式

第一个方程

∂ φ ∂ t − M c Δ w = f 1 \frac{\partial \varphi }{\partial t} - M_c\Delta w = f_1 ∂t∂φ​−Mc​Δw=f1​

  1. Multiply a function a test function
    multiply a test function v ( x , y ) v(x,y) v(x,y) on both side of the the above equation.
    ∂ φ ∂ t v − ( M c Δ w ) v = f 1 v ∫ Ω ∂ φ ∂ t v d x d y − ∫ Ω ( M c Δ w ) v d x d y = ∫ Ω f 1 v d x d y \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 ∂t∂φ​v−(Mc​Δw)v=f1​v∫Ω​∂t∂φ​vdxdy−∫Ω​(Mc​Δw)vdxdy=∫Ω​f1​vdxdy
    φ ( x , y , t ) \varphi(x,y,t) φ(x,y,t) and w ( x , y , t ) w(x,y,t) w(x,y,t) are called trail function and v ( x , y ) v(x,y) v(x,y) is called a test function.

  2. Green’s formula
    using Green’s formula (divergence theory, integration by parts in multi-dimension).
    ∫ Ω ( M c Δ w ) v d x d y = ∫ Ω ( M c ∇ ⋅ ∇ w ) v d x d y = ∫ ∂ Ω ( M c ∇ w ⋅ n ⃗ ) v − ∫ Ω M c ∇ w ⋅ ∇ v d x d y \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 ∫Ω​(Mc​Δw)vdxdy=∫Ω​(Mc​∇⋅∇w)vdxdy=∫∂Ω​(Mc​∇w⋅n )v−∫Ω​Mc​∇w⋅∇vdxdy
    Assume it is Dirichlet boundary condtion. Then v = 0 v=0 v=0 on ∂ Ω \partial \Omega ∂Ω
    ∫ Ω ∂ φ ∂ t v d x d y + ∫ Ω M c ∇ w ⋅ ∇ v d x d y = ∫ Ω f 1 v d x d y \int_\Omega\frac{\partial \varphi }{\partial t}vdxdy +\int_\Omega M_c\nabla w \cdot \nabla vdxdy = \int_\Omega f_1vdxdy ∫Ω​∂t∂φ​vdxdy+∫Ω​Mc​∇w⋅∇vdxdy=∫Ω​f1​vdxdy
    ∇ w ⋅ ∇ v = w x v x + w y v y \nabla w \cdot \nabla v = w_xv_x+w_yv_y ∇w⋅∇v=wx​vx​+wy​vy​
    ∫ Ω ∂ φ ∂ t v d x d y + ∫ Ω M c w x v x d x d y + ∫ Ω M c w y v y d x d y = ∫ Ω f 1 v d x d 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 ∫Ω​∂t∂φ​vdxdy+∫Ω​Mc​wx​vx​dxdy+∫Ω​Mc​wy​vy​dxdy=∫Ω​f1​vdxdy

Let ( φ t , v ) = ∫ Ω ∂ φ ∂ t v d x d y (\varphi_t,v) = \int_\Omega\frac{\partial \varphi }{\partial t}vdxdy (φt​,v)=∫Ω​∂t∂φ​vdxdy
a 1 ( w , v ) = ∫ Ω M c ∇ w ⋅ ∇ v d x d y a1(w,v) =\int_\Omega M_c\nabla w \cdot \nabla vdxdy a1(w,v)=∫Ω​Mc​∇w⋅∇vdxdy
( f 1 , v ) = ∫ Ω f 1 v d x d y (f1,v) = \int_\Omega f_1vdxdy (f1,v)=∫Ω​f1​vdxdy \

H 1 ( 0 , T ; H 1 ( Ω ) ) = { v ( t , ⋅ )   |   ∂ v ∂ t ( t , ⋅ ) ∈ H 1 ( Ω ) ,   ∀ t ∈ [ 0 , T ] } 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\} H1(0,T;H1(Ω))={v(t,⋅) ​∂t∂v​(t,⋅)∈H1(Ω),∀t∈[0,T]}

Weak formulation: find w ∈ H 1 ( 0 , T ; H 1 ( Ω ) ) w\in H^1(0, T; H^1(\Omega)) w∈H1(0,T;H1(Ω)) and φ ∈ H 1 ( 0 , T ; H 1 ( Ω ) ) \varphi \in H^1(0, T; H^1(\Omega)) φ∈H1(0,T;H1(Ω)), such that ：
( φ t , v ) + a 1 ( w , v ) = ( f 1 , v ) (\varphi_t,v) + a1(w,v) = (f1,v) (φt​,v)+a1(w,v)=(f1,v)
for any v ∈ H 0 1 ( Ω ) v\in H^1_0(\Omega) v∈H01​(Ω)

第二个方程

w + r ϵ Δ φ = f 2 w + r\epsilon \Delta \varphi = f_2 w+rϵΔφ=f2​

  1. multiply a test function v v v
    multiply a function a test function v ( x , y ) v(x,y) v(x,y) on both side of the the above equation.
    w + r ϵ Δ φ = f 2 w + r\epsilon \Delta \varphi = f_2 w+rϵΔφ=f2​
    ∫ Ω w v d x d y + ∫ Ω r ϵ Δ φ v d x d y = ∫ Ω f 2 v d x d y \int_\Omega wvdxdy + \int_\Omega r\epsilon \Delta \varphi vdxdy = \int_\Omega f_2vdxdy ∫Ω​wvdxdy+∫Ω​rϵΔφvdxdy=∫Ω​f2​vdxdy
    φ ( x , y , t ) \varphi(x,y,t) φ(x,y,t) and w ( x , y , t ) w(x,y,t) w(x,y,t) are called trail function and v ( x , y ) v(x,y) v(x,y) is called a test function.

  2. Green’s formula
    using Green’s formula (divergence theory, integration by parts in multi-dimension).
    ∫ Ω r ϵ Δ φ v d x d y = ∫ Ω ( r ϵ ∇ ⋅ ∇ φ ) v d x d y = ∫ ∂ Ω ( r ϵ ∇ φ ⋅ n ⃗ ) v − ∫ Ω r ϵ ∇ φ ⋅ ∇ v d x d y \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 ∫Ω​rϵΔφvdxdy=∫Ω​(rϵ∇⋅∇φ)vdxdy=∫∂Ω​(rϵ∇φ⋅n )v−∫Ω​rϵ∇φ⋅∇vdxdy

    Assume it is Dirichlet boundary condtion. Then v = 0 v=0 v=0 on ∂ Ω \partial \Omega ∂Ω
    ∫ Ω w v d x d y − ∫ Ω r ϵ ∇ φ ⋅ ∇ v d x d y = ∫ Ω f 2 v d x d y \int_\Omega wvdxdy- \int_\Omega r\epsilon \nabla \varphi \cdot \nabla vdxdy = \int_\Omega f_2vdxdy ∫Ω​wvdxdy−∫Ω​rϵ∇φ⋅∇vdxdy=∫Ω​f2​vdxdy
    ∇ φ ⋅ ∇ v = φ x v x + φ y v y \nabla \varphi \cdot \nabla v = \varphi_xv_x+\varphi_yv_y ∇φ⋅∇v=φx​vx​+φy​vy​

    ∫ Ω w v d x d y − ∫ Ω r ϵ φ x v x d x d y − ∫ Ω r ϵ φ y v y d x d y = ∫ Ω f 2 v d x d y \int_\Omega w vdxdy -\int_\Omega r\epsilon \varphi_xv_x dxdy - \int_\Omega r\epsilon\varphi_yv_ydxdy = \int_\Omega f_2vdxdy ∫Ω​wvdxdy−∫Ω​rϵφx​vx​dxdy−∫Ω​rϵφy​vy​dxdy=∫Ω​f2​vdxdy
    Let
    a 2 ( w , v ) = ∫ Ω w v d x d y a2(w,v) = \int_\Omega w vdxdy a2(w,v)=∫Ω​wvdxdy
    a 3 ( φ , v ) = − ∫ Ω r ϵ ∇ φ ⋅ ∇ v d x d y a3(\varphi,v) = - \int_\Omega r\epsilon \nabla \varphi \cdot \nabla vdxdy a3(φ,v)=−∫Ω​rϵ∇φ⋅∇vdxdy

( f 2 , v ) = ∫ Ω f 2 v d x d y (f2,v) = \int_\Omega f_2vdxdy (f2,v)=∫Ω​f2​vdxdy \

Weak formulation: find w ∈ H 1 ( 0 , T ; H 1 ( Ω ) ) w\in H^1(0, T; H^1(\Omega)) w∈H1(0,T;H1(Ω)) and φ ∈ H 1 ( 0 , T ; H 1 ( Ω ) ) \varphi \in H^1(0, T; H^1(\Omega)) φ∈H1(0,T;H1(Ω)), such that \
a 2 ( w , v ) + a 3 ( φ , v ) = ( f 2 , v ) a2(w,v) + a3(\varphi,v) = (f2,v) a2(w,v)+a3(φ,v)=(f2,v)
for any v ∈ H 0 1 ( Ω ) v\in H^1_0(\Omega) v∈H01​(Ω)

半离散

第一个方程

Assume there is a finite dimensional subspace U h = span { ϕ j } j = 1 N b U_h = \text{span}\{\phi_j\}_{j=1}^{N_b} Uh​=span{ϕj​}j=1Nb​​,
where { ϕ j } j = 1 N b \{\phi_j\}_{j=1}^{N_b} {ϕj​}j=1Nb​​ are the global finite element basis functions.\

Then the Galerkin formulation is to find φ h ∈ H 1 ( 0 , T ; U h ) \varphi_{h} \in H^1(0,T;U_h) φh​∈H1(0,T;Uh​) and w h ∈ H 1 ( 0 , T ; U h ) w_{h} \in H^1(0,T;U_h) wh​∈H1(0,T;Uh​)
such that
( φ h t , v h ) + a 1 ( w h , v h ) = ( f 1 , v h ) (\varphi_{h_t},v_h) + a1(w_h,v_h) = (f1,v_h) (φht​​,vh​)+a1(wh​,vh​)=(f1,vh​)
for any v h ∈ U h v_h \in U_h vh​∈Uh​

It also can be express as:
∫ Ω φ h t v h d x d y + ∫ Ω M c ∇ w h ⋅ ∇ v h d x d y = ∫ Ω f 1 v h d x d y \int_\Omega\varphi_{h_t}v_hdxdy +\int_\Omega M_c\nabla w_h \cdot \nabla v_hdxdy = \int_\Omega f_1v_hdxdy ∫Ω​φht​​vh​dxdy+∫Ω​Mc​∇wh​⋅∇vh​dxdy=∫Ω​f1​vh​dxdy
or
∫ Ω φ h t v h d x d y + ∫ Ω M c ( w h ) x ( v h ) x d x d y + ∫ Ω M c ( w h ) y ( v h ) y d x d y = ∫ Ω f 1 v h d x d y \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 ∫Ω​φht​​vh​dxdy+∫Ω​Mc​(wh​)x​(vh​)x​dxdy+∫Ω​Mc​(wh​)y​(vh​)y​dxdy=∫Ω​f1​vh​dxdy

Since $\varphi_h \in H^1(0,T;U_h) $ and w h ∈ H 1 ( 0 , T ; U h ) w_{h} \in H^1(0,T;U_h) wh​∈H1(0,T;Uh​), and U h = span { ϕ j } j = 1 N b U_h = \text{span}\{\phi_j\}_{j=1}^{N_b} Uh​=span{ϕj​}j=1Nb​​,then
φ h ( x , y , t ) = ∑ j = 1 N b φ j ( t ) ϕ j ( x , y ) \varphi_{h}(x,y,t) = \sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j(x,y) φh​(x,y,t)=j=1∑Nb​​φj​(t)ϕj​(x,y)

w h ( x , y , t ) = ∑ j = 1 N b w j ( t ) ϕ j ( x , y ) w_{h}(x,y,t) = \sum_{j=1}^{N_b} w_{j}(t)\phi_j(x,y) wh​(x,y,t)=j=1∑Nb​​wj​(t)ϕj​(x,y)
for some coefficients $\varphi_{j}(t)(j=1,…,N_b) $ and w j ( t ) ( j = 1 , . . . , N b ) w_{j}(t)(j=1,...,N_b) wj​(t)(j=1,...,Nb​).\
We choose V h = ϕ i ( x , y ) ( i = 1 , . . . , N b ) V_h = \phi_i(x,y)(i=1,...,N_b) Vh​=ϕi​(x,y)(i=1,...,Nb​).\

Then we can have

∫ Ω ( ∑ j = 1 N b φ j ( t ) ϕ j ) t ϕ i d x d y + ∫ Ω M c ∇ ( ∑ j = 1 N b w j ( t ) ϕ j ) ⋅ ∇ ϕ i d x d y = ∫ Ω f 1 ϕ i d x d y \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 ∫Ω​(j=1∑Nb​​φj​(t)ϕj​)t​ϕi​dxdy+∫Ω​Mc​∇(j=1∑Nb​​wj​(t)ϕj​)⋅∇ϕi​dxdy=∫Ω​f1​ϕi​dxdy

∫ Ω ( ∑ j = 1 N b φ j ( t ) ϕ j ) t ϕ i d x d y + ∫ Ω M c ∑ j = 1 N b w j ( t ) ∇ ϕ j ⋅ ∇ ϕ i d x d y = ∫ Ω f 1 ϕ i d x d y \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 ∫Ω​(j=1∑Nb​​φj​(t)ϕj​)t​ϕi​dxdy+∫Ω​Mc​j=1∑Nb​​wj​(t)∇ϕj​⋅∇ϕi​dxdy=∫Ω​f1​ϕi​dxdy

∑ j = 1 N b φ j ′ ( t ) ∫ Ω ϕ j ϕ i d x d y + ∑ j = 1 N b w j ( t ) ∫ Ω M c ∇ ϕ j ⋅ ∇ ϕ i d x d y = ∫ Ω f 1 ϕ i d x d y \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 j=1∑Nb​​φj′​(t)∫Ω​ϕj​ϕi​dxdy+j=1∑Nb​​wj​(t)∫Ω​Mc​∇ϕj​⋅∇ϕi​dxdy=∫Ω​f1​ϕi​dxdy

∑ j = 1 N b φ j ′ ( t ) ∫ Ω ϕ j ϕ i d x d y + ∑ j = 1 N b w j ( t ) ( ∫ Ω M c ∂ ϕ j ∂ x ∂ ϕ i ∂ x d x d y + ∫ Ω M c ∂ ϕ j ∂ y ∂ ϕ i ∂ y d x d y ) = ∫ Ω f 1 ϕ i d x d y \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 j=1∑Nb​​φj′​(t)∫Ω​ϕj​ϕi​dxdy+j=1∑Nb​​wj​(t)(∫Ω​Mc​∂x∂ϕj​​∂x∂ϕi​​dxdy+∫Ω​Mc​∂y∂ϕj​​∂y∂ϕi​​dxdy)=∫Ω​f1​ϕi​dxdy

Define the stiffness matrix
A 1 = [ ∫ Ω M c ∇ ϕ j ⋅ ∇ ϕ i d x d y ] i , j = 1 N b , N b A_1= \left[\int_\Omega M_c \nabla \phi_j \cdot \nabla \phi_i dxdy\right]_{i,j=1}^{N_b,N_b} A1​=[∫Ω​Mc​∇ϕj​⋅∇ϕi​dxdy]i,j=1Nb​,Nb​​
A 1 = A 11 + A 12 A_1= A_{11} + A_{12} A1​=A11​+A12​
A 11 = [ ∫ Ω M c ∂ ϕ j ∂ x ∂ ϕ i ∂ x d x d y ] i , j = 1 N b , N b 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} A11​=[∫Ω​Mc​∂x∂ϕj​​∂x∂ϕi​​dxdy]i,j=1Nb​,Nb​​
A 12 = [ ∫ Ω M c ∂ ϕ j ∂ y ∂ ϕ i ∂ y d x d y ] 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} A12​=[∫Ω​Mc​∂y∂ϕj​​∂y∂ϕi​​dxdy]i,j=1Nb​,Nb​​
Define the mass matrix
M = [ ∫ Ω ϕ j ϕ i d x d y ] i , j = 1 N b , N b M = \left[\int_\Omega\phi_j\phi_idxdy\right]_{i,j=1}^{N_b,N_b} M=[∫Ω​ϕj​ϕi​dxdy]i,j=1Nb​,Nb​​

Define the load vector
b ⃗ 1 = [ ∫ Ω f 1 ϕ i d x d y ] i = 1 N b \vec b_1 = \left[\int_\Omega f_1\phi_idxdy \right]_{i=1}^{N_b} b 1​=[∫Ω​f1​ϕi​dxdy]i=1Nb​​
Define the unknown vector
X ⃗ 1 = [ φ j ( t ) ] j = 1 N b \vec X_1=\left[\varphi_j(t) \right]_{j=1}^{N_b} X 1​=[φj​(t)]j=1Nb​​
X ⃗ 2 = [ w j ( t ) ] j = 1 N b \vec X_2=\left[w_j(t) \right]_{j=1}^{N_b} X 2​=[wj​(t)]j=1Nb​​

Then we obtain the system
M X ⃗ 1 ′ ( t ) + A 1 X ⃗ 2 = b ⃗ 1 M\vec X_1'(t) + A_1\vec X_2 = \vec b_1 MX 1′​(t)+A1​X 2​=b 1​

第二个方程

The Galerkin formulation is to find φ h ∈ H 1 ( 0 , T ; U h ) \varphi_{h} \in H^1(0,T;U_h) φh​∈H1(0,T;Uh​) and w h ∈ H 1 ( 0 , T ; U h ) w_{h} \in H^1(0,T;U_h) wh​∈H1(0,T;Uh​)
such that
a 2 ( w h , v h ) + a 3 ( φ h , v h ) = ( f 2 , v h ) a2(wh,vh) + a3(\varphi_h,v_h) = (f2,v_h) a2(wh,vh)+a3(φh​,vh​)=(f2,vh​)
for any v h ∈ U h v_h \in U_h vh​∈Uh​

It also can be express as:
∫ Ω w h v h d x d y − ∫ Ω r ϵ ∇ φ h ⋅ ∇ v h d x d y = ∫ Ω f 2 v h d x d y \int_\Omega w_hv_hdxdy- \int_\Omega r\epsilon \nabla \varphi_h \cdot \nabla v_hdxdy = \int_\Omega f_2v_hdxdy ∫Ω​wh​vh​dxdy−∫Ω​rϵ∇φh​⋅∇vh​dxdy=∫Ω​f2​vh​dxdy
or
∫ Ω w h v h d x d y − ∫ Ω r ϵ ( φ h ) x ( v h ) x d x d y − ∫ Ω r ϵ ( φ h ) y ( v h ) y d x d y = ∫ Ω f 2 v h d x d y \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 ∫Ω​wh​vh​dxdy−∫Ω​rϵ(φh​)x​(vh)x​dxdy−∫Ω​rϵ(φh​)y​(vh)y​dxdy=∫Ω​f2​vh​dxdy

Since
φ h ( x , y , t ) = ∑ j = 1 N b φ j ( t ) ϕ j ( x , y ) \varphi_{h}(x,y,t) = \sum_{j=1}^{N_b} \varphi_{j}(t)\phi_j(x,y) φh​(x,y,t)=j=1∑Nb​​φj​(t)ϕj​(x,y)

w h ( x , y , t ) = ∑ j = 1 N b w j ( t ) ϕ j ( x , y ) w_{h}(x,y,t) = \sum_{j=1}^{N_b} w_{j}(t)\phi_j(x,y) wh​(x,y,t)=j=1∑Nb​​wj​(t)ϕj​(x,y)
for some coefficients $\varphi_{j}(t)(j=1,…,N_b) $ and w j ( t ) ( j = 1 , . . . , N b ) w_{j}(t)(j=1,...,N_b) wj​(t)(j=1,...,Nb​).\
We choose V h = ϕ i ( x , y ) ( i = 1 , . . . , N b ) V_h = \phi_i(x,y)(i=1,...,N_b) Vh​=ϕi​(x,y)(i=1,...,Nb​).\

∫ Ω ( ∑ j = 1 N b w j ( t ) ϕ j ) ϕ i d x d y − ∫ Ω r ϵ ∇ ( ∑ j = 1 N b φ j ( t ) ϕ j ) ⋅ ∇ ϕ i d x d y = ∫ Ω f 2 ϕ i d x d y \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 ∫Ω​(j=1∑Nb​​wj​(t)ϕj​)ϕi​dxdy−∫Ω​rϵ∇(j=1∑Nb​​φj​(t)ϕj​)⋅∇ϕi​dxdy=∫Ω​f2​ϕi​dxdy

∑ j = 1 N b w j ( t ) ∫ Ω ϕ j ϕ i d x d y + ∑ j = 1 N b φ j ( t ) ∫ Ω − r ϵ ∇ ϕ j ⋅ ∇ ϕ i d x d y = ∫ Ω f 2 ϕ i d x d y \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 j=1∑Nb​​wj​(t)∫Ω​ϕj​ϕi​dxdy+j=1∑Nb​​φj​(t)∫Ω​−rϵ∇ϕj​⋅∇ϕi​dxdy=∫Ω​f2​ϕi​dxdy

∑ j = 1 N b w j ( t ) ∫ Ω ϕ j ϕ i d x d y + ∑ j = 1 N b φ j ( t ) ( ∫ Ω − r ϵ ∂ ϕ j ∂ x ∂ ϕ i ∂ x d x d y + ∫ Ω − r ϵ ∂ ϕ j ∂ y ∂ ϕ i ∂ y d x d y ) = ∫ Ω f 2 ϕ i d x d y \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 j=1∑Nb​​wj​(t)∫Ω​ϕj​ϕi​dxdy+j=1∑Nb​​φj​(t)(∫Ω​−rϵ∂x∂ϕj​​∂x∂ϕi​​dxdy+∫Ω​−rϵ∂y∂ϕj​​∂y∂ϕi​​dxdy)=∫Ω​f2​ϕi​dxdy

Define the stiffness matrix
A 2 = [ ∫ Ω − r ϵ ∇ ϕ j ⋅ ∇ ϕ i d x d y ] i , j = 1 N b , N b A_2= \left[\int_\Omega - r\epsilon \nabla \phi_j\cdot \nabla \phi_idxdy\right]_{i,j=1}^{N_b,N_b} A2​=[∫Ω​−rϵ∇ϕj​⋅∇ϕi​dxdy]i,j=1Nb​,Nb​​
A 2 = A 21 + A 22 A_2= A_{21} + A_{22} A2​=A21​+A22​
A 21 = [ ∫ Ω − r ϵ ∂ ϕ j ∂ x ∂ ϕ i ∂ x d x d y ] i , j = 1 N b , N b 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} A21​=[∫Ω​−rϵ∂x∂ϕj​​∂x∂ϕi​​dxdy]i,j=1Nb​,Nb​​
A 22 = [ ∫ Ω − r ϵ ∂ ϕ j ∂ y ∂ ϕ i ∂ y d x d y ] 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} A22​=[∫Ω​−rϵ∂y∂ϕj​​∂y∂ϕi​​dxdy]i,j=1Nb​,Nb​​
Define the mass matrix
M = [ ∫ Ω ϕ j ϕ i d x d y ] i , j = 1 N b , N b M = \left[\int_\Omega\phi_j\phi_idxdy\right]_{i,j=1}^{N_b,N_b} M=[∫Ω​ϕj​ϕi​dxdy]i,j=1Nb​,Nb​​

Define the load vector
b ⃗ 2 = [ ∫ Ω f 2 ϕ i d x d y ] i = 1 N b \vec b_2 = \left[\int_\Omega f_2\phi_idxdy \right]_{i=1}^{N_b} b 2​=[∫Ω​f2​ϕi​dxdy]i=1Nb​​
Define the unknown vector
X ⃗ 1 = [ φ j ( t ) ] j = 1 N b \vec X_1=\left[\varphi_j(t) \right]_{j=1}^{N_b} X 1​=[φj​(t)]j=1Nb​​
X ⃗ 2 = [ w j ( t ) ] j = 1 N b \vec X_2=\left[w_j(t) \right]_{j=1}^{N_b} X 2​=[wj​(t)]j=1Nb​​

Then we obtain the system
A 2 X ⃗ 1 + M X ⃗ 2 = b ⃗ 2 A_2\vec X_1+ M\vec X_2=\vec b_2 A2​X 1​+MX 2​=b 2​

全离散

第一个方程

M X ⃗ 1 ′ ( t ) + A 1 X ⃗ 2 = b ⃗ 1 M\vec X_1'(t) + A_1\vec X_2 = \vec b_1 MX 1′​(t)+A1​X 2​=b 1​

X ⃗ 1 ′ ( t ) = b 1 − A 1 X ⃗ 2 M \vec X_1'(t) = \frac{b_1 - A_1\vec X_2}{M} X 1′​(t)=Mb1​−A1​X 2​​

Assume that we have a uniform partition of [ 0 , 1 ] [0,1] [0,1] into M m M_m Mm​ elements with mesh size Δ t \Delta t Δt.

The mesh nodes are t m = m Δ t , m = 0 , 1 , . . . , M m t_m = m \Delta t, m = 0,1,...,M_m tm​=mΔt,m=0,1,...,Mm​.
Assume X ⃗ 1 m \vec X_1^{m} X 1m​ is the numerical solution of X ⃗ 1 ( t m ) \vec X_1(t_m) X 1​(tm​)
Assume X ⃗ 2 m \vec X_2^{m} X 2m​ is the numerical solution of X ⃗ 2 ( t m ) \vec X_2(t_m) X 2​(tm​)

θ \theta θ-scheme (Finite difference (FD) method):
X 1 m + 1 − X 1 m Δ t = θ b 1 m + 1 − A 1 m + 1 X 2 m + 1 M + ( 1 − θ ) b 1 m − A 1 m X 2 m M , m = 0 , 1 , . . . , M m \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 ΔtX1m+1​−X1m​​=θMb1m+1​−A1m+1​X2m+1​​+(1−θ)Mb1m​−A1m​X2m​​,m=0,1,...,Mm​
where:
θ = 0 \theta = 0 θ=0: forward Euler scheme;
θ = 1 \theta = 1 θ=1: backward Euler scheme;
θ = 1 / 2 \theta =1/2 θ=1/2: Crank-Nicolson scheme

Then we can have the ireration scheme:

M Δ t X ⃗ 1 m + 1 + θ A 1 m + 1 X ⃗ 2 m + 1 = M Δ t X ⃗ 1 m + θ b ⃗ 1 m + 1 + ( 1 − θ ) ( b ⃗ 1 m − A 1 m X ⃗ 2 m ) \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}) ΔtM​X 1m+1​+θA1m+1​X 2m+1​=ΔtM​X 1m​+θb 1m+1​+(1−θ)(b 1m​−A1m​X 2m​)

第二个方程

A 2 m + 1 X ⃗ 1 m + 1 + M X ⃗ 2 m + 1 = b ⃗ 2 m + 1 A_2^{m+1}\vec X_1^{m+1}+ M\vec X_2^{m+1}=\vec b_2^{m+1} A2m+1​X 1m+1​+MX 2m+1​=b 2m+1​

So the systems becomes to

M Δ t X ⃗ 1 m + 1 + θ A 1 m + 1 X ⃗ 2 m + 1 = M Δ t X ⃗ 1 m + θ b ⃗ 1 m + 1 + ( 1 − θ ) ( b ⃗ 1 m − A 1 m X ⃗ 2 m ) A 2 m + 1 X ⃗ 1 m + 1 + M X ⃗ 2 m + 1 = b ⃗ 2 m + 1 \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} ΔtM​X 1m+1​+θA1m+1​X 2m+1​=ΔtM​X 1m​+θb 1m+1​+(1−θ)(b 1m​−A1m​X 2m​)A2m+1​X 1m+1​+MX 2m+1​=b 2m+1​

A m + 1 = [ M Δ t θ A 1 m + 1 A 2 m + 1 M ] 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} Am+1= ​ΔtM​A2m+1​​θA1m+1​M​ ​

X ⃗ m + 1 = [ X ⃗ 1 m + 1 X ⃗ 2 m + 1 ] \vec{X}^{m+1} = \begin{bmatrix}\displaystyle \vec{X}_1^{m+1}\\ \vec{X}_2^{m+1}\end{bmatrix} X m+1=[X 1m+1​X 2m+1​​]

b ⃗ m + 1 = [ M Δ t X ⃗ 1 m + θ b ⃗ 1 m + 1 + ( 1 − θ ) ( b ⃗ 1 m − A 1 m X ⃗ 2 m ) b ⃗ 2 m + 1 ] \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} b m+1= ​ΔtM​X 1m​+θb 1m+1​+(1−θ)(b 1m​−A1m​X 2m​)b 2m+1​​ ​

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

结果

The exact solution: φ = e x + y + t \varphi = e^{x+y+t} φ=ex+y+t, w = x y ( 1 − x 2 ) ( 1 − y ) e x + y + t w=xy(1-\frac{x}{2})(1-y)e^{x+y+t} w=xy(1−2x​)(1−y)ex+y+t
Set M c = 1 Mc=1 Mc=1, r = 1 r=1 r=1, ϵ = 1 \epsilon=1 ϵ=1, θ = 1 \theta = 1 θ=1, then compute f 1 f_1 f1​ and f 2 f_2 f2​ 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（ θ = 1 \theta =1 θ=1）
h 1 h 2 d t L ∞  error Convergence  L ∞ L 2  error Convergence  L 2 H 1  semi-error Convergence  H 1 1 / 4 1 / 4 1 / 8 1.3626 e − 01 4.9320 e − 02 8.8899 e − 01 1 / 8 1 / 8 1 / 64 3.6311 e − 02 1.9079 1.2185 e − 02 2.0171 4.4348 e − 01 1.0033 1 / 16 1 / 16 1 / 512 9.3747 e − 03 1.9536 3.0301 e − 03 2.0077 2.2161 e − 01 1.0008 1 / 32 1 / 32 1 / 4096 2.3818 e − 03 1.9767 7.5579 e − 04 2.0033 1.1079 e − 01 1.0002 \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} h1​1/41/81/161/32​h2​1/41/81/161/32​dt1/81/641/5121/4096​L∞​ error1.3626e−013.6311e−029.3747e−032.3818e−03​Convergence L∞​1.90791.95361.9767​L2​ error4.9320e−021.2185e−023.0301e−037.5579e−04​Convergence L2​2.01712.00772.0033​H1 semi-error8.8899e−014.4348e−012.2161e−011.1079e−01​Convergence H11.00331.00081.0002​​

Table2：Results for w w w ，Using linear basis function， Backforward Euler scheme（ θ = 1 \theta =1 θ=1）
h 1 h 2 d t L ∞  error Convergence  L ∞ L 2  error Convergence  L 2 H 1  semi-error Convergence  H 1 1 / 4 1 / 4 1 / 8 2.2236 e − 01 4.4761 e − 02 7.6775 e − 01 1 / 8 1 / 8 1 / 64 6.4205 e − 02 1.7922 1.2208 e − 02 1.8771 3.8709 e − 01 9.8796 e − 01 1 / 16 1 / 16 1 / 512 1.7240 e − 02 1.8969 3.2441 e − 03 2.0104 1.9411 e − 01 9.9579 e − 01 1 / 32 1 / 32 1 / 4096 4.4660 e − 03 1.9487 8.3737 e − 04 2.1018 9.7133 e − 02 9.9886 e − 01 \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} h1​1/41/81/161/32​h2​1/41/81/161/32​dt1/81/641/5121/4096​L∞​ error2.2236e−016.4205e−021.7240e−024.4660e−03​Convergence L∞​1.79221.89691.9487​L2​ error4.4761e−021.2208e−023.2441e−038.3737e−04​Convergence L2​1.87712.01042.1018​H1 semi-error7.6775e−013.8709e−011.9411e−019.7133e−02​Convergence H19.8796e−019.9579e−019.9886e−01​​

Table3：：Results for φ \varphi φ, Using Quadratic basis function， Backforward Euler scheme（ θ = 1 \theta =1 θ=1）
h 1 h 2 d t L ∞  error Convergence  L ∞ L 2  error Convergence  L 2 H 1  semi-error Convergence  H 1 1 / 4 1 / 4 1 / 8 2.7228 e − 03 1.1655 e − 03 2.8968 e − 02 1 / 8 1 / 8 1 / 64 3.6743 e − 04 2.8896 1.5276 e − 04 2.9316 7.1757 e − 03 2.0133 1 / 16 1 / 16 1 / 512 4.6325 e − 05 2.9876 1.9377 e − 05 2.9788 1.7893 e − 03 2.0037 1 / 32 1 / 32 1 / 4096 5.8420 e − 06 2.9873 2.4375 e − 06 2.9909 4.4704 e − 04 2.0009 \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} h1​1/41/81/161/32​h2​1/41/81/161/32​dt1/81/641/5121/4096​L∞​ error2.7228e−033.6743e−044.6325e−055.8420e−06​Convergence L∞​2.88962.98762.9873​L2​ error1.1655e−031.5276e−041.9377e−052.4375e−06​Convergence L2​2.93162.97882.9909​H1 semi-error2.8968e−027.1757e−031.7893e−034.4704e−04​Convergence H12.01332.00372.0009​​

Table4：Results for w w w，Using Quadratic basis function， Backforward Euler scheme（ θ = 1 \theta =1 θ=1）
h 1 h 2 d t L ∞  error Convergence  L ∞ L 2  error Convergence  L 2 H 1  semi-error Convergence  H 1 1 / 4 1 / 4 1 / 8 3.6298 e − 02 1.9126 e − 02 1.1580 e − 01 1 / 8 1 / 8 1 / 64 4.7576 e − 03 2.9316 2.4981 e − 03 2.9366 2.2043 e − 02 2.3932 1 / 16 1 / 16 1 / 512 6.0271 e − 04 2.9807 3.1466 e − 04 2.9889 4.9292 e − 03 2.1609 1 / 32 1 / 32 1 / 4096 7.5641 e − 05 2.9942 3.9415 e − 05 2.9970 1.1932 e − 03 2.0465 \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} h1​1/41/81/161/32​h2​1/41/81/161/32​dt1/81/641/5121/4096​L∞​ error3.6298e−024.7576e−036.0271e−047.5641e−05​Convergence L∞​2.93162.98072.9942​L2​ error1.9126e−022.4981e−033.1466e−043.9415e−05​Convergence L2​2.93662.98892.9970​H1 semi-error1.1580e−012.2043e−024.9292e−031.1932e−03​Convergence H12.39322.16092.0465​​

其实 可以先不看那些空间定义什么的，步骤就是很简单，先乘 v v v, 然后积分化简， 得到离散格式，组装矩阵，求解 就完了。

标签：有限元,phi,int,varphi,相场,Hilliard,dxdy,Omega,vec
From： https://blog.csdn.net/wv0uue68/article/details/139891876

相关文章

  • ANSYS有限元网格划分的基本原则
    一、引言ANSYS有限元网格划分是进行数值模拟分析至关重要的一步，它直接影响着后续数值计算分析结果的精确性。网格划分涉及单元的形状及其拓扑类型、单元类型、网格生成器的选择、网格的密度、单元的编号以及几何体素。从几何表达上讲，梁和杆是相同的，从物理和数值求解上讲则是......
  • 基于matlab的动力学有限元期刊文章（关节接触界面的随机建模和更新）---论文复现
      stochasticmodellingandupdatingofajointcontactinterface是发表在mechanicalsystemsandsignalprocessing上的一篇较高质量文章。笔者成功复现该文章，效果优异，配备原文，方便学习使用，适合该方向的学习者。  接头和机械连接中接触界面的动态性能对装配结构的......
  • 什么是晶体塑性有限元？
    晶体塑性有限元是一种结合了晶体塑性理论和有限元方法的数值模拟技术。它主要用于分析和预测晶体材料的塑性变形行为，特别是在微观尺度上的变形机制。这种方法考虑了晶体材料的各向异性、滑移系统的开动和相互作用、以及变形过程中的硬化效应。在晶体塑性有限元中，材料被建模为包含......
  • 有限元法之有限元法的实现
    目录一、单元刚度矩阵及单元荷载二、总刚度矩阵及总荷载的合成图1三角形剖分三、边界条件处理四、算例实现4.1C++代码4.2计算结果五、结论        前三节我们介绍了有限元的基本概念、变分理论及有限元空间的构造，本节我们探讨如何实现有限元法。我们继续......
  • 二维有限元，线性插值
    设置u=-（x*x+y*y），c=(x+y),可得f=6*（x+y），设置所有边界条件为dirichlet边界条件，其他条件应该也不复杂。boundaryedge矩阵是自己对着生成网格给出来的。感觉最难的地方就是在计算单元刚度矩阵的时候，因为使用了坐标变换，变成平面的标准三角形。（xi,yi),(xj,yj),(xm,ym)分别对应到（0，0），（1，0），（0，1）......
  • 有限元分析与应用 | Finite Element Method (FEM) Analysis and Applications
    第1讲引论/1.2变形体力学的要点https://learning.edx.org/course/course-v1:TsinghuaX+70120073x+1T2024/block-v1:TsinghuaX+70120073x+1T2024+type@sequential+block@5c00cb7f61af4dc8abb857abadc46151/block-v1:TsinghuaX+70120073x+1T2024+type@vertical+block@579410847......
  • 有限元方法[Matlab]-笔记
    <<MATLABCodesforFiniteElementAnalysis-SolidsandStructures(Ferreira)>>笔记chapter01matlabbasic略第二章:离散系统笔记、例题Matlab代码problem1.m%MATLABcodesforFiniteElementAnalysis%Problem1:3springsproblem%clearme......
  • 有限元仿真 有限体积仿真
    0前言声明：此篇博客仅用于个人学习记录之用，并非是分享代码。HornortheCode本来想仔细学一下有限元再进行LAB3的制作，之前也是这么做的。问题是咱是做游戏的，不是做cae的，养家糊口之后真没有时间去研究那些东西。所以后面的像张量分析之类的也就没仔细去看。只是做了矩阵微积......
  • 有限元入门教程
    Introductiontofiniteelementmethod?有限元与卷积的区别和联系：如果把有限元网格当成图像，把xxx当成卷积核，对有限元进行积分等于对图像进行卷积操作吗？对图像进行卷积操作之后得到的是特征图像，对有限元进行积分得到的是？。参考：https://mp.weixin.qq.com/s?__biz=MzU2NDgzNjQxMw=......
  • abaqus6.14下载 Abaqus(有限元分析软件) 安装包下载方式
    Abaqus2020是一款广泛使用的有限元分析软件，它可以帮助用户解决各种复杂的工程和科学问题。作为一款先进的计算工具，Abaqus2020拥有强大的建模与仿真能力，可以用于各种不同领域的应用，包括汽车工程、航空航天、土木工程、生物医学、电子、机械工程等等。软件地址：看置顶贴abaqus2016特......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99