使用 FiPy 在一维圆柱坐标系上求解 PDE

Solving PDE on 1D cylindrical coordinates with FiPy

首先让我说一下,我在 NARKIVE FiPy 邮件列表存档中发现了与我类似的问题,但由于方程式无法加载,因此它们不是很有用。例如 Convection-diffusion problem on a 1D cylindrical grid, or on another mailing list archive Re: FiPy Heat Transfer Solution。在第二封 linked 邮件中丹尼尔说:

There are two ways to solve on a cylindrical domain in FiPy. You can either use the standard diffusion equation in Cartesian coordinates (2nd equation below) and with a mesh that is actually cylindrical in shape or you can use the diffusion equation formulated on a cylindrical coordinate system (1st equation below) and use a standard 2D / 1D grid mesh.

方程式不存在。在这种情况下它实际上很好,因为我理解第一个解决方案并且我想使用它。

我想在一维圆柱网格上求解以下方程(抱歉,我还没有 10 个声望,所以我不能 post 漂亮的渲染方程):

有边界条件:

其中 rho_core 是网格的左侧,rho_edge 是网格的右侧。 rho 是归一化半径,J 是雅可比行列式:

R是以米为单位的真实半径,所以雅可比矩阵的量纲是距离。初始条件并不重要,但在我的代码示例中,我将在 R=0.8.

处使用数字 Dirac-delta

我有一个没有(!)雅可比行列式的工作示例,但它很长,而且它没有使用 FiPy 的查看器所以我将 link 一个要点:https://gist.github.com/leferi99/142b90bb686cdf5116ef5aee425a4736

问题主要部分如下:

    import fipy as fp  ## finite volume PDE solver
    from fipy.tools import numerix  ## requirement for FiPy, in practice same as numpy
    import copy  ## we need the deepcopy() function because some FiPy objects are mutable
    import numpy as np
    import math

    ## numeric implementation of Dirac delta function
    def delta_func(x, epsilon, coeff):
        return ((x < epsilon) & (x > -epsilon)) * \
            (coeff * (1 + numerix.cos(numerix.pi * x / epsilon)) / (2 * epsilon))

    rho_from = 0.7  ## normalized inner radius
    rho_to = 1.  ## normalized outer radius
    nr = 1000  ## number of mesh cells
    dr = (rho_to - rho_from) / nr  ## normalized distance between the centers of the mesh cells
    duration = 0.001  ## length of examined time evolution in seconds
    nt = 1000  ## number of timesteps
    dt = duration / nt  ## length of one timestep
    
    ## 3D array for storing the density with the correspondant normalized radius values
    ## the density values corresponding to the n-th timestep will be in the n-th line 
    solution = np.zeros((nt,nr,2))
    ## loading the normalized radial coordinates into the array
    for j in range(nr):
        solution[:,j,0] = (j * dr) + (dr / 2) + rho_from
    
    mesh = fp.CylindricalGrid1D(dx=dr, nx=nr)  ## 1D mesh based on the normalized radial coordinates 
    mesh = mesh + (0.7,)  ## translation of the mesh to rho=0.7
    n = fp.CellVariable(mesh=mesh)  ## fipy.CellVariable for the density solution in each timestep
    
    diracLoc = 0.8  ## location of the middle of the Dirac delta
    diracCoeff = 1.  ## Dirac delta coefficient ("height")
    diracPercentage = 2  ## width of Dirac delta (full width from 0 to 0) in percentage of full examined radius
    diracWidth = int((nr / 100) * diracPercentage)
    
    ## diffusion coefficient
    diffCoeff = fp.CellVariable(mesh=mesh, value=100.)
    ## convection coefficient - must be a vector
    convCoeff = fp.CellVariable(mesh=mesh, value=(1000.,))
    
    ## applying initial condition - uniform density distribution
    n.setValue(1)

    ## boundary conditions
    gradLeft = (0.,)  ## density gradient (at the "left side of the radius") - must be a vector
    valueRight = 0.  ## density value (at the "right end of the radius")
    n.faceGrad.constrain(gradLeft, where=mesh.facesLeft)  ## applying Neumann boundary condition
    n.constrain(valueRight, mesh.facesRight)  ## applying Dirichlet boundary condition
    convCoeff.setValue(0, where=mesh.x<(R_from + dr))  ## convection coefficient 0 at the inner edge
    
    ## the PDE
    eq = (fp.TransientTerm() == fp.DiffusionTerm(coeff=diffCoeff)
          - fp.ConvectionTerm(coeff=convCoeff))
    
    ## Solving the PDE and storing the data
    for i in range(nt):
        eq.solve(var=n, dt=dt)
        solution[i,0:nr,1]=copy.deepcopy(n.value)

我的代码可以使用与上述相同的边界条件求解以下方程:

为了简单起见,我使用了空间无关的系数,唯一的例外是在内边缘,对流系数为 0,扩散系数几乎为 0。在 linked 代码中,我使用均匀分布的初始条件。

我的第一个问题是为什么在使用 fipy.Grid1Dfipy.CylindricalGrid1D 时得到完全相同的结果?我应该得到不同的结果,对吧?我应该如何重写我的代码才能区分简单的一维网格和一维圆柱网格?

我的实际问题不在于这个确切的代码,我只是想简化我的问题,但正如评论中所指出的那样,这段代码不会对不同的网格产生相同的结果。所以我只会 post 一个 GitHub link 到一个 Jupyter Notebook,它可能会在未来停止工作。 The Jupyter Notebook 如果你想 运行 它,第一个代码单元格应该是 运行,然后只有最后一个单元格是相关的。忽略参考图像。线图显示扩散和对流系数。当我 运行 最后一个带有 Grid1DCylindricalGrid1D 的单元格时,我得到了相同的结果(我非常精确地比较了这些图)

抱歉,我无法重命名我的所有变量,所以我希望根据我的评论和上面更改的代码(我也更改了代码中的注释)您可以理解我正在尝试做什么。

我的另一个问题是关于雅可比矩阵的。我该如何实施?我查看了文档中唯一使用雅可比行列式的示例,但该雅可比行列式是一个矩阵,并且它还使用了 scipy.optimize.fsolve() 函数。

[根据评论中的讨论拼凑答案]

a Grid1D 和 a CylindricalGrid1D 的结果相似,尤其是在早期步骤中,但它们并不相同。随着问题的发展,它们是完全不同的。

FiPy 不喜欢散度以外的东西,但你应该可以将方程乘以 J 并将其放入 TransientTerm 的系数中,例如,

eq = fp.TransientTerm(J) == fp.DiffusionTerm(coeff=J * diffCoeff) - fp.ConvectionTerm(coef=J * convCoeff)

对于雅可比矩阵,您可以根据归一化半径为实际半径创建 CellVariable,然后取其梯度:

real_radius = fp.CellVariable(mesh=mesh, value=...)
J = real_radius.grad.dot([[1]])

.grad returns一个向量,即使是一维的,但系数必须是标量,所以取点积得到x分量。