用python求解一维水分运移方程

Question

我试图用python求解一维水分传输方程，但结果并不稳定。也许你们可以帮我。 我有以下圆柱形样品。外壳表面不透水，因此水分只能通过前表面进入： Experimental setup

方程式如下： Used Equation

对于偏方程的空间离散化，使用了有限体积技术。时间离散化采用隐式公式： Explaining solving schemeMy 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, 108.87296769) # adjusted free water saturation (a_ad, a_de)
    b = (1.55379379, 2.4556701)  # correction factor (b_ad, 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

我对每一个暗示都感到高兴。 我目前的假设是：

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

谢谢大家!