如何解决将 FiPy 的 Flow.stokesCavity 示例改编为一维场景
我已经按照建议使用 FiPy 中的 Stokes Cavity 示例攻击了我的 1d 可压缩流(我真的想将 FiPy 用于该模型)。但它仍然不起作用。在扫描期间,我收到错误“对于 0 级解变量,系数必须为 0 级。”。我非常想知道如何解决这个问题,但也很想知道为什么会出现这个错误。
使用的代码:
from fipy import CellVariable,FaceVariable,Grid1D,DiffusionTerm,Viewer
from fipy.tools import numerix
import matplotlib.pyplot as plt
L = 1.0
N = 50
dL = L / N
viscosity = 1
#0.8 for pressure and 0.5 for veLocity are typical relaxation values for SIMPLE
pressureRelaxation = 0.8
veLocityRelaxation = 0.5
if __name__ == '__main__':
sweeps = 300
else:
sweeps = 5
mesh = Grid1D(nx=N,dx=dL)
pressure = CellVariable(mesh=mesh,name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVeLocity = CellVariable(mesh=mesh,name='X veLocity')
veLocity = FaceVariable(mesh=mesh,rank = 1)
xVeLocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([1.])
ap = CellVariable(mesh=mesh,value=1.)
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._celldistances
pressureCorrectionEq = DiffusionTerm(coeff=coeff) - veLocity.divergence
from fipy.variables.faceGradVariable import _FaceGradVariable
volume = CellVariable(mesh=mesh,value=mesh.cellVolumes,name='Volume')
contrvolume=volume.arithmeticFaceValue
xVeLocity.constrain(10.,mesh.facesRight | mesh.facesLeft )
X = mesh.faceCenters
pressureCorrection.constrain(0.,mesh.facesLeft)
from builtins import range
for sweep in range(sweeps):
## solve the Stokes equations to get starred values
xVeLocityEq.cacheMatrix()
xres = xVeLocityEq.sweep(var=xVeLocity,underRelaxation=veLocityRelaxation)
xmat = xVeLocityEq.matrix
## update the ap coefficient from the matrix diagonal
ap[:] = -numerix.asarray(xmat.takeDiagonal())
## update the face veLocities based on starred values with the
## Rhie-Chow correction.
## cell pressure gradient
presgrad = pressure.grad
## face pressure gradient
facepresgrad = _FaceGradVariable(pressure)
veLocity[0] = xVeLocity.arithmeticFaceValue \
+ contrvolume / ap.arithmeticFaceValue * \
(presgrad[0].arithmeticFaceValue-facepresgrad[0])
veLocity[...,mesh.exteriorFaces.value] = 0.
## solve the pressure correction equation
pressureCorrectionEq.cacheRHsvector()
## left bottom point must remain at pressure 0,so no correction
pres = pressureCorrectionEq.sweep(var=pressureCorrection)
rhs = pressureCorrectionEq.RHsvector
## update the pressure using the corrected value
pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
## update the veLocity using the corrected pressure
xVeLocity.setValue(xVeLocity - pressureCorrection.grad[0] / \
ap * mesh.cellVolumes)
if __name__ == '__main__':
if sweep%10 == 0:
print('sweep:',sweep,',x residual:',xres,\
',p residual:',pres,continuity:',max(abs(rhs)))
viewer.plot()
解决方法
这是一个微妙的事情,但是 fipy 没有将单位向量识别为向量。试试:
xVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([[1.]])
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。