Solution does not fulfil homogeneous Neumann conditions

[MPenaR]

First of all sorry for the long question. I could not find a smaller MWE. If someone reads it I would be really appreciate it. I tried to explain it as much as I could.

The Physical problem

Let \mathcal{R} denote the rectangular domain:

\mathcal{R}:=\left\lbrace\left(x,y\right)\in\mathbb{R}^2, \; 0<x<L_x, 0<y<L_y\right\rbrace

whose boundary will be split in:

\Gamma_\mathrm{front}:= \left\lbrace \left(0,y\right),\,0<y<L_y\right\rbrace

\Gamma_\mathrm{back}:= \left\lbrace \left(L_x,y\right),\,0<y<L_y\right\rbrace

\Gamma_\mathrm{sides}:= \left\lbrace \left(x,0\right),\,0<x<L_x\right\rbrace\cup\left\lbrace \left(x,L_y\right),\,0<x<L_x\right\rbrace

I’m trying to solve the heat equation in time-harmonic regime:

\begin{cases} -\kappa\Delta T-i\omega\rho cT=0 & \mathcal{R}\\ -\kappa\nabla T\cdot\mathbf{n}=0 & \Gamma_{\mathrm{sides}}\\ -\kappa\nabla T\cdot\mathbf{n}=h_{\mathrm{eff}}T-\alpha Q & \Gamma_{\mathrm{front}}\\ -\kappa\nabla T\cdot\mathbf{n}=h_{\mathrm{eff}}T & \Gamma_{\mathrm{back}} \end{cases}

where \kappa, \omega,\rho,c,\alpha and h_\mathrm{eff} are positive constants, and

Q\left(\mathbf{x}\right) = \frac{I}{2\pi}\frac{\left(\mathbf{x}-\mathbf{s}\right)\cdot\mathbf{i}}{\left\Vert\mathbf{x}-\mathbf{s}\right\Vert^2}

where I is a complex constant, \mathbf{i} is the unit vector in the x direction, and \mathbf{s}\in\mathbb{R}^2 such that \mathbf{s}\cdot\mathbf{i}<0.

I construct the complex weak formulation, and I get that T\in H^1\left(\mathcal{R}\right) must satisfy:

a\left(T,\Psi\right)=L\left(\Psi\right)\quad\forall\Psi\in H^1\left(\mathcal{R}\right)

with

a\left(\Phi,\Psi\right)=\int_{\mathcal{R}}\kappa\nabla\Phi\cdot\overline{\nabla\Psi}\,\mathrm{d}\mathbf{x}-\int_{\mathcal{R}}i\omega\rho c\Phi\overline{\Psi}\,\mathrm{d}\mathbf{x}+\int_{\Gamma_{\mathrm{front}}}h_{\mathrm{eff}}\Phi\overline{\Psi}\,\mathrm{d}s+\int_{\Gamma_{\mathrm{back}}}h_{\mathrm{eff}}\Phi\overline{\Psi}\,\mathrm{d}s

L\left(\Psi\right)=\int_{\Gamma_\mathrm{front}} \alpha Q \overline{\Psi}\,\mathrm{d}s

where the \overline{\Psi} stands for complex conjugation.

Dolfin solution

To implement that in dolfin I’m splitting the complex variational form in its real and imaginary parts. Defining \Phi=\Phi_R + i\Phi_I and \Psi=\Psi_R + i\Psi_I the forms are rewritten as:

\Re\left(a\left(\left(\Phi_{R},\Phi_{I}\right),\left(\Psi_{R},\Psi_{I}\right)\right)\right)=\\ \int_{\mathcal{R}}\kappa\left(\nabla\Phi_{R}\cdot\nabla\Psi_{R}+\nabla\Phi_{I}\cdot\nabla\Psi_{I}\right)\,\mathrm{d}\mathbf{x}+\\ \int_{\mathcal{R}}\omega\rho c\left(-\Phi_{R}\Psi_{I}+\Phi_{I}\Psi_{R}\right)\,\mathrm{d}\mathbf{x}+\\ \int_{\Gamma_{\mathrm{front}}}h_{\mathrm{eff}}\left(\Phi_{R}\Psi_{R}+\Phi_{I}\Psi_{I}\right)\,\mathrm{d}s+\int_{\Gamma_{\mathrm{back}}}h_{\mathrm{eff}}\left(\Phi_{R}\Psi_{R}+\Phi_{I}\Psi_{I}\right)\,\mathrm{d}s

