计算科学 - 双调和方程 - 吾爱随笔录

双调和方程

计算科学 pde 椭圆pde

2021-12-09 12:49:31

我想用数值求解双调和方程，即：

Δ^{2} u = f i n Ω

$\Delta^2 u=f~~in~~\Omega$

u = g_{1} o n \partial Ω

$u=g_1~~on ~~\partial \Omega$

\frac{\partial u}{\partial n} = g_{2} o n \partial Ω

$\frac {\partial u}{\partial n}=g_2~~on ~~\partial \Omega$

使用格林公式，我们有如果 $u$ 足够光滑

\int_{\partial Ω} \frac{\partial Δ u}{\partial n} v - \int_{Ω} \nabla (Δ u) . \nabla v d x = \int_{\partial Ω} \frac{\partial Δ u}{\partial n} v - \int_{\partial Ω} \frac{\partial v}{\partial n} Δ u + \int_{Ω} Δ u Δ v = \int_{Ω} f v

$\int_{\partial \Omega}\frac {\partial \Delta u}{\partial n}v-\int_{\Omega}\nabla(\Delta u).\nabla v dx =\int_{\partial \Omega}\frac {\partial \Delta u}{\partial n}v-\int_{\partial \Omega}\frac {\partial v }{\partial n}\Delta u+\int_{\Omega}\Delta u \Delta v=\int_{\Omega}fv$

测试空间是 $H^2(\Omega)$ , 我怎么能使用给定的边界条件 $g_1,g_2$ .

我怎样才能找到新条款 $\frac {\partial v }{\partial n}$ 和 $\frac {\partial \Delta u}{\partial n}$

4个回答

双调和方程是拉普拉斯能量的欧拉-拉格朗日方程 $\frac{1}{2} \langle \Delta u,\Delta u \rangle$ . 离散化高阶问题的系统方法是将无约束问题转换为有约束问题：最小化 $\frac{1}{2} \langle v,v\rangle$ 英石 $\Delta u=v$ ; 那是，

\frac{1}{2} ⟨ v, v ⟩ + ⟨ λ, Δ u - v ⟩ \frac{1}{2} ⟨ v, v ⟩ + ⟨ λ, \frac{\partial u}{\partial n} ⟩ - ⟨ λ, v ⟩ - ⟨ \nabla λ, \nabla v ⟩

$\begin{equation} \frac{1}{2} \langle v,v\rangle +\langle \lambda, \Delta u -v\rangle \\ \frac{1}{2} \langle v,v\rangle + \langle \lambda,\frac{\partial u}{\partial n}\rangle -\langle \lambda,v\rangle -\langle \nabla \lambda, \nabla v\rangle \end{equation}$ 我们在第二个方程中使用了格林的恒等式，并且

\frac{\partial}{\partial n}

$\frac{\partial }{\partial n}$ 是正规导数。

新的优化问题可以使用分段线性元素离散化，这通常用于二阶问题。这是混合有限元离散化的共同点。

我决定将我之前的评论扩展为答案。我建议使用 Morley 元素 $P_2$ 每个元素的基础和自由度是

顶点的值
边缘中点的正态导数

它可能是双谐波问题实施明智的最简单的元素。

这是使用sp.fem 库（版本 8019aa7）的示例：

from spfem.mesh import MeshTri                                                                                                                                                                                                           
from spfem.asm import AssemblerGlobal                                                                                                                                                                                                    
from spfem.element import ElementGlobalMorley                                                                                                                                                                                            
from spfem.utils import direct,const_cell                                                                                                                                                                                                
import numpy as np                                                                                                                                                                                                                       

# define mesh                                                                                                                                                                                                                            
m=MeshTri()                                                                                                                                                                                                                              
m.refine(5)                                                                                                                                                                                                                              

# define element and assemble the stiffness matrix and load vector                                                                                                                                                                       
e1=AbstractElementMorley()                                                                                                                                                                                                                 
a=AssemblerAbstract(m,e1)                                                                                                                                                                                                                  

def bilinear(ddu,ddv):                                                                                                                                                                                                                   
    return ddu[0][0]*ddv[0][0]+\                                                                                                                                                                                                         
           ddu[0][1]*ddv[0][1]+\                                                                                                                                                                                                         
           ddu[1][0]*ddv[1][0]+\                                                                                                                                                                                                         
           ddu[1][1]*ddv[1][1]                                                                                                                                                                                                           

def F(x):                                                                                                                                                                                                                                
    return 50                                                                                                                                                                                                                            

A=a.iasm(bilinear)                                                                                                                                                                                                                       
f=a.iasm(lambda v,x: F(x)*v)                                                                                                                                                                                                             

