How to impose Neumann BC for Biharmonic equation

variational formulation
1

Hello, everyone

I am a newcomer to FEniCS, and I am working hard to learn the relevant tutorials, but it is still somewhat challenging for me, especially regarding the demos related to the Biharmonic equation: Biharmonic equation — DOLFIN documentation (fenicsproject.org)

My doubts mainly focus on whether the boundary conditions for the equation are sufficient. In this tutorial, it points out that boundary conditions u=0 on ∂Ω and ∇2u=0 on ∂Ω are applied, but I did not find in the code how the boundary condition ∇2u=0 on ∂Ω is implemented. Could you please guide me on this?

At the same time, if I want to apply a boundary condition like n⋅∇u=0 to the Biharmonic equation, how should I do that?

Finally, I am curious about how the penalty parameter is chosen; is the value of 8 suitable for any situation?

I noticed that there are some tutorials in the demo for correctly applying Neumann boundary conditions, but at my current level, it may be difficult for me to transfer this knowledge to solving the Biharmonic equation. Therefore, can someone teach me how to apply boundary conditions such as n⋅∇u=0 or ∇2u=0 for the Biharmonic equation? I would be very grateful.

Best Regards.

2

If you compute the integration by parts of

\begin{align} \int_\Omega \nabla^4(u):v~\mathrm{d}x = -\int_{\Omega} \nabla (\nabla^2 u): \nabla v~\mathrm{d}x + \int_{\partial\Omega} \nabla (\nabla^2u) v~\mathrm{d} x = -\int_{\Omega} \nabla (\nabla^2 u): \nabla v~\mathrm{d}x \end{align}

when v\in V^h_0.

Then you can integrate this inner product by parts using (Integration by parts - Wikipedia)

\begin{align} -\int_{\Omega} \nabla (\nabla^2 u): \nabla v~\mathrm{d}x &= \int_\Omega \nabla^2u: \nabla \cdot (\nabla v)~\mathrm{d}x - \int_{\partial\Omega} \nabla^2u \cdot \nabla v \cdot n~\mathrm{d}s\\ &= \int_\Omega \nabla^2u: \nabla^2 v~\mathrm{d}x - \int_{\partial\Omega} \nabla^2u \cdot \nabla v \cdot n~\mathrm{d}s \end{align}

Thus you use the natural boundary condition \nabla^2u=0 to get the variational form

\int_\Omega \nabla^4(u):v~\mathrm{d}x =\int_\Omega \nabla^2u: \nabla^2 v~\mathrm{d}x

Note that in my derivation I have skipped the weak enforcement internally, strictly for clarity. See: Biharmonic equation — DOLFINx 0.8.0 documentation for the demo in DOLFINx

3

Thank you for your formula derivation; it’s very clear, and even I can handle it without any pressure. It perfectly resolves my confusion about boundary conditions, so I really appreciate it.

However, I still have questions regarding boundary conditions for Biharmonic equation. Specifically, when the Neumann boundary condition is first-order, do you have any ideas on how to address it? For instance, with the natural boundary condition n⋅∇u=0, how should I approach this? I’m currently at a loss.

By the way, the problem I want to solve can be relatively more complex. The computational domain has two boundaries, which you can imagine as resembling a ring. The Dirichlet boundary conditions on the inner and outer boundaries are known. I want to satisfy n⋅∇u=0 on the inner boundary and ∇2u=0 on the outer boundary. For these slightly more complex boundary conditions, I might have already imagined the complexities of its weak form for Biharmonic equation. I’ll keep learning and hope to master it soon.

Best wishes!

4

If I am not mistaken, these should be treated similarly to the jump terms of \nabla u over internal facets in the DG weak formulation except that the jump terms is now the difference of the solution gradient and the boundary condition gradient. See for instance: Clamped Kirchhoff-Love plate — fenics-shells

1 Like
5

