0

I want to solve a fluid flow problem with some source terms. But I am not able to write the equation at the end with the source term as desired by me. Is there a way to convert the FaceVariable to cellVariable.

I am getting this error

TypeError: The coefficient can not be a FaceVariable.

the code is as follows

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
jww
  • 97,681
  • 90
  • 411
  • 885
psnbaba
  • 91
  • 3

1 Answers1

0

There is no simple method to interpolate from face to cell centers in FiPy because it is so rarely required. In your case, I don't think it is required either. The source_p face variable can be formulated as a cell variable.

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

If we define coeff to be a cell variable to start with then we don't need to interpolate back to the cells. This does change the discretization for coeff slightly. If you're worried about that then you can define two versions of coeff

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

where the coeff_face is used as the diffusion coefficient and coeff_cell is used in the expression for the source term.

wd15
  • 1,068
  • 5
  • 8