如何将ODE系统与有限元系统相结合

Question

我有一个动态模型，被建立为一个(刚性)的ODEs系统。目前，我使用CVODE (来自Assimulo包中的SUNDIALS包)来解决这个问题，所有这些都是好的。

我现在想添加一个新的三维散热器(与温度相关的热参数)来解决这个问题。我的想法不是为3D热方程从头到尾写出所有方程，而是使用现有的FEM或FVM框架向我提供一个接口，使我能够轻松地将3D块的(t, y)提供给一个例程，并返回剩余的y'。其原理是使用有限元系统中的方程，而不是求解者。CVODE可以利用稀疏性，但是组合系统的求解速度要比有限元系统自己解决的要慢，是针对这种情况而量身定做的。

# pseudocode of a residuals function for CVODE
def residual(t, y):

    # ODE system of n equations 
    res[0] = <function of t,y>;
    res[1] = <function of t,y>;
    ...
    res[n] = <function of t,y>;

    # Here we add the FEM/FVM residuals
    for i in range(FEMcount):
        res[n+1+i] = FEMequations[FEMcount](t,y)

    return res

我的问题是：(a)这种方法是否明智；(b)是否有一个有限元法或FVM库，可以很容易地让我把它当作一个方程组，以便我可以将它“按”到我现有的一组ODE方程。

如果不能让两个系统共享相同的时间轴，那么我将不得不以步进模式运行它们，其中我运行一个模型很短时间，更新另一个模型的边界条件，运行那个模型，更新第一个模型的BCs，等等。

我对优秀的库FiPy有一些经验，我期望最终以上述方式使用该库。但我想知道在这种性质的问题上与其他系统打交道的经验，以及我错过的其他方法。

编辑:我现在有了一些示例代码，它显示了如何使用CVODE解决FiPy网格扩散残值。然而，这只是一种方法(使用FiPy)，其余的问题和关注仍然存在。欢迎任何建议。

from fipy import *
from fipy.solvers.scipy import DefaultSolver
solverFIPY = DefaultSolver()

from assimulo.solvers import CVode as solverASSIMULO
from assimulo.problem import Explicit_Problem as Problem

# FiPy Setup - Using params from the Mesh1D example
###################################################
nx = 50; dx = 1.; D = 1.
mesh = Grid1D(nx = nx, dx = dx)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
valueLeft, valueRight = 1., 0.

phi.constrain(valueRight, mesh.facesRight)
phi.constrain(valueLeft, mesh.facesLeft)

# Instead of eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D),
# Rather just operate on the diffusion term. CVODE will calculate the
# Transient side
edt = ExplicitDiffusionTerm(coeff=D)

timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100

# For comparison with an analytical solution - again,
# taken from the Mesh1D.py example
phiAnalytical = CellVariable(name="analytical value", mesh=mesh)
x = mesh.cellCenters[0]
t = timeStepDuration * steps
from scipy.special import erf
phiAnalytical.setValue(1 - erf(x / (2 * numerix.sqrt(D * t))))
if __name__ == '__main__':
    viewer = Viewer(vars=(phi, phiAnalytical))#, datamin=0., datamax=1.)
    viewer.plot()

raw_input('Press a key...')

# Now for the Assimulo/Sundials solver setup
############################################

def residual(t, X):
    # Pretty straightforward, phi is the unknown
    phi.value = X # This is a vector, 50 elements
    # Can immediately return the residuals, CVODE sees this vector
    # of 50 elements as X'(t), which is like TransientTerm() from FiPy
    return edt.justResidualVector(var=phi, solver=solverFIPY)

x0 = phi.value
t0 = 0.
model = Problem(residual, x0, t0)
simulation = solverASSIMULO(model)
tfinal = steps * timeStepDuration # s,

cell_tol = [1.0e-8]*50
simulation.atol = cell_tol
simulation.rtol = 1e-6
simulation.iter = 'Newton'

t, x = simulation.simulate(tfinal, 0)

print x[-1]
# Write back the answer to compare
phi.value = x[-1]
viewer.plot()
raw_input('Press a key...')

这将产生一个显示完美匹配的图表：

​

​

Answer 1

ODE是一维微分方程。

有限元模型是针对空间问题，即高维问题建立的。你想要一个有限差分法。从一个来自ODE世界的人的角度来看，它更容易解决和理解。

我认为您真正应该问的问题是，如何将您的ODE转换到一个三维问题空间。

多维偏微分方程是很难解决的，但我将介绍一种CFD方法来做这件事，但是只在2D中。http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/

你得花上整整一个下午的时间才能渡过难关！祝好运!