Issues with Biharmonic Equation Decomposition in FEniCS

variational formulation
1

Hello, everyone. I’m encountering some issues while learning FEniCS, specifically focused on solving the biharmonic equation.

To be more specific, I studied the official biharmonic equation example Biharmonic equation — DOLFIN documentation, but it doesn’t fully match the problem I want to solve, so I need to modify it. I thought of testing the problem using a decomposition method commonly used in engineering.

The original problem is: ΔΔu = 0 in Ω, with boundary conditions u = 0 on ∂Ω and Δu = 0 on ∂Ω.
I attempted to decompose it into: Δu = v in Ω, Δv = 0 in Ω, with boundary conditions u = 0 on ∂Ω and v = 0 on ∂Ω (derived from Δu = 0 on ∂Ω).

However, in this case, since Δv = 0 in Ω and v = 0 on ∂Ω, v is identically zero throughout the domain, reducing the biharmonic equation to a Laplace equation. But when I set the right-hand term ΔΔu = f in Ω with f = 0 in the biharmonic equation demo, the computed results are clearly different from the Laplace equation solution. Where did my analysis go wrong? I would be extremely grateful if anyone could point out the flaw in my reasoning.

By the way, if you have the time, could you also take a look at my other post How to impose Neumann BC for Biharmonic equation? This one is about boundary conditions and choosing the penalty parameter for the biharmonic equation, and any insights there would be incredibly helpful.

Best regards.

2

At first glance, it seems to me that your analysis is correct. The example code uses zero dirichlet on u, though, so I guess you had to implement non-homogeneous BC’s to compare with the pure Poisson equation? Could it be that a mistake slipped in there? Could you share the code snippets that show the unexpected different results.

3

I hope you are familiar with some concepts of mesh deformation.

I am trying to solve the Laplace or biharmonic equation for mesh deformation. The mesh is an NACA0012 airfoil, with the deformation prescribed as a rotation about the quarter-chord point at (-0.25, 0). In the example, only the displacement contour in the X-direction after rotation is shown.

It might seem a bit complex, but it is simplified to a simple mathematical expression:
Laplace: Δu = 0 in Ω, u = fun_rotate on ∂Ω_in, and u = 0 on ∂Ω_out.
Biharmonic: ΔΔu = 0 in Ω, u = fun_rotate on ∂Ω_in, u = 0 on ∂Ω_out, and Δu = 0 on ∂Ω_in and ∂Ω_out.

I am currently unsure whether the result for the biharmonic equation is correct, as according to my previous analysis, it should be identical to the Laplace result.
So either the biharmonic result is incorrect, or my previous analysis is flawed.

If you could provide any assistance, I would be very grateful.
My code is:

from dolfin import *
import ufl
import matplotlib.pyplot as plt

file = 'naca0012'
x0 = -0.25
y0 = 0.0
angle0 = 0.2

def run_fenics_biharm():
    parameters["form_compiler"]["cpp_optimize"] = True
    parameters["form_compiler"]["optimize"] = True
    parameters["ghost_mode"] = "shared_facet"
    grid = Mesh((file + ".xml"))
    boundaries = MeshFunction("size_t", grid, (file + "_facet_region.xml"))

    h = CellDiameter(grid)
    h_avg = (h('+') + h('-')) / 2.0
    n = FacetNormal(grid)
    f = Constant(0.0)
    alpha = Constant(8.0)

    V1 = FunctionSpace(grid, 'CG', 2)
    u_D1 = Expression('x[0] * cos(angle) + x[1] * sin(angle) + x0 * (1 - cos(angle)) - y0 * sin(angle) - x[0]',
                      degree=2, x0=x0, y0=y0, angle=angle0)
    bc_in1 = DirichletBC(V1, u_D1, boundaries, 2)
    bc_out1 = DirichletBC(V1, Constant(0.0), boundaries, 3)
    bc1 = [bc_in1, bc_out1]
    u1 = TrialFunction(V1)
    v1 = TestFunction(V1)
    try:
        a1 = (inner(div(grad(u1)), div(grad(v1))) * ufl.dx
              - inner(avg(div(grad(u1))), jump(grad(v1), n)) * dS
              - inner(jump(grad(u1), n), avg(div(grad(v1)))) * dS
              + alpha / h_avg * inner(jump(grad(u1), n), jump(grad(v1), n)) * dS)
        L1 = f * v1 * ufl.dx
    except TypeError as e:
        print(f"\nCaught TypeError: {e}")
        a1 = (inner(div(grad(u1)), div(grad(v1))) * dx
              - inner(avg(div(grad(u1))), jump(grad(v1), n)) * dS
              - inner(jump(grad(u1), n), avg(div(grad(v1)))) * dS
              + alpha / h_avg * inner(jump(grad(u1), n), jump(grad(v1), n)) * dS)
        L1 = f * v1 * dx
    u1 = Function(V1)
    solve(a1 == L1, u1, bc1)
    return u1


def run_fenics():
    parameters["form_compiler"]["cpp_optimize"] = True
    parameters["form_compiler"]["optimize"] = True
    parameters["ghost_mode"] = "shared_facet"
    grid = Mesh((file + ".xml"))
    boundaries = MeshFunction("size_t", grid, (file + "_facet_region.xml"))

    V1 = FunctionSpace(grid, 'CG', 2)
    u_D1 = Expression('x[0] * cos(angle) + x[1] * sin(angle) + x0 * (1 - cos(angle)) - y0 * sin(angle) - x[0]',
                      degree=2, x0=x0, y0=y0, angle=angle0)
    bc_in1 = DirichletBC(V1, u_D1, boundaries, 2)
    bc_out1 = DirichletBC(V1, Constant(0.0), boundaries, 3)
    bc1 = [bc_in1, bc_out1]
    u1 = Function(V1)
    v1 = TestFunction(V1)
    f1 = Constant(0.0)
    try:
        F1 = dot(grad(u1), grad(v1)) * ufl.dx - f1 * v1 * ufl.dx
    except TypeError as e:
        print(f"\nCaught TypeError: {e}")
        F1 = dot(grad(u1), grad(v1)) * dx - f1 * v1 * dx
    solve(F1 == 0, u1, bc1)
    return u1


u_laplace = run_fenics()
u_biharm = run_fenics_biharm()
plot(u_laplace)
plt.show()
plot(u_biharm)
plt.show()

My mesh file is: