usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
504 stars 148 forks source link

How can I solve this error 'The coefficient must be a vector value' in the following code? #938

Closed ghost closed 1 year ago

ghost commented 1 year ago

`def Flow():

Import fipy

from fipy import CellVariable, FaceVariable, Grid2D, DiffusionTerm, Viewer, TransientTerm, ImplicitSourceTerm, ConvectionTerm
############################################################################################################################

# Define the initial & physical parameters
viscosity = 0.1
diffusivity = 1.0
interfacialAreaCoeff = 0.1
u = 1.0
v = 0.0
k = 1.0
############################################################################################################################

# Set up the domain
L1 = 5.0    #cm
L2 = 20    #cm
x = 10.0
y = 2.0
dy = L1 / y
dx = L2 / x
mesh = Grid2D(nx=x, ny=y, dx=dx, dy=dy)
x, y = mesh.cellCenters
############################################################################################################################

# Define the fluid velocity
velocity = FaceVariable(name="Velocity", mesh=mesh, value=(u,v))

#Define the pressure
pressure = CellVariable(name="Pressure", mesh=mesh)

# Define the species concentration
C = CellVariable(name="Concentration", mesh=mesh, hasOld=True)

# Define the volume fraction
phi = CellVariable(name="Volume Fraction", mesh=mesh, hasOld=True)

# Define the interfacial area
A = CellVariable(name="Interfacial Area", mesh=mesh, hasOld=True)
############################################################################################################################

# Define the initial and boundary conditions
C.setValue(0.0)
C.constrain(1.0, mesh.facesBottom)
phi.setValue(0.5)
phi.constrain(0.5, mesh.facesRight)
A.setValue(0.0)
velocity.constrain(0.0, mesh.exteriorFaces)
velocity.constrain(1.0, mesh.facesLeft)
############################################################################################################################

# Define the chemical reaction
reaction = -k * C

# Define the transport coefficients
fluidCoeff = FaceVariable(mesh=mesh, value=viscosity)
speciesCoeff = FaceVariable(mesh=mesh, value=diffusivity)
interfacialAreaCoeff = FaceVariable(mesh=mesh, value=interfacialAreaCoeff)
############################################################################################################################

# Define the Navier-Stokes equations
equ1 = TransientTerm(var=velocity) == - ConvectionTerm(coeff=fluidCoeff, var=velocity) - pressure.diverdence
equ2 = TransientTerm(var=C) == DiffusionTerm(coeff=speciesCoeff, var=C) + reaction
equ3 = ImplicitSourceTerm(coeff=1./fluidCoeff.divergence, var=pressure) == (velocity * mesh.faceNormals).divergence
equ4 = TransientTerm(var=phi) == -velocity.divergence
equ5 = TransientTerm(var=A) == - ConvectionTerm(coeff=interfacialAreaCoeff, var=phi) - A / phi * (velocity * mesh.faceNormals).divergence
############################################################################################################################

# Solve the equations
eqns = equ1 & equ2 & equ3 & equ4 & equ5
timeStepDuration = 0.1
steps = 10
for step in range(steps):
    eqns.solve(dt=timeStepDuration)
############################################################################################################################

# Visualize the results
viewer = Viewer(vars=(C, velocity), datamin=0.0, datamax=1.0)
viewer.plot()

Flow() `