0

我想用一些源项来解决流体流动问题。但是我无法按照我的要求在最后用源项写出方程。有没有办法将 FaceVariable 转换为 cellVariable。

我收到此错误

TypeError：系数不能是 FaceVariable。

代码如下

from fipy import CellVariable, FaceVariable, Grid2D, DiffusionTerm, Viewer, ExponentialConvectionTerm
from fipy.tools import numerix
from fipy.meshes.gmshMesh import Gmsh2D
import numpy as np

L = 1.0
N = 100
dL = L / N
viscosity = 1
U = 1.

meshsize = 1./N
geo = '''
cl__1 = 1;
Point(1) = {1, 1, 0, %f};
Point(2) = {2, 1, 0, %f};
Point(3) = {2, 2, 0, %f};
Point(4) = {1, 2, 0, %f};
Line(1) = {4, 3};
Line(2) = {3, 2};
Line(3) = {2, 1};
Line(4) = {1, 4};
Line Loop(6) = {1, 2, 3, 4};
Plane Surface(6) = {6};
Transfinite Line {4, 1, 2, 3} = %d Using Progression 1;
Transfinite Surface {6};
Recombine Surface {6};
'''

geo = geo %(meshsize ,meshsize, meshsize, meshsize , 1/meshsize )

mesh = Gmsh2D(geo , background=None)

#0.8 for pressure and 0.5 for velocity are typical relaxation values for SIMPLE
pressureRelaxation = 0.1
velocityRelaxation = 0.1
if __name__ == '__main__':
    sweeps = 100
else:
    sweeps = 5

x = mesh.cellCenters.value[0]
y = mesh.cellCenters.value[1]

sourcevx = CellVariable(mesh=mesh , name = 'sourcevx')
sourcevy = CellVariable(mesh=mesh , name = 'sourcevy')
source_p = FaceVariable(mesh=mesh , rank=1)
source_pcorr = CellVariable(mesh=mesh , name = 'source_pcor')

pressure = CellVariable(mesh=mesh, name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVelocity = CellVariable(mesh=mesh, name='X velocity')
yVelocity = CellVariable(mesh=mesh, name='Y velocity')

velocity = FaceVariable(mesh=mesh, rank=1)

sourcevx.setValue(-x**2*numerix.sin(x*y) - 2*x*numerix.cos(x**2 + y**2) - \
            y**2*numerix.sin(x*y))
sourcevy.setValue(-x*y*numerix.sin(x*y) - 2*y*numerix.cos(x**2 + y**2) \
           + numerix.cos(x*y) - y**3*numerix.sin(x*y)/x - \
           3*y**2*numerix.cos(x*y)/x**2 + 6*y*numerix.sin(x*y)/x**3 \
           + 6*numerix.cos(x*y)/x**4)
xVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([1.,0.]) \
               == sourcevx
yVelocityEq = DiffusionTerm(coeff=viscosity) - pressure.grad.dot([0.,1.]) \
               == sourcevy



ap = CellVariable(mesh=mesh, value=1.)
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances

x , y = mesh.faceCenters

source_p.setValue(coeff * (-4*x**2*numerix.sin(x**2 + y**2) \
            - 4*y**2*numerix.sin(x**2 + y**2) + 4*numerix.cos(x**2 + y**2))\
                  - 2*y*numerix.cos(x*y) )

pressureCorrectionEq = DiffusionTerm(coeff=coeff) - velocity.divergence \
                       == source_p
4

1 回答 1

0

在 FiPy 中，没有简单的方法可以从面到单元中心进行插值，因为它很少需要。就您而言，我认为也不需要。source_p面变量可以表述为单元变量。

source_p = CellVariable(mesh=mesh, rank=1)
ap = CellVariable(mesh=mesh, value=1.)
coeff = 1 / ap * mesh.cellVolumes
x ,y = mesh.cellCenters

source_p.setValue(coeff * (-4*x**2*numerix.sin(x**2 + y**2) \
            - 4*y**2*numerix.sin(x**2 + y**2) + 4*numerix.cos(x**2 + y**2))\
                  - 2*y*numerix.cos(x*y) )

pressureCorrectionEq = DiffusionTerm(coeff=coeff) - velocity.divergence \
                       == source_p

如果我们定义coeff为开始的单元格变量，那么我们不需要插回单元格。coeff这确实会稍微改变离散化。如果您对此感到担心，那么您可以定义两个版本coeff

coeff_cell = 1 / ap * mesh.cellVolumes
coeff_face = 1 / ap.arithmeticFaceValue * mesh._faceAreas * mesh._cellDistances

其中coeff_face用作扩散系数并coeff_cell用于源项的表达式中。

于 2018-08-17T15:03:03.730 回答