【课程】03 Richards方程数值解 |
您所在的位置：网站首页 › ola方程 › 【课程】03 Richards方程数值解 |
【课程】03 Richards方程数值解
|
Z Z Z方向Richards方程 ∂ θ ∂ t = ∂ ∂ z [ D ( θ ) ∂ θ ∂ z − K ( θ ) ] \frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z}[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)] ∂t∂θ=∂z∂[D(θ)∂z∂θ−K(θ)] 式中， D ( θ ) D(\theta) D(θ)为扩散率， K ( θ ) K(\theta) K(θ)为导水率 1.2 差分格式∂ θ ∂ z = θ i + 1 2 j + 1 − θ i − 1 2 j + 1 Δ z \frac{\partial \theta}{\partial z} = \frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z} ∂z∂θ=Δzθi+21j+1−θi−21j+1 ∂ θ ∂ t = θ i j + 1 − θ i j Δ t \frac{\partial \theta}{\partial t} = \frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} ∂t∂θ=Δtθij+1−θij 1.3 方程离散θ i j + 1 − θ i j Δ t = [ D ( θ ) ∂ θ ∂ z − K ( θ ) ] i + 1 2 j + 1 − [ D ( θ ) ∂ θ ∂ z − K ( θ ) ] i − 1 2 j + 1 Δ z \frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} = \frac{[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i+\frac{1}{2}}^{j+1} - [D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i-\frac{1}{2}}^{j+1}}{\Delta z} Δtθij+1−θij=Δz[D(θ)∂z∂θ−K(θ)]i+21j+1−[D(θ)∂z∂θ−K(θ)]i−21j+1 [ D ( θ ) ∂ θ ∂ z − K ( θ ) ] i + 1 2 j + 1 = D ( θ ) i + 1 2 j + 1 ( θ i + 1 2 j + 1 − θ i − 1 2 j + 1 Δ z ) i + 1 2 j + 1 − K ( θ ) i + 1 2 j + 1 = D ( θ ) i + 1 2 j + 1 θ i + 1 j + 1 − θ i j + 1 Δ z − K ( θ ) i + 1 2 j + 1 [D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i+\frac{1}{2}}^{j+1} = D(\theta)_{i+\frac{1}{2}}^{j+1} (\frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z})_{i+\frac{1}{2}}^{j+1}- K(\theta)_{i+\frac{1}{2}}^{j+1} = D(\theta)_{i+\frac{1}{2}}^{j+1} \frac{\theta_{i+1}^{j+1}-\theta_{i}^{j+1}}{\Delta z}- K(\theta)_{i+\frac{1}{2}}^{j+1} [D(θ)∂z∂θ−K(θ)]i+21j+1=D(θ)i+21j+1(Δzθi+21j+1−θi−21j+1)i+21j+1−K(θ)i+21j+1=D(θ)i+21j+1Δzθi+1j+1−θij+1−K(θ)i+21j+1 [ D ( θ ) ∂ θ ∂ z − K ( θ ) ] i − 1 2 j + 1 = D ( θ ) i − 1 2 j + 1 ( θ i + 1 2 j + 1 − θ i − 1 2 j + 1 Δ z ) i − 1 2 j + 1 − K ( θ ) i − 1 2 j + 1 = D ( θ ) i − 1 2 j + 1 θ i j + 1 − θ i − 1 j + 1 Δ z − K ( θ ) i − 1 2 j + 1 [D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i-\frac{1}{2}}^{j+1} = D(\theta)_{i-\frac{1}{2}}^{j+1} (\frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z})_{i-\frac{1}{2}}^{j+1}- K(\theta)_{i-\frac{1}{2}}^{j+1} = D(\theta)_{i-\frac{1}{2}}^{j+1} \frac{\theta_{i}^{j+1}-\theta_{i-1}^{j+1}}{\Delta z}- K(\theta)_{i-\frac{1}{2}}^{j+1} [D(θ)∂z∂θ−K(θ)]i−21j+1=D(θ)i−21j+1(Δzθi+21j+1−θi−21j+1)i−21j+1−K(θ)i−21j+1=D(θ)i−21j+1Δzθij+1−θi−1j+1−K(θ)i−21j+1 则 θ i j + 1 − θ i j Δ t = 1 Δ z [ D ( θ ) i + 1 2 j + 1 θ i + 1 j + 1 − θ i j + 1 Δ z − K ( θ ) i + 1 2 j + 1 − D ( θ ) i − 1 2 j + 1 θ i j + 1 − θ i − 1 j + 1 Δ z + K ( θ ) i − 1 2 j + 1 ] \frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} = \frac{1}{\Delta z}[D(\theta)_{i+\frac{1}{2}}^{j+1} \frac{\theta_{i+1}^{j+1}-\theta_{i}^{j+1}}{\Delta z}- K(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1} \frac{\theta_{i}^{j+1}-\theta_{i-1}^{j+1}}{\Delta z}+ K(\theta)_{i-\frac{1}{2}}^{j+1}] Δtθij+1−θij=Δz1[D(θ)i+21j+1Δzθi+1j+1−θij+1−K(θ)i+21j+1−D(θ)i−21j+1Δzθij+1−θi−1j+1+K(θ)i−21j+1] Δ t D ( θ ) i − 1 2 j + 1 Δ z 2 θ i − 1 j + 1 + { 1 + Δ t Δ z 2 [ D ( θ ) i + 1 2 j + 1 − D ( θ ) i − 1 2 j + 1 ] } θ i j + 1 − Δ t D ( θ ) i + 1 2 j + 1 Δ z 2 θ i + 1 j + 1 = θ i j − Δ t Δ z K ( θ ) i + 1 2 j + 1 + Δ t Δ z K ( θ ) i − 1 2 j + 1 \frac{\Delta t D(\theta)_{i-\frac{1}{2}}^{j+1}}{{\Delta z}^2}\theta_{i-1}^{j+1}+\{1+\frac{\Delta t}{{\Delta z}^2}[D(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1}]\}\theta_i^{j+1}-\frac{\Delta t D(\theta)_{i+\frac{1}{2}}^{j+1}}{{\Delta z}^2}\theta_{i+1}^{j+1} = \theta_i^j-\frac{\Delta t}{\Delta z}K(\theta)_{i+\frac{1}{2}}^{j+1}+\frac{\Delta t}{\Delta z}K(\theta)_{i-\frac{1}{2}}^{j+1} Δz2ΔtD(θ)i−21j+1θi−1j+1+{1+Δz2Δt[D(θ)i+21j+1−D(θ)i−21j+1]}θij+1−Δz2ΔtD(θ)i+21j+1θi+1j+1=θij−ΔzΔtK(θ)i+21j+1+ΔzΔtK(θ)i−21j+1 即 A i θ i − 1 j + 1 + B i θ i j + 1 + C i θ i + 1 j + 1 = D i A_i\theta_{i-1}^{j+1} + B_i\theta_{i}^{j+1}+C_i\theta_{i+1}^{j+1}=D_i Aiθi−1j+1+Biθij+1+Ciθi+1j+1=Di 式中 A i = Δ t Δ z 2 D ( θ ) i − 1 2 j + 1 A_i = \frac{\Delta t }{{\Delta z}^2}D(\theta)_{i-\frac{1}{2}}^{j+1} Ai=Δz2ΔtD(θ)i−21j+1 B i = 1 + Δ t Δ z 2 [ D ( θ ) i + 1 2 j + 1 − D ( θ ) i − 1 2 j + 1 ] B_i = 1+\frac{\Delta t}{{\Delta z}^2}[D(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1}] Bi=1+Δz2Δt[D(θ)i+21j+1−D(θ)i−21j+1] C i = − Δ t Δ z 2 D ( θ ) i + 1 2 j + 1 C_i = -\frac{\Delta t }{{\Delta z}^2}D(\theta)_{i+\frac{1}{2}}^{j+1} Ci=−Δz2ΔtD(θ)i+21j+1 D i = θ i j − Δ t Δ z [ K ( θ ) i + 1 2 j + 1 − K ( θ ) i − 1 2 j + 1 ] D_i = \theta_i^j-\frac{\Delta t}{\Delta z}[K(\theta)_{i+\frac{1}{2}}^{j+1}-K(\theta)_{i-\frac{1}{2}}^{j+1}] Di=θij−ΔzΔt[K(θ)i+21j+1−K(θ)i−21j+1] 采用半隐式解法，用 j j j时刻的 D D D和 K K K代替 j + 1 j+1 j+1时刻的值，采用上游迁移法或 K i ± 1 2 = K i + K i ± 1 2 K_{i \pm \frac{1}{2}}=\frac{K_i+K_{i\pm1}}{2} Ki±21=2Ki+Ki±1 D i ± 1 2 = D i + D i ± 1 2 D_{i \pm \frac{1}{2}}=\frac{D_i+D_{i\pm1}}{2} Di±21=2Di+Di±1 2 模型条件 2.1 边界条件上下边界可以是土壤含水率固定值或时间序列 θ ( t ) \theta(t) θ(t)，在三对角矩阵的第一行和最后一行添加。 2.2 初始条件θ x 0 = θ 0 \theta_x^0 = \theta_0 θx0=θ0 2.3 K ( θ ) K(\theta) K(θ)非饱和土壤导水率 K K K与土壤水吸力 s s s、或与土壤含水率 θ \theta θ的关系，常根据实测结果拟合为经验公式。常用的经验公式形式有 K ( θ ) = K t ( θ θ s ) m 或 K ( θ ) = K t ( θ − θ 0 θ s − θ 0 ) m K(\theta) = K_t(\frac{\theta}{\theta_s})^m 或 K(\theta)=K_t(\frac{\theta - \theta_0}{\theta_s - \theta_0})^m K(θ)=Kt(θsθ)m或K(θ)=Kt(θs−θ0θ−θ0)m K ( θ ) = K t e − β ( θ s − θ ) K(\theta) = K_t e^{-\beta(\theta_s-\theta)} K(θ)=Kte−β(θs−θ) 式中， K t K_t Kt为饱和导水率， θ s \theta_s θs为饱和含水率， θ 0 \theta_0 θ0为某一特征含水率， m m m为经验常数。 2.4 D ( θ ) D(\theta) D(θ)定义非饱和土壤水的扩散率 D ( θ ) D(\theta) D(θ)为导水率 K ( θ ) K(\theta) K(θ)和比水容量 C ( θ ) C(\theta) C(θ)的比值，即 D ( θ ) = K ( θ ) C ( θ ) = K ( θ ) / d θ d ψ m D(\theta)=\frac{K(\theta)}{C(\theta)}=K(\theta)/\frac{d\theta}{d\psi_m} D(θ)=C(θ)K(θ)=K(θ)/dψmdθ 显然，非饱和土壤水的扩散率 D D D同样是土壤含水率 θ \theta θ的或基质势 ψ m \psi_m ψm的函数。其函数关系必须通过试验测定，常用的经验公式形式有 D ( θ ) = a e b θ D(\theta) = ae^{b\theta} D(θ)=aebθ D ( θ ) = D 0 ( θ θ s ) m D(\theta) = D_0(\frac{\theta}{\theta_s})^m D(θ)=D0(θsθ)m D ( θ ) = D 0 e − β ( θ 0 − θ ) D(\theta) = D_0e^{-\beta(\theta_0-\theta)} D(θ)=D0e−β(θ0−θ) 式中， θ 0 \theta_0 θ0为某一特征含水率， a a a， b b b， D 0 D_0 D0、 m m m、 β \beta β为经验常数，取决于土壤质地和结构。 3 模型求解 3.1 三对角矩阵三对角矩阵 A x = B Ax=B Ax=B形式如下： ( 1 0 0 0 … … … … … … 0 A 2 B 2 C 2 A 3 B 3 C 3 A 4 B 4 C 4 A n − 2 B n − 2 C n − 2 A n − 1 B n − 1 C n − 1 0 … … … … … … 0 0 0 1 ) ⏟ A ( θ 1 θ 2 θ 3 θ 4 θ 5 θ n − 2 θ n − 1 θ n ) ⏟ x = ( θ 1 ( t ) D 2 D 3 D 4 D 5 D n − 2 D n − 1 θ n ( t ) ) ⏟ B \underbrace{\begin{pmatrix} 1&0&0&0&\dots&\dots&\dots&\dots&\dots&\dots&0\\[4pt] A_{2}&B_{2}&C_{2}\\[4pt] & A_{3}&B_{3}&C_{3}\\[4pt] & & A_{4}&B_{4}&C_{4} \\[4pt] \\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & &&&&&& A_{n-2}&B_{n-2}&C_{n-2} \\[4pt] & & &&&&&&A_{n-1}&B_{n-1}&C_{n-1} \\[4pt] 0&\dots &\dots &\dots&\dots&\dots&\dots&0 & 0 &0 &1 \\[4pt] \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} \theta_{1} \\[4pt] \theta_{2} \\[4pt] \theta_{3} \\[4pt] \theta_{4} \\[4pt] \theta_{5} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] \theta_{n-2} \\[4pt] \theta_{n-1} \\[4pt] \theta_{n} \\[4pt] \end{pmatrix}}_{x}=\underbrace{\begin{pmatrix} \theta_{1}(t) \\[4pt] D_{2} \\[4pt] D_{3} \\[4pt] D_{4} \\[4pt] D_{5} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] D_{n-2} \\[4pt] D_{n-1} \\[4pt] \theta_{n}(t) \\[4pt] \end{pmatrix}}_{B} A 1A200B2A3…0C2B3A4…0C3B4……C4………………An−20…Bn−2An−10…Cn−2Bn−100Cn−11 x θ1θ2θ3θ4θ5θn−2θn−1θn =B θ1(t)D2D3D4D5Dn−2Dn−1θn(t) 3.2 各参数取值 参数含义取值Z土柱深度1mT模拟时间1hdz空间步长1cmdt时间步长0.0001h θ 1 \theta_1 θ1上边界土壤含水量0.30 θ s \theta_s θs饱和土壤含水量0.45 θ 0 \theta_0 θ0初始土壤含水量0.17 K t K_t Kt饱和导水率16.72m/h m k m_k mk导水率经验常数13 a a a扩散率经验常数20.019 b b b扩散率经验常数212.3804 3 求解结果 3.1 计算结果
[1]土壤水动力学，1987.3 [2]李映强.赤红壤非饱和土壤水扩散率及其影响因素[J].华南农业大学学报,1998(02):3-5. 作者简介很高兴认识您！ 我叫卢家波，河海大学水文学及水资源博士研究生，研究兴趣为高效洪水淹没预测、洪水灾害预警、机器学习、替代模型和降阶模型。 变化环境下，极端洪水事件多发，我希望能通过研究为水灾害防御做出贡献，为人民服务。 欢迎交流讨论和研究合作，vx Jiabo_Lu。 主页 https://lujiabo98.github.io 简历 https://lujiabo98.github.io/file/CV_JiaboLu_zh.pdf 博客 https://blog.csdn.net/weixin_43012724?type=blog 来信请说明博客标题及链接，谢谢。 |
【本文地址】