How to implement inhomogeneous Neumann boundary conditions for a system of PDEs (Coupled)?

variational formulation
1

Hi everyone,

I am able to implement Dirichlet boundary conditions for coupled PDEs from the tutorial : https://fenicsproject.org/pub/tutorial/html/._ftut1010.html
However, the tutorial assumes homogeneous Neumann boundary conditions however for my couple expression:
Equations
I have input flux for the first equation and zero flux for the second expression. I’m not sure if my variational expression applies a flux to both u1 and u2:

#Define Boundary Conditions
Flux = Expression("t<30+1e-9 ? ((k1*A1)/Z_w)*exp(-1*k1*t) : 0", degree=2, k1=k1, A1=A1, Z_w=Z_w, t=0)
u1_D = Constant(0)

boundaryconditions = {
    1 : {"Neumann" : Flux},
    2 : {"Dirichlet" : u1_D}
}

ds = Measure('ds', domain=mesh, subdomain_data=mf)

bcs = []
for i in boundaryconditions:
    if "Dirichlet" in boundaryconditions[i]:
        # bc = DirichletBC(V, boundaryconditions[i]["Dirichlet"], boundaries, i)
        bc = DirichletBC(V.sub(0), boundaryconditions[i]["Dirichlet"], boundaries, i) # Dirichlet is appied only to 1 subspace
        bcs.append(bc)

integrals_N = []
for i in boundaryconditions:
    if "Neumann" in boundaryconditions[i]:
        if boundaryconditions[i]["Neumann"] != 0: #if not zero flux
            g = boundaryconditions[i]["Neumann"]
            integrals_N.append(v1*g*ds(i))

#Define Variational Problem
F = (((u1-u_n1 + u2-u_n2)*v1)/dt)*dx + D_T*dot(grad(u1),grad(v1))*dx + (((u2-u_n2)*v2)/dt)*dx - dt*sum(integrals_N) \
    - v2*ka*u1*B_M*dx + ka*B_M*v2*u2*dx - kd*v2*u2*dx

This does not generate any error but the flux seems to be implemented to both u1 and u2.
Is this the right way to implement flux for just u1 and not u2?
Any help/advice will be greatly appreciated

Thanks

2

I think you are right in the variational form, I assume you are using mixed fem, with u2 is zero flux, you write Constant(0)v2ds, or just don’t write anything

3

Hi, thanks for getting back to me. Yup I’m using mixed fem through the following:

    element = MixedElement([P1,P1])
    V = FunctionSpace(mesh, element)
    v1,v2 = TestFunction(V)
    u = Function(V)
    u_n = Function(V)
    u1,u2 = split(u)
    u_n1,u_n2 = split(u_n)

And I’ve elected to not to write anything to implement a zero flux for u2.