Transient heat equation in a two layered medium: comparison to theory


235
views
0
7 months ago by
I am new in FEniCS. I have implemented the transient heat equation in a two layered medium (with different conductivity and diffusivity) modifying this code. Boundary conditions are, at one side of the box 310.15 K and 315.15 K on the other side. I have also implemented the heat flow at stationary state using the electromagnetic analogy but I have different results. am I doing something wrong with the code?

Thank you very much in advance
from __future__ import print_function
from dolfin import *
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 3.0e-3)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 3.0e-3)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.0, 1.2e-3)) and between(x[0], (0.0, 3.0e-3)))

class K(Expression):
    def __init__(self, subdomains, k_0, k_1, **kwargs):
        self.subdomains = subdomains
        self.k_0 = k_0
        self.k_1 = k_1        
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1

class Alfa(Expression):
    def __init__(self, subdomains, alfa_0, alfa_1, **kwargs):
        self.subdomains = subdomains
        self.alfa_0 = alfa_0
        self.alfa_1 = alfa_1        
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 1:
            values[0] = self.alfa_0
        else:
            values[0] = self.alfa_1

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
# Define mesh
mesh = RectangleMesh(Point(0.0, 0.0),Point(3.0e-3, 3.0e-3), 10, 10)
# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
# Define input paramters
T = 51.0            # final time
num_steps = 201     # number of time steps
dt = T / num_steps # time step size
t = 0.0
k0 = 0.65#Thermal conductivity medium 1
k1 = 0.43#Thermal conductivity medium 2
a0 = Constant(0.65/(1000.0*4180.0))#Thermal diffusivity medium 1
a1 = Constant(0.43/(1178.0*2274.0))#Thermal diffusivity medium 2
g_L = Constant(0.0)
g_R = Constant(0.0)
# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)
u_n = interpolate(Constant(310.15), V)
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 310.15+5.0, boundaries, 2),
       DirichletBC(V, 310.15, boundaries, 4)]
# Define new measures associated with the interior domains and
# exterior boundaries
dx = dx(subdomain_data=domains)
ds = ds(subdomain_data=boundaries)
kappa = K(domains, k0, k1, degree=2)
alfa = Alfa(domains, a0, a1, degree=2)

F = (u*v*dx+inner(dt*alfa*grad(u), grad(v))*dx 
     - g_L*v*ds(1) - g_R*v*ds(3)
     - u_n*v*dx)

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)
# Solve problem
u = Function(V)
flujo_t = []

n = FacetNormal(mesh)
for nn in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(a == L, u, bcs)
    # Update previous solution
    u_n.assign(u)
    # Compute the flux
    flux = -kappa*dot(grad(u),n)*ds(4)
    total_flux = assemble(flux)
    flujo_t.append(total_flux)

#Compute steady numerical flux
flux = -(kappa)*dot(grad(u),n)*ds(4)
total_flux = assemble(flux)
size = 3e-3
print('numerical flux:',total_flux/size)

k1 = 0.43
k2 = 0.65

width1 = 1.2e-3
width2 = 1.8e-3

Resistence = width1/(k1)+width2/(k2) #Thermal resistence
DeltaT = 315.15-310.15 #Thermal difference
print('theoretical flux:',DeltaT/Resistence)

# Plot solution and gradient
plot(u, title="u")
plot(grad(u), title="Projected grad(u)")
plot(domains)
plot(mesh)

interactive()

import numpy as np
import matplotlib.pyplot as plt

flujo = np.asarray(flujo_t)/size
tiempo = np.arange(0.0,len(flujo)*dt,dt)

plt.subplot(111)
plt.plot(tiempo,flujo, '.k', linewidth=2)  # magnify w
plt.plot(tiempo,(DeltaT/Resistence)*np.ones(len(flujo)), 'k', linewidth=2)
plt.grid(True)
plt.xlabel('$time$')
plt.ylabel('$flux$')
plt.legend(['Numerical','Steady theory'], loc='upper left')

plt.show()
​
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »