将 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()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

fipy python