The Dirichlet set value is inconsistent with the solution value

General
1

Hello everyone, I am doing a heat transfer analysis, I apply a first class boundary condition at the bottom of 298K, but the solution is not the temperature value, what should I do?

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

L = 0.01    
W = 0.01  
nx = 50    
ny = 50    
mesh = RectangleMesh(Point(0, 0), Point(L, W), nx, ny)
V = FunctionSpace(mesh, 'P', 1)

def k(T):   
    return 174.9274 - 0.1067 * T + 5.0067e-5 * T**2 - 7.8349e-9 * T**3
mu = Constant(19300)    
c = Constant(137)     
q = Expression('1e7', degree=0)    
h = Expression('70000', degree=0)    
T_ini = Expression('298', degree=0)    
f = Constant(0.0) 
t_total = 10  
num_steps = 500    
dt = t_total / num_steps  


class BoundaryX0(SubDomain):   
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0]-0.00)<DOLFIN_EPS

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0]-0.01)<DOLFIN_EPS
    
class BoundaryY0(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1]-0.00)<DOLFIN_EPS

class BoundaryY1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1]-0.01)<DOLFIN_EPS

bx0 = BoundaryX0()    
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()

bc_bottom = DirichletBC(V, T_ini, by0)

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)    

bx0.mark(boundary_markers, 0)    
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)

boundary_conditions = {0: {'Neumann': Constant(0)}, 1: {'Neumann': Constant(0)}, 2: {'Dirichlet': T_ini}, 3: {'Neumann': q}}

bcs = []    
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
        
T_n = interpolate(T_ini, V)

T = Function(V)
v = TestFunction(V)

integrals_N = []    
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        #if boundary_conditions[i]['Neumann'] != 0:
        q = boundary_conditions[i]['Neumann']
        integrals_N.append(q*v*ds(i))

integrals_R = []   
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        h, T_envi = boundary_conditions[i]['Robin']
        integrals_R.append(h*(T - T_envi)*v*ds(i))

F = mu*c*T*v*dx + dt*k(T)*dot(grad(T), grad(v))*dx - mu*c*T_n*v*dx - dt*f*v*dx - dt*sum(integrals_N)

progress = Progress('Time-stepping', num_steps)

t = 0
plt.figure(figsize=(12,8))
for n in range(num_steps):
    t += dt
    solve(F == 0, T, bcs)
    p = plot(T, cmap='coolwarm')
    T_n.assign(T)
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)
plt.colorbar(p)
plt.show()
3

This is my result, the temperature set at the bottom is 298K, but the solution is 280K.
1686832999669

4

This is likely a visualisation issue with matplotlib applying some default setting of ticks on the deca-values. What do you see if you evaluate the finite element function along the bottom boundary?

5

Thank you very much for your answer, I added the following code before and after solving:

t = 0
temperatures=[]
plt.figure(figsize=(12,8))
for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(F == 0, T, bcs)
    p = plot(T, cmap='coolwarm')
    
    point = Point(0.005, 0)
    temperature = T(point)
    
    temperatures.append(temperature)
    
    # Update previous solution
    T_n.assign(T)
    
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)

plt.colorbar(p)
plt.show()

The code collects the temperature value at the bottom,Turns out the bottom is 298K.This really surprised me, but how do I visualize it accurately?

6

Just now, I output the temperature change at a certain point at the bottom boundary with the time step, and the result is consistent with my setting. Then I go through the bottom node and find that it is also consistent with the setting, which seems to be a visualization problem.