Hello, and thanks again for your derivation! After diving into it more deeply, I have some new questions, especially regarding the choice of penalty parameter in solving the biharmonic equation. I’m referring specifically to https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/biharmonic/demo_biharmonic.py.html and Clamped Kirchhoff-Love plate — fenics-shells.

In the first case, α is set to 8.0, which feels a bit arbitrary to me. Why was it chosen this way? In the second case, α is set to (E*t)^3, which seems to have a stronger justification. But if we disregard the physical interpretation and look at it purely as a mathematical problem, it seems challenging to determine a suitable α. Are there any general guidelines for selecting a penalty parameter?

Moreover, I’m a bit confused. Shouldn’t the problem we’re solving have a unique solution? At least, that was my impression. But when I change α in these examples, the solution changes. How should I understand this apparent non-uniqueness? How do I ensure that the solution I obtain is correct (i.e., that my choice of α is appropriate)? Or, if I’m using FEniCS to solve my problem, how can I be sure that the solution is accurate, and how close it is to the true analytical solution? It seems difficult to judge.

Any suggestions would be greatly appreciated. Thanks a lot! And if anyone else can chime in, I’d be just as grateful.

6

Hello, I have studied the Kirchhoff-Love example in detail, and it’s really great—thanks for recommending it!

However, new questions have arisen. I noticed that there seem to be significant differences between the Kirchhoff-Love example and the biharmonic equation example, in various aspects.

Specifically, the Kirchhoff-Love example consists of Dirichlet boundary conditions combined with ∂w/∂n = 0, using a weak form of boundary conditions, whereas the biharmonic equation example uses Dirichlet boundary conditions with Δu = 0, which is essentially a natural boundary condition. The solution methods for the two examples also differ significantly.

In fact, the problem I want to solve is a biharmonic equation with given Dirichlet boundary conditions, where part of the boundary should satisfy Δu = 0 and another part should satisfy ∂u/∂n = 0. How should I approach this? Are there any examples I could refer to?

Alternatively, should I modify the Kirchhoff-Love example code (to add Δu = 0) or the biharmonic equation example code (to add ∂u/∂n = 0) to achieve my goal? The solution methods seem quite different.

Any suggestions would be greatly appreciated. And if anyone else can help, I’d be just as grateful. Thanks a lot!

7

Hello,
\Delta u=0 is a natural boundary condition for the biharmonic equation. It is therefore accounted for automatically in the weak form when ignoring the corresponding term on the outer edges. For \partial u/\partial n =0, it is an essential boundary condition but which cannot be enforced directly as a dirichletbc since we use a DG scheme. Thus you should include it in the weak form exactly as handled as the L_{BC} term in the Kirchhoff demo. Finally, please note that the biharmonic demo requires 2 different BCs for a given point on the boundary. Here you seem to be missing some info since you only enforce either \Delta u=0 or \partial u/\partial n=0 at a given point. Do you also have u=0 on the whole boundary ?

8

Thank you for your prompt reply. I haven’t forgotten to impose the Dirichlet boundary condition u = 0 on the entire boundary, thanks for the reminder.

So, if I understand correctly, in the biharmonic equation demo, the boundary conditions are Dirichlet boundary conditions + natural boundary condition Δu = 0, while in the Kirchhoff-Love demo, the boundary conditions are Dirichlet boundary conditions + natural boundary condition Δu = 0 + additional ∂u/∂n = 0? Is there any misunderstanding on my part?

For my problem, I want ∂u/∂n = 0 to hold only on part of the boundary, so I need to modify the following code to apply the weak form boundary condition on just part of the boundary instead of the whole boundary, as in the original code. Is this correct?

facet_function = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
facet_function.set_all(0)
all_boundary.mark(facet_function, 1)
ds = Measure("ds")(subdomain_data=facet_function)

By the way, would you be able to help me with another question regarding the penalty parameter?

Thanks again for your help, best wishes!

9

Hello everyone, just bumping this thread.

Any suggestions would be greatly appreciated

10