\Re\left(L\left(\Psi_{R},\Psi_{I}\right)\right)=\int_{\Gamma_{\mathrm{front}}}\alpha\left(Q_{R}\Psi_{R}+Q_{I}\Psi_{I}\right)\,\mathrm{d}s

where Q_R and Q_I are respectively the real and imaginary parts of Q.

The code for implementing it in dolfin is:

#imports
from dolfin import  RectangleMesh, MeshFunction, SubDomain, Measure, FunctionSpace, MixedFunctionSpace, Constant
from dolfin import TrialFunctions, TestFunctions, inner, grad, dx, Function, solve, plot, Point, near, Expression


#geometrical constants and mesh generation
Lx = 0.01
aspect_ratio = 100
Ly = aspect_ratio*Lx
Nx = 5
Ny = aspect_ratio*Nx 

mesh = RectangleMesh( Point(0,0), Point(Lx,Ly), Nx, Ny )

boundary_markers = MeshFunction("size_t", mesh, 1, 0)

tol = 1E-14

class Gamma_front(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, tol)

GF = Gamma_front()
GF.mark(boundary_markers, 0)

class Gamma_back(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], Lx, tol)

GB = Gamma_back()
GB.mark(boundary_markers, 1)

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


#Physical constants
kappa = 200
h_conv = 15 
epsilon = 0.8
sigma = 5.67E-8
T_air = 298
rho = 2700
c_sp = 900
omega = 0.8
k = Constant(kappa)
wrc = Constant(omega*rho*c_sp)
h_eff = Constant(h_conv + 4*epsilon*sigma*T_air**3)

alpha = 0.4

I = 600*(1 + 0.0j)

sx = -15*Lx
sy = 0.5*Ly

Q_R = Expression("alpha*I_R/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_R = I.real, sx=sx, sy=sy)
Q_I = Expression("alpha*I_I/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_I = I.imag, sx=sx, sy=sy)


# Funtion spaces 
P1 = FunctionSpace(mesh, "Lagrange", 1 )
H = MixedFunctionSpace( P1, P1  ) 
Phi_R, Phi_I = TrialFunctions(H)
Psi_R, Psi_I = TestFunctions(H)


# variational forms
a =  k*inner( grad(Phi_R), grad(Psi_R) )*dx + k*inner( grad(Phi_I), grad(Psi_I) )*dx  \
                                          - wrc*Phi_R*Psi_I*dx +  wrc*Phi_I*Psi_R*dx  \
                                          + h_eff*Phi_R*Psi_R*ds(0) + h_eff*Phi_I*Psi_I*ds(0) \
                                          + h_eff*Phi_R*Psi_R*ds(1) + h_eff*Phi_I*Psi_I*ds(1)
L = Q_R*Psi_R*ds(0) + Q_I*Psi_I*ds(0)


#problem solution
T = Function(H)
solve(a==L, T)
T_R, T_I = T.split()

The question

The code runs without any problem, however when I plot the solution T\left(x_i,y\right) for several values of x_i:

#plotting of the solution
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,Lx,Nx)
for x_i in x:
    y = np.linspace(0,Ly,Ny)
    T = np.array([T_R(x_i,y_i) for y_i in y])
    plt.plot(y,T)
plt.xlabel('$y$ [m]')
plt.legend([f'$T({100*x_i:.2f} \\times 10^{{-2}},y)$' for x_i in x])

I realize that the homogeneous Neumann conditions are not being met, that is, the slope both at y=0 and y=L_y should be 0.

Am I missing something?

[conpierce8]

You should choose different markers for the front and back of your domain, e.g. 1 and 2. In your current code, both the front and sides of your domain are marked as 0 (because the default value of boundary_markers is 0), so the homogeneous Neumann condition is not enforced anywhere.

[MPenaR]

Thanks: