I'm trying to simulate the heat diffusion in a 3D piston. I marked the boundaries on GMSH.

I have used a Dirichlet BC of 300 on the top face of piston. But the results look abnormal. There is a change in temperature along the edge on the top face, which should not be case.

Here is code snippet

from __future__ import print_function
from fenics import *
import numpy as np

mesh = Mesh("piston-3d.xml");

V = FunctionSpace(mesh, 'P', 1)

subdomains = MeshFunction("size_t", mesh, "piston-3d_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "piston-3d_facet_region.xml")

#Define neumann bcs
g_1= Constant(300.0)
g_2 = Constant(100.0)
g_3 = Constant(30.0)
g_4 = Constant(0.0) // On the cross-sectional faces, Neumann
g_5 = Constant(0.0) // On the cross-sectional faces, Neumann

bc1 = DirichletBC(V, g_1, boundaries, 1) // Top face of piston
bc2 = DirichletBC(V, g_2, boundaries, 2) // Liner face of the piston
bc3 = DirichletBC(V, g_3, boundaries, 3) // Inner side of piston and bottom

bcs = [bc1, bc2, bc3]

#define measures
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

u_D = Constant(100.0)

# Define initial value
u_n = interpolate(u_D, V)

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

a = dot(grad(u), grad(v))*dx 

L = f*v*dx - (g_4*v*ds(4) + g_5*v*ds(5))

# Time-stepping
u = Function(V)

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

vtkfile = File("solution_steady.pvd")   
vtkfile << u

The result looks like this

My question is why is the temperature of Top face not remaining constant (300)? Is it because of the boundary sharing?

  • 2
    $\begingroup$ It is hard to tell without knowing the mesh marked boundaries. Also, it is more likely to get an answer at the FEniCS discourse group fenicsproject.discourse.group $\endgroup$
    – balborian
    Jan 9, 2020 at 18:13
  • 1
    $\begingroup$ But concretely, you need to decide which boundary conditions should hold on the nodes at the edge of the quarter circle. It can't be both the top and side boundary conditions. $\endgroup$ Jan 9, 2020 at 18:23
  • $\begingroup$ @balborian I have added the BC details $\endgroup$ Jan 9, 2020 at 18:30
  • 3
    $\begingroup$ The result looks as expected, since the surface temperature is set to 100. The transition zone from 300 to 100 is visible in the figure. $\endgroup$ Jan 10, 2020 at 17:08


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.