Temperature jump near boundry condition

1

Dear colleagues,
I am using the heat equation based upon the FENICS heat_eq tutorial as per following code:

from __future__ import print_function
from fenics import *
import time

T = 1.0            # final time
num_steps = 1000     # number of time steps
dt = T / num_steps # time step size

# Import mesh and define function space
mesh = UnitCubeMesh(12,12,12);
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
def left_boundary(x, on_boundary):
    if (x[1]<DOLFIN_EPS):
        return on_boundary

def right_boundary(x, on_boundary):
    if (abs(x[1] - 10) < DOLFIN_EPS):
        return on_boundary

def bottom_boundary(x, on_boundary):
    if (x[1] < DOLFIN_EPS):
        return on_boundary

        
bcs = []
#bc_left = DirichletBC(V, Constant(0), left_boundary)
#bcs.append(bc_left)
#bc_right = DirichletBC(V, Constant(1000), right_boundary)
#bcs.append(bc_right)
bc_bottom = DirichletBC(V, Constant(400), bottom_boundary)
bcs.append(bc_bottom)

#bc = [bc_left,bc_right]


# Define initial value
#u_0 = Expression('0*1',
#                 degree=2, a=5)

#u_n = interpolate(u_0, V)
u_n = interpolate(Constant(800), V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Create VTK file for saving solution
vtkfile = File('heat_equation/solution.pvd')

# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bcs)

    # Save to file and plot solution
    #if (n%20 == 0):
        #print('Printing the result at time step number =', n)
    vtkfile << (u, t)
        #plot(u)

    # Update previous solution
    u_n.assign(u)
    print('Time step number  =', n)
    print('Time  =', t)

The boundary condition is Dirichlet BC of fixed temperature at 400 at bottom. The initial solution declared is 800. When I start the analysis, in early time steps there I observe an unusual jump near BC where temperature crosses the value of 800 as shown in attached graph and contour plot.
I tried refining the mesh and it vanished but for smaller time steps, it appears again. Did anyone faced the same issue in FEM analysis? Any reasons of this behaviour and suggestions to resolve it? Thank you!

jump_graph
jump_graph1437Ă—807 9.73 KB

temperature_jump_neastrong textr_bottom
temperature_jump_neastrong textr_bottom1344Ă—768 148 KB

2

Having a smooth transition from your Dirichlet boundary condition will probably yield a physical result.

2 Likes
3

Thanks @nate for your informative response. Other than reducing value of Dirichlet BC, in terms of FEM can one achieve this smooth transition using higher order element or any other traditional methodology?

4

Changing the value imposed in the DirichletBC isn’t what I had in mind. Rather that your initial state:

u_n = interpolate(Constant(800), V)

would benefit from being a smooth transition from your DirichletBC value of 400, to your maximum value of 800. You could pick any smooth transition function that you like. This transition can even occur over a very small spatial region.

5

Dear @nate, Thanks a lot for your suggestion. I tried with the smoothing function of expr = Expression(' 400 + (0.5 + 0.5 * tanh( (180/pi) * ((x[2]-a)/b)) ) * 400 ', a=0.8333, b=3, degree = FE_degree) near the Dirchlet BC to smooth out the initial solution. In my problem, I also involve the Robin BC on the boundaries of the cube. So during initial time steps I also get kind of a jump near the Robin BC faces. As now the simulation is already on the run and I cannot predict what would be the solution because of Robin BC near the edges unlike case of Dirchlet BC therefore I also cannot define smoothing function earlier.
In your opinion what would be the suggestion to apply it when the simulation is running i.e. applying smoothing near the boundary on the go when simulation is up and running? My mwe is following:

Dirchlet BC
"""

from __future__ import print_function
from fenics import *
import time

T = 1.0            # final time
num_steps = 1000     # number of time steps
dt = T / num_steps # time step size

# Import mesh and define function space
mesh = UnitCubeMesh(12,12,12);
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
def left_boundary(x, on_boundary):
    if (x[1]<DOLFIN_EPS):
        return on_boundary

def right_boundary(x, on_boundary):
    if (abs(x[1] - 1) < DOLFIN_EPS):
        return on_boundary

def bottom_boundary(x, on_boundary):
    if (x[1] < DOLFIN_EPS):
        return on_boundary

#def top_boundary(x, on_boundary):
#    if (abs(x[1] - 1) < DOLFIN_EPS):
#        return on_boundary
class top_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1-1e-14, 1)

        
bcs = []
#bc_left = DirichletBC(V, Constant(0), left_boundary)
#bcs.append(bc_left)
#bc_right = DirichletBC(V, Constant(1000), right_boundary)
#bcs.append(bc_right)
bc_bottom = DirichletBC(V, Constant(800), bottom_boundary)
bcs.append(bc_bottom)

#bc = [bc_left,bc_right]


# Define initial value
#u_0 = Expression('0*1',
#                 degree=2, a=5)

#u_n = interpolate(u_0, V)
u_n = interpolate(Constant(800), V)

## Define convection BC

# parameters
conv_coeff = 20
thermal_cond = 1
ambient_temp = 400

r = Constant(conv_coeff/thermal_cond)
s = Constant(ambient_temp)

# Mark boundaries for convection
boundary_marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 9999)
boundary_class = top_boundary(); boundary_class.mark(boundary_marker, 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

# surface integral for convection    
ds = Measure('ds', domain = mesh, subdomain_data = boundary_marker)
local_factor = dt*thermal_cond*r*u*v*ds(1)
ambient_factor = dt*thermal_cond*r*s*v*ds(1)


F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx + local_factor - ambient_factor
a, L = lhs(F), rhs(F)

# Create VTK file for saving solution
vtkfile = File('heat_equation/solution.pvd')

# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bcs)

    # Save to file and plot solution
    #if (n%20 == 0):
        #print('Printing the result at time step number =', n)
    vtkfile << (u, t)
        #plot(u)

    # Update previous solution
    u_n.assign(u)
    print('Time step number  =', n)
    print('Time  =', t)
    min=u.vector().get_local().min()
    max=u.vector().get_local().max()
    print('minima as well as the maxima in domain', min, max, n, t)

Screenshot from 2021-06-09 19-46-19
Screenshot from 2021-06-09 19-46-191453Ă—820 189 KB