在 FiPy 中求解总体平衡 PDE 时遇到问题

Question

我是使用FIPY以及通常使用有限体积方法的初学者。因此，如果我的问题很幼稚，请原谅我

，我试图解决形式的人口平衡方程 pde 边界条件如下 边界条件

x范围为0.5e-6到1000e-6的值，我想使用均匀的网眼。

这是我当前拥有的代码，并且无法正常工作。在代码中，我认为术语g和b_0是时间不变，可以使问题更容易。有人可以指导我如何定义上述问题中的术语以获取工作解决方案吗？

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

您的脚本有几个问题：

---------------------------------------------------------------------------
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'

您必须在使用它们之前导入特定的Scipy模块，例如导入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在单元中心定义。 g在单元格之间的面上定义。 Fipy不知道如何在两个不同的位置繁殖变量。

如书面，您的第一个边界条件可以重新排列到dirichlet条件

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

​=“ https://i.sstatic.net/uihij.png” rel =“ nofollow noreferrer”> ，或

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.

coeff对对流必须是向量（rank-1），但是cons_coeff = g是等级 - 0。您可以制作g rank-1（可能应该），但是您的边界条件没有意义，因为 表示“ vector * scalar = scalar = scalar”。因为 和 只是时间的函数，您可以通过将定义更改为

G = Variable(value = g0)

，而

conv_coefficient = G * [[1]]

完整的代码

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

仍然存在疑问，因为您的表达方式不一致关于标量和矢量数量。我强烈建议您尝试解决1D表达的PDE。几乎总是误导。您需要知道什么是标量，什么是向量，这很难在1D中跟踪。