Hello, sorry to bother you again.

Due to my limited knowledge of finite element methods, I’m facing significant challenges using FEniCS to solve my problem, even though I have the excellent Clamped Kirchhoff-Love plate demo as a reference.

In simple terms, my problem is straightforward: I need to solve a biharmonic equation with fully known Dirichlet boundary conditions, where part of the boundary has a Neumann boundary condition Δu = 0, and another part has a Neumann condition ∇u ⋅ n = 0. In this case, the equation should have a unique solution.

However, in the biharmonic equation demo, I cannot apply the ∇u ⋅ n = 0 boundary condition, and in the Clamped Kirchhoff-Love plate demo, the problem becomes more complicated due to the external forces, whereas I want to solve a purely mathematical problem without involving physical parameters like elastic modulus, thickness, or Poisson’s ratio.

For days, I’ve been feeling lost, so I am reaching out to you and others for help. How should I approach solving this problem?

I’ve been trying to combine the two demos by analogy to achieve my goal, but it hasn’t been successful so far. I think I need help.

Thanks again to all of you for your kindness.

11

The strong condition \nabla v\cdot n = g is quite tricky in FEM, for typical elements. These do not have normal gradients as degrees of freedom, so setting them manually is difficult. You have two options:

  • You use a mixed form of the PDE. You already imply this is something you’re considering in your other post Issues with Biharmonic Equation Decomposition in FEniCS If you do that, the condition actually becomes a natural condition, and merely requires integration over the boundary. I believe that the old fenics versions had some demos on that. Those might help you on your way.
  • You leverage the ‘DG’-type technology that is already in the biharmonic example. In that example, these terms are added that integrate over the internal facets:
  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
  + alpha/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS

Those terms effectively enforce continuity of \nabla v\cdot n from element to element (i.e, (\nabla v\cdot n)^{left} = -(\nabla v\cdot n)^{right}. You can do the same on the domain boundary, to enforce \nabla v\cdot n = 0 there:

  - inner(div(grad(u)), dot(grad(v), n))*ds \
  - inner(dot(grad(u), n), div(grad(v)))*ds \
  + alpha/h_avg*inner(dot(grad(u),n), dot(grad(v),n))*ds

Note the small letter ds as opposed to dS, changing the domain of integration from interior facets to domain boundary.

That is called Nitsche’s method, which may give you a handle on further search.

1 Like
12

Thank you for your analysis. However, it seems a bit complex for me, and I will study it in more detail. If I have further questions, I will ask again in the future.

That said, I am still not entirely clear on how to choose the penalty parameter. In this demo Biharmonic equation — DOLFIN documentation, setting it to 8 feels somewhat arbitrary to me.

In my own problem, I am unsure about the proper selection of the penalty parameter. I tried different values, and when it was too small, the solution seemed unstable. However, with larger values, the results varied, which contradicts my understanding of the uniqueness of the finite element solution.

13

The value of 8 might seem arbitrary, but has quite some mathematical basis. In the end it is still an estimate though, and hence a nice round number. To prove stability you require a minimal value. But, choosing it too large unnecessarily constraints the solution space, leading to poor solution behavior. You’ve observed both of these results.

So then, what value to choose? Mathematically, the penalty must be large enough to ensure positive definiteness of the resulting matrix. More specifically, it must be large enough to ensure the penalty boundary integral alpha/h_avg*inner(dot(grad(u),n), dot(grad(v),n))*ds always produces larger values than the other boundary integrals - inner(div(grad(u)), dot(grad(v), n))*ds - inner(dot(grad(u), n), div(grad(v)))*ds for equalling u and v. The value that ensures this follows from an ‘inverse estimate’; a specially constructed matrix eigenvalue problem. You can find the proof in Lemma 4.3 in this work: Redirecting That is certainly not the best reference to familiarize yourself with this theory, I just know the proof is in there.

Useful references for actually determining/estimating what that value should then be are: Redirecting , Redirecting.

1 Like