Having trouble solving Population balance PDE in FiPy

Question

I am a beginner in using FiPy and also in working with Finite Volume Methods in general. So please forgive me if my questions are naïve

I am trying to solve a population balance equation of the form PDE with the boundary conditions as follows Boundary Conditions

The value of x ranges from 0.5e-6 to 1000e-6 and I want to use a uniform mesh.

This is the code that I currently have and it is not working. In the code I have assumed that the terms G and b_0 are time invariant to make the problem easier. Can someone please guide me about how I should define the terms in the above problem to get a working solution?

from fipy import *
import scipy

b0 = 11.78*10**7
g0 = 2.549e-7

#Define Mesh
nx = 1000
dx = 10e-7
mesh = Grid1D(nx = nx, dx = dx)
x = mesh.cellCenters[0]

#Define coefficients
G = FaceVariable(mesh=mesh, value = g0)

#Defining the solution Variable
phi = CellVariable(name = 'Solution Variable', mesh = mesh, 
                    hasOld = True)

#Set initial condition
g = (scipy.stats.norm(loc = x[500], scale = 5e-6).pdf(x))*b0
phi.setValue(g)


#Set Boundary Condition
(phi*G).constrain(b0, where = mesh.facesLeft)
(phi*G).constrain(0, where = mesh.facesRight)


#Define the convective coefficient
conv_coefficient = G

#Define the Equation
eqn = TransientTerm() == -ConvectionTerm(coeff = conv_coefficient)


steps = 1
for step in range(steps):
    dt = 10
    eqn.solve(var = phi, dt = dt)


plt.plot(phi.value)

Answer 1

Your script has several problems:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [1], in <cell line: 21>()
     17 phi = CellVariable(name = 'Solution Variable', mesh = mesh, 
     18                     hasOld = True)
     20 #Set initial condition
---> 21 g = (scipy.stats.norm(loc = x[500], scale = 5e-6).pdf(x))*b0
     22 phi.setValue(g)
     25 #Set Boundary Condition

AttributeError: module 'scipy' has no attribute 'stats'

You must import specific scipy modules before using them, e.g., import scipy.stats.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [2], in <cell line: 26>()
     22 phi.setValue(g)
     25 #Set Boundary Condition
---> 26 (phi*G).constrain(b0, where = mesh.facesLeft)
     27 (phi*G).constrain(0, where = mesh.facesRight)
     30 #Define the convective coefficient

TypeError: unsupported operand type(s) for *: 'CellVariable' and 'FaceVariable'

phi is defined at cell centers. G is defined at the faces between cells. FiPy doesn't know how to multiply variables in two different locations.

As written, your first boundary condition can be rearranged to the Dirichlet condition , or

phi.constrain(b0/G, where=mesh.facesLeft)

Your second boundary condition can be rearranged to the Dirichlet condition , or

phi.constrain(0, where=mesh.facesRight)

---------------------------------------------------------------------------
VectorCoeffError                          Traceback (most recent call last)
Input In [3], in <cell line: 34>()
     31 conv_coefficient = G
     33 #Define the Equation
---> 34 eqn = TransientTerm() == -ConvectionTerm(coeff = conv_coefficient)
     37 steps = 1
     38 for step in range(steps):

VectorCoeffError: The coefficient must be a vector value.

The coeff argument to ConvectionTerm must be a vector (rank-1), but conv_coeff = G is rank-0. You can make G rank-1 (and probably should), but then your boundary condition doesn't make sense, as expresses "vector * scalar = scalar". Since and are only functions of time, you can make this work by changing the definitions to

G = Variable(value = g0)

and

conv_coefficient = G * [[1]]

The full code is

from fipy import *
import scipy.stats

b0 = 11.78*10**7
g0 = 2.549e-7

#Define Mesh
nx = 1000
dx = 10e-7
mesh = Grid1D(nx = nx, dx = dx)
x = mesh.cellCenters[0]

#Define coefficients
G = Variable(value = g0)

#Defining the solution Variable
phi = CellVariable(name = 'Solution Variable', mesh = mesh, 
                    hasOld = True)

#Set initial condition
g = (scipy.stats.norm(loc = x[500], scale = 5e-6).pdf(x))*b0
phi.setValue(g)


#Set Boundary Condition
phi.constrain(b0/G, where = mesh.facesLeft)
phi.constrain(0, where = mesh.facesRight)


#Define the convective coefficient
conv_coefficient = G * [[1]]

#Define the Equation
eqn = TransientTerm() == -ConvectionTerm(coeff = conv_coefficient)


steps = 1
for step in range(steps):
    dt = 10
    eqn.solve(var = phi, dt = dt)


Viewer(vars=phi).plot()

I still have doubts as your expressions aren't consistent about scalar and vector quantities. I strongly recommend against trying to solve PDEs expressed in 1D. It's almost always misleading. You need to know what is scalar and what is vector and this is hard to keep track of in 1D.