使用SYSPY计算速度Jacobian

Question

我想在多孔介质中获得对流方程的通量雅各布的分析表达，这在执行计算流体动力学时很有用。我以为我可以轻松地在sympy做到这一点，但我被困在某个地方。

这是我的方程式：

，

phi呈孔隙率，c元素的浓度和v流体速度。

在我的情况下，jacobian的通量很简单：

我可以使用darcy的定律来计算v：

具有k的渗透性，nf流体粘度和pf流体压力。

渗透性取决于这种非常简单的关系：

通过做这种非常简单的变量更改，我可以简化我的问题：

”>

alt = “ 我尝试以这种方式解决的方法：

from sympy import * 

# Create symbols
k0, phi, n, Pn, Pc, Ps, C = symbols('k0 phi n P_{i+1} P_{i} P_{i-1} C')
C_tild = phi * C
dx = Symbol(r"\Delta x")
nf = Symbol(r"\mu_f")

# Darcy's law while using finite differences for the Gradient operator in the pressure term
qf = (k0 * phi**n / nf) * (((Pn + Pc)/2 - (Ps + Pc)/2) / dx)
vf = qf / phi

flux = vf * C_tild
display(flux)
derivative = diff(flux, C_tild)

此返回错误，因为我试图同时使用两个变量（PHI和C）进行衍生物。 我在这里做错了什么？我如何使Sympy了解变量的变化？还是可以用更简单的表达来重写雅各布？

编辑

错误如下：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/var/folders/h6/3p69g_d146x78r6b75_h_td40000gp/T/ipykernel_64334/485359947.py in <module>
     11 flux = vf * C_tild
     12 display(flux)
---> 13 derivative = diff(flux, C_tild)

~/opt/anaconda3/envs/Two-phase_flow/lib/python3.9/site-packages/sympy/core/function.py in diff(f, *symbols, **kwargs)
   2497     """
   2498     if hasattr(f, 'diff'):
-> 2499         return f.diff(*symbols, **kwargs)
   2500     kwargs.setdefault('evaluate', True)
   2501     return _derivative_dispatch(f, *symbols, **kwargs)

~/opt/anaconda3/envs/Two-phase_flow/lib/python3.9/site-packages/sympy/core/expr.py in diff(self, *symbols, **assumptions)
   3526     def diff(self, *symbols, **assumptions):
   3527         assumptions.setdefault("evaluate", True)
-> 3528         return _derivative_dispatch(self, *symbols, **assumptions)
   3529 
   3530     ###########################################################################

~/opt/anaconda3/envs/Two-phase_flow/lib/python3.9/site-packages/sympy/core/function.py in _derivative_dispatch(expr, *variables, **kwargs)
   1921         from sympy.tensor.array.array_derivatives import ArrayDerivative
   1922         return ArrayDerivative(expr, *variables, **kwargs)
-> 1923     return Derivative(expr, *variables, **kwargs)
   1924 
   1925 

~/opt/anaconda3/envs/Two-phase_flow/lib/python3.9/site-packages/sympy/core/function.py in __new__(cls, expr, *variables, **kwargs)
   1346             if not v._diff_wrt:
   1347                 __ = ''  # filler to make error message neater
-> 1348                 raise ValueError(filldedent('''
   1349                     Can't calculate derivative wrt %s.%s''' % (v,
   1350                     __)))

ValueError: 
Can't calculate derivative wrt C*phi.

Answer 1

仅供参考，这是isympy会话中的代码，

In [7]: k0, phi, n, Pn, Pc, Ps, C = symbols('k0 phi n P_{i+1} P_{i} P_{i-1} C')
   ...: C_tild = phi * C
   ...: dx = Symbol(r"\Delta x")
   ...: nf = Symbol(r"\mu_f")
   ...: 
   ...: # Darcy's law while using finite differences for the Gradient operator i
   ...: n the pressure term
   ...: qf = (k0 * phi**n / nf) * (((Pn + Pc)/2 - (Ps + Pc)/2) / dx)
   ...: vf = qf / phi
   ...: 
   ...: flux = vf * C_tild

In [8]: flux
Out[8]: 
      n ⎛P_{i+1}   P_{i-1}⎞
C⋅k₀⋅φ ⋅⎜─────── - ───────⎟
        ⎝   2         2   ⎠
───────────────────────────
       \Delta x⋅\mu_f      

In [9]: C_tild
Out[9]: C⋅φ

因为c是一个简单的乘数，我认为它不会为计算添加任何内容。只需使用flux没有c，才能获得：

In [18]: diff(vf*phi, phi)
Out[18]: 
      n ⎛P_{i+1}   P_{i-1}⎞
k₀⋅n⋅φ ⋅⎜─────── - ───────⎟
        ⎝   2         2   ⎠
───────────────────────────
      \Delta x⋅\mu_f⋅φ