2D-модель теплопроводно-диффузионного тепла с непостоянной емкостью

Question

Я потратил довольно много времени на разработку двумерной модели теплопроводности-диффузии для стационарного приближения.

http://mathurl.com/p3rfqer.png

Для простоты рассмотрим ламинарную сдвиговую пленку, то есть нулевую скорость на дне и постоянную скорость увеличения.

http://mathurl.com/p6e3hau.png

Теплоемкость может быть либо постоянной, либо линейной с увеличением температуры.

Граничными условиями являются постоянная температура на входе (слева) и постоянный входной поток (вверху), в то время как все внешние грани не имеют градиента.

Смотрите код здесь.

При использовании постоянной теплоемкости входная мощность равна выходной мощности.

input = 50.00e3 W
ouput = 50.00e3 W

При использовании непостоянной теплоемкости они существенно различаются. Чем больше теплоемкость изменяется в зависимости от температуры, тем больше расход и расход.

input = 50.00e3 W
ouput = 33.78e3 W

Введение коэффициента переменной скорости (здесь v * c * rho) было сделано, как указано в fipy FAQ (который только явно показывает пример для диффузионных членов). Разрешение сетки не меняет выходной мощности. Так что я бы сказал, что это не проблема сетки. Я также попытался добавить переходный термин и решение для очень высокого временного шага, который не меняет решение.

Боюсь, что я сделал что-то ужасно неправильное при определении термина конвекции, но не могу найти ошибку. Также я запутался, если фипи умеет смешивать theta (переменная ячейки ранга =0) с velocity (переменная ячейки ранга =1), а затем приведите их в переменную лица, что необходимо для термина конвекции.

Answer 1

По теореме о расходимости

http://mathurl.com/oz3wd3c.png

Я бы рассчитал поверхностный поток с

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)

я получил 4.958e+04 J/s received. Ответ улучшается с разрешением х.

Обратите внимание, что, поскольку это использует вектор скорости, L поток в и R это поток, и поэтому мы добавляем их, чтобы получить разницу. _orientedAreaProjections вектор указывает на домен для всех внешних граней, поэтому скалярное произведение положительно, когда поток входит в домен, и отрицательно, когда оно выходит. Поскольку мы интегрируем по всей внешней границе, вы можете просто написать

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)

Точно так же я бы рассчитал входной тепловой поток с (mymodel.flux * mesh._faceAreas * mesh.facesTop).sum(),

То, что вы рассчитали, я думаю,

http://mathurl.com/qzn7wbk.png

Если вы хотите вычислить объемную интегральную форму теоремы о расхождении, вы можете сделать это, но это будет

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 а также _orientedAreaProjections достаточно часто, чтобы мы могли сделать их частью публичного API.

[Отредактировано для ясности, чтобы ответить на вопросы, которые возникли в комментариях]