用python求解一维水分输送方程

如何解决用python求解一维水分输送方程

我正在尝试用 python 求解一维水分输送方程，但我的结果不稳定。或许你们可以帮帮我。我有以下圆柱形样品。外壳表面不透水和蒸汽，因此水分只能从前表面进入： Experimental setup

对于部分方程的空间离散化，使用了有限体积技术。时间离散是通过隐式公式完成的： Explaining solving scheme My solution

这是我的python代码：

import numpy as np
import math
import scipy.linalg 
import matplotlib.pyplot as plt
import functionCollection as fC
import time
import random
"""
To do:
1) Define and update parameters
2) Initial and boundary conditions
3) Create the tridiag
4) Solve system of linear equations
After 1-4 works,do ADI method to update parameters in tridiag with step = 1/2(phi_n + phi_n-1)
"""
#1) Define the parameters
probeLength = 0.15
spatialSteps = 23
timePeriod = 3600*24*36
timesteps = 24*36
# create meshpoints in space
x = np.linspace(0,probeLength,spatialSteps+1)
# define stepsize
dx = x[1]-x[0]
# create meshpoints in time
t = np.linspace(0,timePeriod,timesteps+1)
# define stepsize
dt = t[1]-t[0]

# define initial and boundary conditions [-]
ic_phi = 0.5
bc_phi = 0.9
# define vectors for humidity
phi_n = phi_n1 = np.zeros(spatialSteps+2)
# fill ibc
phi_n1 = phi_n1+ic_phi
phi_n1[0] = phi_n1[-1] = bc_phi

# define fixed parameters 
p_sat = 2814.6337703571003
delta_pw = 6.579*10**-11

# define list of results for humidity and water content
lstPhi = []
lstW = []
lstPhi.append(phi_n1)
# for loop to solve for timesteps
for i in range(0,timesteps+1):
    
    # update variable parameters
    D_phi = fC.calc_D_phi(phi_n1)
    dw_dphi = fC.derivative_dw_dphi(phi_n1)
    # calculate F
    F = (D_phi + delta_pw * p_sat) * (1 / dw_dphi ) * (dt / (dx**2))
    b = phi_n1
    b[0] = b[-1] = bc_phi
    tridiag = fC.make_tridiag(F,spatialSteps)
    phi_n = scipy.linalg.solve(tridiag,b)

    lstPhi.append(phi_n)
    phi_n1 = phi_n

plt.figure()
for phi in lstPhi:
    plt.plot(phi)
plt.show()

我的辅助函数是：

import numpy as np
import math


def rf_to_h2o_content(phi,is_increasing=True):
    """
    Maps relative humidity [-] to water content in [kg/m³] using the curve fitting described in Künzel et. al.
    correction factor: b_ad = 1.55379379,b_de = 2.4556701
    adjusted free water saturation only for fitting "Sorptionsisotherme": a_ad = 109.60291049,a_de = 108.87296769
    """
    a = (109.60291049,108.87296769) # adjusted free water saturation (a_ad,a_de)
    b = (1.55379379,2.4556701)  # correction factor (b_ad,b_de)
    if is_increasing:
        w = a[0]*(b[0]-1)*np.array(phi)/(b[0]-np.array(phi))
    else:
        w = a[1]*(b[1]-1)*np.array(phi)/(b[1]-np.array(phi))


    return w


def derivative_dw_dphi(phi,is_increasing=True):
    """
    Derivative of "Feuchtespeicherfunktion" [kg/m³].
    Returns the gradient for dw/dphi at given phi.
    """
    a = (109.60291049,b_de)
    if is_increasing:
        dw_dphi = a[0]*((b[0]-1) * (b[0]-phi) + (b[0]-1)*phi) / (b[0]-phi)**2
    else:
        dw_dphi = a[1]*((b[1]-1) * (b[1]-phi) + (b[1]-1)*phi) / (b[1]-phi)**2


    return dw_dphi




def _calc_D_ws(w):
    """
    Calculate "Kapillartransportkoeffizient für den Saugvorgang" [m²/s]
    Needs water content of last/current step.
    """
    w_f = 274.175 # free water saturation [kg/m3]
    A = 0.0247126 # Wasseraufnahmekoeffizient [kg/m2s0.5]
    D_ws = 3.8 * (A / w_f)**2 * 1000**(w/w_f-1) 


    return D_ws


def calc_D_phi(phi,is_increasing=True):
    """
    Calculate "Flüssigleitkoeffzient" D_phi = D_w * dw/dphi [kg/ms]
    """
    w = rf_to_h2o_content(phi,is_increasing)
    D_ws = _calc_D_ws(w)
    dw_dphi = derivative_dw_dphi(phi,is_increasing)
    D_phi = D_ws * dw_dphi


    return D_phi


def make_tridiag(F,tridiagDim):
    """tridiagDim excludes the extension."""
    A = np.zeros((tridiagDim+2,tridiagDim+2))
    for i in range(1,tridiagDim+1):
        A[i,i-1] = -F[i-1]
        A[i,i+1] = -F[i+1]
        A[i,i] = 1 + 2*F[i]
    A[0,0] = A[tridiagDim+1,tridiagDim+1] = 1
    return A

我的结果应该是这样的： How the results should look like 但它们看起来像这样： How the wrong results look like at the moment

我对每一个提示都很满意。我目前的假设是：

也许我无法像以前那样实现边界条件
也许我的求解类型不适用于水相关系数

谢谢！