System of two equations: setting both Neumann and Dirichlet condition at the same boundary

variational formulation
1

Hello,

I am relatively new to FEniCS and finite elements in general, and I just came across an equation I do not know how to solve. This issue is that I do not understand how to set the boundary conditions given the weak forms of the equations.

I have a system of two PDEs, time-dependant and in 1D space (spatial coordinate r), I am solving for h and p, with associated test functions u and v. One of my equation is something like p \sim \partial^2{h}/\partial{r^2}. When writing the weak formulation for this equation, integration by parts leads to a boundary integral of the form \int_{\partial \Omega} v ~\partial{h}/\partial{r}~{\rm d}s.

So far so good. But it seems like one of my boundary conditions can´t be applied through that integral. I set Neumann conditions for h at r=0 and r=A (the boundaries of my 1D domain), a Neumann condition for p at r=0, and a Dirichlet condition for p at r=A.

Going back to this term: \int_{\partial \Omega} v ~\partial{h}/\partial{r}~{\rm d}s
At r=A from the Dirichlet condition imposed on p there, the test function v is 0 there regardless of the value I want to impose on \partial{h}/\partial{r}: one of my Neumann cannot be imposed. When trying to implement this in FEniCS I can see that this boundary condition is not fulfilled, but this is not surprising.

Is there something I misunderstand, or is this a general problem when mixing Dirichlet and Neumann conditions? If so, how is it usually solved?

Thank you.

P.S. I think this is a general/conceptual problem but in case this matters, here is my specific system of equation:
\partial{h}/\partial{t}=\frac{1}{12 r}\frac{\partial}{\partial{r}}\left(r h^3 \partial{p}/\partial{r}\right)
p=\frac{1}{\delta} \left(2-\frac{1}{r} \frac{\partial}{\partial{r}}\left(r\partial h /\partial r\right)\right)
for r\in\Omega=[0,A] with A and \delta some constants, and with boundary conditions:
\partial{p}/\partial{r}= \partial{h}/ \partial{r}= 0 at r=0
\partial{h}/ \partial{r} = B (a constant) at r=A
p=0 at r=A
and some initial condition for h(t=0,r)
I find the weak form by multiply the equations by r\times{\rm(test~function)} and integrating over \Omega.

2

Someone suggested me to use a third intermediate function: introducing q=\partial h/ \partial r and solving for the system of 3 PDEs allows to have the right boundary conditions.

It works well but it seems a bit magical to me that introducing an intermediate function changes the behaviour so much; I can see it in the weak form of the equations; but if anyone has insights to understand more deeply why this works that would still be welcomed.

1 Like
3

Doesn’t a Robin Boundary Condition work for this?

u = g \quad \text{on } \Gamma \\ \nabla u \cdot \vec{n} = h \quad \text{on } \Gamma \\

Can be combined as

g \nabla u \cdot \vec{n} = uh \quad \text{on } \Gamma

Which is a specific form of a Robin Boundary Condition

\nabla u \cdot \vec{n} = \beta - ru \quad \text{on } \Gamma