Heat balance equations with Neuman conditions on the adjacent boundary (one-dimensional problem)

variational formulation
1

Hi everyone,

I am trying to solve a system of two heat balance equations:

\frac{1}{A} \nabla^2 u_1(x) - u_1(x) = - (T_b + p),

in one dimensional domain [0, x_0], where A = 2000 is a constant, T_b = 2.0 is a constant, p is the unknown heat source, and u_1 is the temperature distribution over the domain [0, x_0].

\frac{1}{A} \nabla^2 u_2(x) - u_2(x) = - T_b,

in one dimensional domain [x_0, L], where L is the length (a constant), and u_2 is temperature distribution over the domain [x_0, L]. In both cases x_0 is unknown and should be found from the second boundary condition.

The boundary conditions are:

  • at x = 0, \nabla u_1(0) = 0;
  • at x = x_0, \nabla u_1(x_0) = \nabla u_2(x_0), and u_1(x_0) = u_2(x_0) = T_c, where T_c = 10 is a constant;
  • at x = L, u_2(L) = T_b.

The problem here is that u_1, u_2, and p depend on the boundary x_0 which is not defined.

Alternatively, the problem can be treated as for subdomains made of two different materials:

\frac{1}{A} \nabla^2 u(x) - u(x) = - (T_b + p),

in one dimensional domain [0, L], where p = 0 in the subdomain x > x_0, and p is unknown in the subdomain x < x_0. Here, x_0 is unknown and should be found from the second boundary condition. The boundary conditions are:

  • at x = 0, \nabla u(0) = 0;
  • at x = x_0, \nabla u(x_0 - \epsilon) = \nabla u(x_0 + \epsilon), and u(x_0) = T_c, where T_c = 10 is a constant, \epsilon is the smallest positive;
  • at x = L, u(L) = T_b.

So far I came up with this code, where I simply define x_0 and p as some constants and apply Dirichlet BC at x_0 for both equations. But I really don’t understand how to solve it in the case when x_0 and p are unknown.


from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
L = 10
nx = 750        #The number of cells
xL = 0
x0 = 2         #The minimum point (inclusive)
xR = L         #The maximum point (inclusive)
mesh1 = IntervalMesh(nx, xL, x0)
V1 = FunctionSpace(mesh1, 'CG', 2)
mesh2 = IntervalMesh(nx, x0, xR)
V2 = FunctionSpace(mesh2, 'CG', 2)
Tc = 10.0
Tb = 2.0
def DirichletBoundary_xR(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], Constant(xR), tol)
def DirichletBoundary_x0(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], Constant(x0), tol)
bc_xR = DirichletBC(V2, Constant(Tb), DirichletBoundary_xR)
bc_x0 = DirichletBC(V2, Constant(Tc), DirichletBoundary_x0)
bc_x01 = DirichletBC(V1, Constant(Tc), DirichletBoundary_x0)
bcs2 = [bc_x0, bc_xR]
bc1 = [bc_x01]
u1 = TrialFunction(V1)
u2 = TrialFunction(V2)
v1 = TestFunction(V1)
v2 = TestFunction(V2)
p = Tc + 7.0
f1 = Constant(-Tb - p)
f2 = Constant(-Tb)
A = Constant(2e3)
a1 = - 1/A * dot(grad(u1), grad(v1)) * dx - dot(u1, v1) * dx
L1 = f1 * v1 * dx
a2 = - 1/A * dot(grad(u2), grad(v2)) * dx - dot(u2, v2) * dx
L2 = f2 * v2 * dx
u1 = Function(V1)
solve(a1 == L1, u1, bc1)
u2 = Function(V2)
solve(a2 == L2, u2, bcs2)
plot(u1, label='u1')
plot(u2, label='u2')
plt.legend()
plt.show()

Thanks.