0
$\begingroup$

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?

$\endgroup$
  • 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 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$ – Wolfgang Bangerth Jan 9 at 18:23
  • $\begingroup$ @balborian I have added the BC details $\endgroup$ – Black Heart Jan 9 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$ – Alex Trounev Jan 10 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.