问题描述
我正在尝试用 python 求解一维水分输送方程,但我的结果不稳定。或许你们可以帮帮我。 我有以下圆柱形样品。外壳表面不透水和蒸汽,因此水分只能从前表面进入: Experimental setup
等式如下: Used Equation
对于部分方程的空间离散化,使用了有限体积技术。时间离散是通过隐式公式完成的: 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
- 也许我无法像以前那样实现边界条件
- 也许我的求解类型不适用于水相关系数
谢谢!
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)