——最近感覺生活陷入了低谷，已迷失，不知道目前這兒的生活意義在哪，把握不清楚未來走勢。當然目前的重大挫折也會帶給自己成長，以後學會自我調節、耐心處理困難才能有所成就、有所收穫。
馬上要到520了，祝大家當天過得愉快！
上一篇寫了一個簡單的非飽和土壤下滲的代碼，有些問題，需要勘誤下，擴散係數D，忘了加進去，雖然上一篇取的是1，並不影響計算結果。但是如果更改D的值，就完全沒效果了。因為變量D並沒有加到代碼中。
這次進階篇，我們把含水率θ當做變量，利用土壤下滲方程，計算土壤非飽和帶的含水率，然後再計算濃度。
其實跟水質模擬過程基本上一樣的，計算出水動力，再計算水質。有水質模擬興趣的可以拿此代碼的思路去試試模擬水質過程。
溶質下滲公式依據《環境影響評價技術導則 土壤環境》（2018）：
含水率θ公式為：
然後一切從簡，各種參數均為定值，變量僅為c、θ（懶人，簡單來做--，
）開始貼代碼了：
離散後形式：
c[i]= c_α*cn[i+1]+c_β*cn[i]+c_γ*cn[i-1];n為時間、i為土壤剖分網格
θ[i]= θ_α*θn[i+1]+θ_β*θn[i]+θ_γ*θn[i-1];n為時間、i為土壤剖分網格
1# -*- coding: utf-8 -*- 2""" 3Created on Sat May 9 16:08:40 2020 4@author: zz 5""" 6import numpy as np 7import matplotlib.pyplot as plt 8from pylab import * 9rcParams['font.sans-serif']=['SimHei']10D=0.5; #彌散係數,m2/s11q=1; #滲透速率,12Z=10;#沿Z軸的距離，m13nx=10; #縱向網格個數14nt=500;#時間步數15dx=Z/nx;#j為距離，m16dt=0.01;#i為時間步長，s17i=0;#i為土壤剖面網格編號數，s18c = np.zeros(nx)#初始條件19#土壤含水率的邊界等設置情況20θ= np.ones(nx)#給出一個數組存儲不同時間步長的數據21θ[0] = 20#上邊界22θ[nx-1]=10#下邊界23θ_α=D*dt/dx**2;#合併後的係數24θ_β=1-2*D*dt/dx**2;#合併後的係數25θ_γ=D*dt/dx**2;#合併後的係數26#土壤濃度的一些邊角條件27c[0] = 10#上邊界28c[nx-1]=5#下邊界,其實下邊界並不是必需的，但是還是給了。29cn = np.ones(nx)#給出一個數組存儲不同時間步長的數據30c_α=np.zeros(nx)31c_β=np.zeros(nx)32c_α[i]=D*dt/dx**2-q/θ[i]*dt/dx;#合併後的係數33c_β[i]=1+q/θ[i]*dt/dx-2*D*dt/dx**2;#合併後的係數34c_γ=D*dt/dx**2;#合併後的係數35#求土壤含水量36for n in range(1,nt):37 θn = θ.copy() ##複製循環數組38for i in range(1,nx-1):39 θ[i]= θ_α*θn[i+1]+θ_β*θn[i]+θ_γ*θn[i-1];40plt.plot(θ,label="含水率")41# plt.title("θ(%)")42plt.xlabel('X(m)')43# plt.ylabel('θ(%)') 44plt.legend()45#求土壤濃度46for n in range(1,nt):47 cn = c.copy() ##複製循環數組48for i in range(1,nx-1):49 c_α[i]=D*dt/dx**2-q/θ[i]*dt/dx;#合併後的係數50 c_β[i]=1+q/θ[i]*dt/dx-2*D*dt/dx**2;#合併後的係數51 c[i]= c_α[i]*cn[i+1]+c_β[i]*cn[i]+c_γ*cn[i-1];52plt.plot(c,linestyle='--',label="汙染物濃度")53# plt.title("C(mg/l)")54# plt.xlabel('X(m)')55# plt.ylabel('C(mg/l)')56plt.legend()57plt.show() 5859
1：計算結果：
2：大家可以試試代碼，運算不同的結果！
3：可以增加擴散係數的非常量過程等，自己玩玩，深入一點！