biharmonic equation

Question

I want to solve the biharmonic equation numerically, that is:

$$\Delta^2 u=f~~in~~\Omega$$ $$u=g_1~~on ~~\partial \Omega$$ $$\frac {\partial u}{\partial n}=g_2~~on ~~\partial \Omega$$

Using Green's formula, we have that if $u$ is smooth enough $$\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$$

the test space is $H^2(\Omega)$, How could I use the given boundary conditions$g_1,g_2$ .

how could I find the new terms $\frac {\partial v }{\partial n}$ and $\frac {\partial \Delta u}{\partial n}$

Answer 1

The biharmonic equation is the Euler-Lagrange equation of the Laplacian energy $\frac{1}{2} \langle \Delta u,\Delta u \rangle$. A systematic approach to discretize higher order problems is to convert the unconstrained problem to a constrained problem: Minimize $\frac{1}{2} \langle v,v\rangle$ s.t. $\Delta u=v$; that is, \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} where we used Green's identity for the second equation and $\frac{\partial }{\partial n}$ is the normal derivative.

The new optimization problem can be discretized using piecewise linear elements as is commonly done for second order problems. This is the common thread in mixed finite-element discretizations.

Answer 2

I decided to expand my earlier comment into an answer. I'd suggest using the Morley element that uses $P_2$ basis in each element and the degrees of freedom are

1. values at vertices
2. normal derivatives at edge midpoints

It is probably the easiest element for biharmonic problem implementation wise.

Here is an example using sp.fem library (version 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()    

The solution looks like this:

Note that both $g_1$ and $g_2$ are easily defined as Dirichlet-type boundary conditions.

Answer 3

Using finite element methods, you also have the choice of discretizing the weak formulation directly, using two types of approaches:

1. using H²-conforming elements like the Argyris element or the Hsieh-Clough-Tocher on triangles or the Bogner-Fox-Schmidt element on quadrilaterals.

2. Using regular H¹-conforming elements in conjunction with the so called C⁰-interior penalty method.

If you take the weak formulation you derived through Green's formula, choosing test spaces with $v = \partial_n v = 0$ on the boundary makes the boundary terms go away.

Another approach very well suited to high order differential equations is isogeometric analysis.

Answer 4

The variational or weak formulation of the biharmonic problem reads as : find $u \in V_D$ such that $$ a(u,v) = \ell(v), \quad \forall v \in V_0,$$ where the bilinear and linear forms are given by $$a(u,v) = \int_\Omega \Delta u \Delta v \,d\Omega \quad \text{and} \quad \ell(v) = \int_\Omega f v\, d\Omega,$$ and the hyperplane and test spaces are given by $V_D :=\{ v \in H^2(\Omega), v = g_1, \frac{\partial v}{\partial n} = g_2 \quad on \quad \partial \Omega\}$ and $V_0 :=\{ v \in H^2(\Omega), v = 0, \frac{\partial v}{\partial n} = 0 \quad on \quad \partial \Omega\}.$

The higher order derivatives on the $\partial \Omega$ are actually the Neumann boundary conditions. These are not given hence your choice of test function is enough. This is a biharmonic problem with Dirichlet boundary conditions.