# find boundary and interior DOFs                                                                                                                                                                                                        
bnd_nodes=m.boundary_nodes()                                                                                                                                                                                                             
bnd_facets=m.boundary_facets()                                                                                                                                                                                                           
bnd_facets_mid=0.5*(m.p[:,m.facets[0,bnd_facets]]+\                                                                                                                                                                                      
                    m.p[:,m.facets[1,bnd_facets]])                                                                                                                                                                                       

d1=a.dofnum_u.getdofs(N=bnd_nodes)                                                                                                                                                                                                       
d2=a.dofnum_u.getdofs(F=bnd_facets)                                                                                                                                                                                                      

I=np.setdiff1d(np.arange(a.dofnum_u.N),np.union1d(d1,d2))                                                                                                                                                                                

# set boundary conditions                                                                                                                                                                                                                
x=np.zeros(A.shape[0])                                                                                                                                                                                                                   

def G1(x,y):                                                                                                                                                                                                                             
    return 0.05*np.sin(4*np.pi*y)                                                                                                                                                                                                        

def G2(x,y):                                                                                                                                                                                                                             
    return np.sin(4*np.pi*x)                                                                                                                                                                                                             

# check orientation of boundary normals from sign of jacobian determinant                                                                                                                                                                
tris=np.ones(m.t.shape[1])                                                                                                                                                                                                               
flip=np.nonzero((a.mapping.detDF(np.array([[0]]))<0).flatten())[0]                                                                                                                                                                       
tris[flip]=-1.0                                                                                                                                                                                                                          
ori=tris[m.f2t[0,bnd_facets]]                                                                                                                                                                                                            

x[d1]=G1(m.p[0,bnd_nodes],m.p[1,bnd_nodes])                                                                                                                                                                                              
x[d2]=ori*G2(bnd_facets_mid[0,:],bnd_facets_mid[1,:])                                                                                                                                                                                    

# solve the system                                                                                                                                                                                                                       
x=direct(A,f,I=I,x=x)                                                                                                                                                                                                                    

# draw the solution                                                                                                                                                                                                                      
xverts=x[a.dofnum_u.n_dof[0,:]]                                                                                                                                                                                                          
m.plot3(xverts)                                                                                                                                                                                                                          
m.show()

解决方案如下所示：

请注意，两者 $g_1$ 和 $g_2$ 很容易定义为狄利克雷型边界条件。

使用有限元方法，您还可以选择使用两种方法直接离散化弱公式：

使用符合 H² 的元素，例如三角形上的 Argyris 元素或 Hsieh-Clough-Tocher 或四边形上的 Bogner-Fox-Schmidt 元素。
结合所谓的 C ⁰内部惩罚方法使用常规 H¹ 符合元素。

如果您采用通过格林公式得出的弱公式，则选择测试空间 $v = \partial_n v = 0$ on the boundary 使边界项消失。

另一种非常适合高阶微分方程的方法是等几何分析。

双谐波问题的变分或弱公式如下： $u \in V_D$ 这样

a (u, v) = ℓ (v), \forall v \in V_{0},

$a(u,v) = \ell(v), \quad \forall v \in V_0,$ 其中双线性和线性形式由下式给出

a (u, v) = \int_{Ω} Δ u Δ v d Ω and ℓ (v) = \int_{Ω} f v d Ω,

$a(u,v) = \int_\Omega \Delta u \Delta v \,d\Omega \quad \text{and} \quad \ell(v) = \int_\Omega f v\, d\Omega,$ 超平面和测试空间由下式给出

V_{D} := {v \in H^{2} (Ω), v = g_{1}, \frac{\partial v}{\partial n} = g_{2} o n \partial Ω}

$V_D :=\{ v \in H^2(\Omega), v = g_1, \frac{\partial v}{\partial n} = g_2 \quad on \quad \partial \Omega\}$ 和

V_{0} := {v \in H^{2} (Ω), v = 0, \frac{\partial v}{\partial n} = 0 o n \partial Ω} .

$V_0 :=\{ v \in H^2(\Omega), v = 0, \frac{\partial v}{\partial n} = 0 \quad on \quad \partial \Omega\}.$

上的高阶导数 $\partial \Omega$ 实际上是 Neumann 边界条件。这些没有给出，因此您选择的测试功能就足够了。这是具有狄利克雷边界条件的双调和问题。

其它你可能感兴趣的问题

上一篇用于运行长时间和计算密集型仿真的外部硬件资源下一篇收敛到错误值的 Tanh-sinh 正交数值积分方法