How to implement inhomogeneous Neumann boundary condition in FeniCSX?

dolfinx
1

Hi dear FeniCSX users,

I am using FeniCSX to solve the following 1D coupled system:

\partial_{\tilde{t}} \tilde{m}= \mathrm{Lu} \partial_{xx} \tilde{m} - \mathrm{Lu}\,\mathrm{Pn} \partial_{xx} \tilde{T}
\partial_{\tilde{t}} \tilde{T}= \partial_{xx} \tilde{T} - \varepsilon\mathrm{Ko} \partial_{\tilde{t}} \tilde{m}

with the boundary conditions

\partial_{\tilde{x}} \tilde{m} (0, \tilde{t}) - Pn \partial_{\tilde{x}} \tilde{T} (0, \tilde{t}) = 0
\partial_{\tilde{x}} \tilde{T} (0, \tilde{t}) = -Q

- \partial_{\tilde{x}} \tilde{m} ( 1, \tilde{t}) + Pn \partial_{\tilde{x}} \tilde{T} (1, \tilde{t}) + B_{im} (1 - \tilde{m}) = 0

\partial_{\tilde{x}} \tilde{T} ( 1, \tilde{t})= B_{it} \bigl(1 - \tilde{T}(1, \tilde{t}) \bigr) + (1 - \varepsilon) Ko Lu B_{im} [ \tilde{m} (1, \tilde{t}) - 1]

and the initial conditions

\tilde{m}(\tilde{x}, ,0)=0
\tilde{T}(\tilde{x} ,0)=0

The weak formulation of this problem is given by:

F0 = \int \partial_{\tilde{t}} \tilde{m} q + Lu \int \partial_{\tilde{x}} \tilde{m} \partial_{\tilde{x}} q - Lu Pn \int \partial_{\tilde{x}} \tilde{T} \partial_{\tilde{x}} q + Lu \int_{x=1} B_{im} [m-1] q

F1 = \int \partial_{\tilde{t}} \tilde{T} v + \int \partial_{\tilde{x}} \tilde{T} \partial_{\tilde{x}} v - \varepsilon Ko \int \partial_{\tilde{t}} \tilde{m} v + \int_{x=0} Q v + \int_{x=1} B_{it} \bigl( \tilde{T}(1, \tilde{t}) - 1 \bigr) - (1 - \varepsilon) Ko Lu B_{im} [ \tilde{m} (1, \tilde{t}) - 1] v

Below is my MWE with FeniCSX:
First attempt:


F0 =  inner(m, q) * dx - inner(m_n, q) * dx + dt * Lu_true * inner(grad(m), grad(q)) * dx - dt * ( Lu_true * Pn_true ) * inner(grad(T), grad(q)) * dx

Term_BCs_m = Lu_true * dt * inner(Bim_true*(m-1), q) * ds(2)    

F0_ =  F0 + Term_BCs_m

F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt * inner(grad(T), grad(v)) * dx + (epsilon_true * Ko_true) * inner(m, v) * dx - (epsilon_true * Ko_true) * inner(m_n, v) * dx 

Term_BCs_T_c =  dt * inner(Q, v) * ds(1)  
Term_BCs_T_d =  dt * inner(Bit_true*(T-1) - (1-epsilon_true)*Ko_true*Lu_true*Bim_true(m-1), v) * ds(2)

F1_ = F1 + Term_BCs_T_c + Term_BCs_T_d

Second attempt:

n = FacetNormal(domain)
Q = -dot(n, grad(T))

F0 =  inner(m, q) * dx - inner(m_n, q) * dx + dt * Lu_true * inner(grad(m), grad(q)) * dx - dt * ( Lu_true * Pn_true ) * inner(grad(T), grad(q)) * dx

Term_BCs_m = Lu_true * dt * inner(Bim_true*(m-1), q) * ds(2)    

F0_ =  F0 + Term_BCs_m

F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt * inner(grad(T), grad(v)) * dx + (epsilon_true * Ko_true) * inner(m, v) * dx - (epsilon_true * Ko_true) * inner(m_n, v) * dx 

Term_BCs_T_c =  dt * inner(-Q, v) * ds(1)  
Term_BCs_T_d =  dt * inner(Bit_true*(T-1) - (1-epsilon_true)*Ko_true*Lu_true*Bim_true(m-1), v) * ds(2)

F1_ = F1 + Term_BCs_T_c + Term_BCs_T_d

Both gave the same result. But the temperature should be between [-4, 2] but in my exceeds 2 and equal to 13.

I am struggling to understand how to define the inhomogeneous Neumann boundary condition in my case.

Any help would be appreciate.

Thank you.