This notebook will setup and solve the steady state heat equation:
$$ \nabla(k\nabla)T = h $$where $k$ is the diffusivity, T the temperature field and $h$ the source term.
Keywords: initial conditions, boundary conditions, heat equation
import underworld as uw
import glucifer
# Set box size.
boxHeight = 1.0
boxLength = 2.0
# Set the resolution.
resx = 16
resy = 8
mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"),
elementRes = (resx, resy),
minCoord = (0., 0.),
maxCoord = (boxLength, boxHeight))
# gLucifer visualisation of mesh
fig = glucifer.Figure( figsize=(800,400) )
fig.append(glucifer.objects.Mesh( mesh ))
fig.show()
# Create mesh variables for the temperature field & initialise.
temperatureField = mesh.add_variable( nodeDofCount=1 )
temperatureField.data[:] = 0.
We first wish to determine which vertices will be flagged as boundary conditions. To do this, we create set objects which contain the vertex indices for which the condition will apply. The specialSets dictionary (on the mesh object) contains the sets we usually require. Note however that you may construct sets which contain indices for any vertex on the mesh.
mesh.specialSets.keys()
The vertices along the bottom wall is given by 'MinJ_VertexSet', the top wall given by the 'MaxJ_VertexSet'.
Construct sets for the horizontal walls:
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
Create Dirichlet boundary conditions for the jWalls
and provide values.
tempBC = uw.conditions.DirichletCondition( variable=temperatureField, indexSetsPerDof=(jWalls,) )
# set bottom wall temperature bc
for index in mesh.specialSets["MinJ_VertexSet"]:
temperatureField.data[index] = 1.0
# set top wall temperature bc
for index in mesh.specialSets["MaxJ_VertexSet"]:
temperatureField.data[index] = 0.0
# gLucifer visualisation of temperature field & mesh
fig.append( glucifer.objects.Surface( mesh, temperatureField, colours="blue white red" ) )
fig.show()
Temperature field, diffusivity and boundary conditions are passed to the SteadyStateHeat system function.
heatequation = uw.systems.SteadyStateHeat(temperatureField = temperatureField,
fn_diffusivity = 1.0,
conditions = tempBC)
Solve the heat equation.
# get the default heat equation solver
heatsolver = uw.systems.Solver(heatequation)
# solve
heatsolver.solve()
# gLucifer visualisation of temperature field & mesh
fig.show()