CODEWITHSUNDEEP
pythonfipy
DavidT
DavidT

Reputation: 79

conduction-diffusion heat 2D model with non constant capacity

I've spent quite some time on developing a two-dimensional heat conduction-diffusion model for steady state approximation.

PDE

For simplicity consider a laminar shear film, i.e. zero velocity at the bottom and a constant velocity increase.

shear

Heat capacity might either be constant or linear increasing over temperature.

Boundary conditions are a constant inlet temperature (left) and constant input flux (top) while all exterior faces are forced to have no gradient.

See the code here.

When using a constant heat capacity input power is equal to output power.

input = 50.00e3 W
ouput = 50.00e3 W

When using a non-constant heat capacity they differ significantly. The more heat capacity varies over temperature, the more input and output differ.

input = 50.00e3 W
ouput = 33.78e3 W

Introducing a variable velocity coefficient (here v * c * rho) was done as stated in the fipy FAQ (which only explicitly shows an example for diffusion terms). Grid resolution does not change output power. So I'd say it is not a grid-problem. I've also tried adding a transient term and solving for a very high time step, which does not change the solution.

I fear I've done something terribly wrong when defining the convection term, but cannot find the error. Also I'm confused if fipy is able to mix theta (a rank=0 cell variable) with the velocity (a rank=1 cell variable) and then cast them into a face variable, which is necessary for the convection term.

Upvotes: 1

Views: 1218

Answers (1)

jeguyer
jeguyer

Reputation: 2484

Per the divergence theorem

divergence theorem

I would calculate the surface flux with

Cp = mymodel.fluid.capacity(solution, use_constant_cp)
veloc = fipy.CellVariable(mesh=mesh, value=0., rank=1, name='velocity')
veloc[0] = mymodel.shear * mesh.y
R = ((Cp * solution * veloc).faceValue.dot(mesh._orientedAreaProjections) * mesh.facesRight * mymodel.fluid.rho).sum()
L = ((Cp * solution * veloc).faceValue.dot(mesh._orientedAreaProjections) * mesh.facesLeft * mymodel.fluid.rho).sum()
print "{:.3e} J/s received.".format((R+L).value)

I get 4.958e+04 J/s received. The answer improves with x resolution.

Note that because this uses the velocity vector, L is the flux in and R is the flux out and so we add them to get the difference. The _orientedAreaProjections vector points into the domain for all exterior faces, so the dot product is positive when the flux enters the domain and negative when it exits. Because we're integrating over the entire exterior boundary, you could just write

J_dot_n = ((Cp * solution * veloc).faceValue.dot(mesh._orientedAreaProjections) * (mesh.facesLeft + mesh.facesRight) * mymodel.fluid.rho).sum()
print "{:.3e} J/s received.".format(J_dot_n.value)

Similarly, I would calculate the input heat flux with (mymodel.flux * mesh._faceAreas * mesh.facesTop).sum().

What you've calculated, I think, is

If you want to calculate the volume integral form of the divergence theorem, you can do that, but it would be

velocF = fipy.FaceVariable(mesh=mesh, value=0., rank=1, name='velocity')
velocF[0] = mymodel.shear * mesh.faceCenters[1]
((Cp * solution).faceValue * velocF * mymodel.fluid.rho).divergence.cellVolumeAverage * mesh.cellVolumes.sum()

_faceAreas and _orientedAreaProjections come up often enough that we should make them part of the public API.

[Edited for clarity to address questions that came up in the comments]

Upvotes: 1

Related Questions