計算傳熱學第一次大作業
1 Taylor級數展開法
1.1 網格劃分
Taylor級數展開發的網格劃分使用外點法
1.2 離散方程表達式
ρ u d ϕ d x = Γ d 2 ϕ d x 2 (1) \rho u\frac{d\phi}{d x}=\Gamma \frac{d^{2}\phi}{dx^{2}}\tag{1} ρudxdϕ=Γdx2d2ϕ(1)
将 ϕ \phi ϕ分别在 i i i+1點和 i i i-1點在 i i i點展開
ϕ ( i + 1 ) = ϕ ( i ) + d ϕ d x ∣ i δ x + d 2 ϕ d x 2 ∣ i δ x 2 2 ! + … … (2) \phi(i+1) = \phi(i)+\left.\frac{d\phi}{dx}\right|_{i}\delta x+\left.\frac{d^{2}\phi}{dx^{2}}\right|_{i}\frac{\delta x^{2}}{2!}+…… \tag{2} ϕ(i+1)=ϕ(i)+dxdϕ∣∣∣∣iδx+dx2d2ϕ∣∣∣∣i2!δx2+……(2)
ϕ ( i − 1 ) = ϕ ( i ) − d ϕ d x ∣ i δ x + d 2 ϕ d x 2 ∣ i δ x 2 2 ! + … … (3) \phi(i-1) = \phi(i)-\left.\frac{d\phi}{dx}\right|_{i}\delta x+\left.\frac{d^{2}\phi}{dx^{2}}\right|_{i}\frac{\delta x^{2}}{2!}+…… \tag{3} ϕ(i−1)=ϕ(i)−dxdϕ∣∣∣∣iδx+dx2d2ϕ∣∣∣∣i2!δx2+……(3)
聯立(2)、(3)式得
d ϕ d x ∣ i = ϕ i + 1 − ϕ i − 1 2 δ x , o ( δ x 2 ) (4) \left.\frac{d\phi}{dx}\right|_{i}=\frac{\phi_{i+1}-\phi_{i-1}}{2\delta x},o(\delta x^{2})\tag{4} dxdϕ∣∣∣∣i=2δxϕi+1−ϕi−1,o(δx2)(4)
同樣的，聯立(2)、(3)式得
d 2 ϕ d x 2 ∣ i = ϕ i + 1 − 2 ϕ i + ϕ i − 1 δ x 2 , o ( δ x 2 ) (5) \left.\frac{d^{2}\phi}{dx^{2}}\right|_{i}=\frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{\delta x^{2}},o(\delta x^{2})\tag{5} dx2d2ϕ∣∣∣∣i=δx2ϕi+1−2ϕi+ϕi−1,o(δx2)(5)
将(4)、(5)式帶入(1)式
ρ u ϕ i + 1 − ϕ i − 1 2 δ x = Γ ϕ i + 1 − 2 ϕ i + ϕ i − 1 δ x 2 (6) \rho u \frac{\phi_{i+1}-\phi_{i-1}}{2\delta x}=\Gamma \frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{\delta x^{2}}\tag{6} ρu2δxϕi+1−ϕi−1=Γδx2ϕi+1−2ϕi+ϕi−1(6)
化簡得
4 Γ ϕ i = ( ρ u δ x + 2 Γ ) ϕ i − 1 + ( 2 Γ − ρ u δ x ) ϕ i + 1 , o ( δ x 2 ) (7) 4\Gamma \phi_{i}=(\rho u \delta x+2\Gamma)\phi_{i-1}+(2\Gamma-\rho u \delta x)\phi_{i+1},o(\delta x^{2})\tag{7} 4Γϕi=(ρuδx+2Γ)ϕi−1+(2Γ−ρuδx)ϕi+1,o(δx2)(7)
2 控制容積積分法
2.1 網格劃分
控制容積積分法使用内節點法劃分網格
2.2 離散方程表達式
ρ u d ϕ d x = Γ d 2 ϕ d x 2 (8) \rho u\frac{d\phi}{d x}=\Gamma \frac{d^{2}\phi}{dx^{2}}\tag{8} ρudxdϕ=Γdx2d2ϕ(8)
将方程在空間内積分
∫ w e ρ u d ϕ d x d x = ∫ w e Γ d 2 ϕ d x 2 d x (9) \int_{w}^{e}\rho u\frac{d\phi}{d x}dx=\int_{w}^{e}\Gamma \frac{d^{2}\phi}{dx^{2}}dx\tag{9} ∫weρudxdϕdx=∫weΓdx2d2ϕdx(9)
其中
∫ w e d ϕ d x d x = ϕ ∣ e − ϕ ∣ w (10) \int_{w}^{e}\frac{d\phi}{d x}dx=\left.\phi\right|_{e}-\left.\phi\right|_{w}\tag{10} ∫wedxdϕdx=ϕ∣e−ϕ∣w(10)
假設 ϕ \phi ϕ對空間呈分段線性變化
ϕ ∣ e − ϕ ∣ w = ϕ E + ϕ P 2 − ϕ P + ϕ W 2 = ϕ E − ϕ W 2 (11) \left.\phi\right|_{e}-\left.\phi\right|_{w}=\frac{\phi_{E}+\phi_{P}}{2}-\frac{\phi_{P}+\phi_{W}}{2}=\frac{\phi_{E}-\phi_{W}}{2}\tag{11} ϕ∣e−ϕ∣w=2ϕE+ϕP−2ϕP+ϕW=2ϕE−ϕW(11)
對于擴散項
∫ w e d 2 ϕ d x 2 d x = d ϕ d x ∣ e − d ϕ d x ∣ w (12) \int_{w}^{e}\frac{d^{2}\phi}{dx^{2}}dx=\left.\frac{d\phi}{dx}\right|_{e}-\left.\frac{d\phi}{dx}\right|_{w}\tag{12} ∫wedx2d2ϕdx=dxdϕ∣∣∣∣e−dxdϕ∣∣∣∣w(12)
假設 d ϕ d x \frac{d\phi}{dx} dxdϕ在空間呈分段線性變化
d ϕ d x ∣ e − d ϕ d x ∣ w = ϕ E − ϕ P ( δ x ) e − ϕ P − ϕ W ( δ x ) w (13) \left.\frac{d\phi}{dx}\right|_{e}-\left.\frac{d\phi}{dx}\right|_{w}=\frac{\phi_{E}-\phi_{P}}{(\delta x)_{e}}-\frac{\phi_{P}-\phi_{W}}{(\delta x)_{w}}\tag{13} dxdϕ∣∣∣∣e−dxdϕ∣∣∣∣w=(δx)eϕE−ϕP−(δx)wϕP−ϕW(13)
若使用均分網格
= ϕ E − ϕ P Δ x − ϕ P − ϕ W Δ x (14) ~~~~~~~~~~~~~~~~~~~~~~~~~=\frac{\phi_{E}-\phi_{P}}{\Delta x}-\frac{\phi_{P}-\phi_{W}}{\Delta x}\tag{14} =ΔxϕE−ϕP−ΔxϕP−ϕW(14)
将式(11)、(14)帶入(8)
ρ u ( ϕ E + ϕ P 2 − ϕ P + ϕ W 2 ) = Γ ( ϕ E − ϕ P Δ x − ϕ P − ϕ W Δ x ) (15) \rho u (\frac{\phi_{E}+\phi_{P}}{2}-\frac{\phi_{P}+\phi_{W}}{2})=\Gamma (\frac{\phi_{E}-\phi_{P}}{\Delta x}-\frac{\phi_{P}-\phi_{W}}{\Delta x})\tag{15} ρu(2ϕE+ϕP−2ϕP+ϕW)=Γ(ΔxϕE−ϕP−ΔxϕP−ϕW)(15)
整理得
4 Γ ϕ P = ( 2 Γ − ρ u Δ x ) ϕ E + ( 2 Γ + ρ u Δ x ) ϕ W (16) 4\Gamma \phi_{P}=(2\Gamma-\rho u \Delta x)\phi_{E}+(2\Gamma+\rho u \Delta x)\phi_{W}\tag{16} 4ΓϕP=(2Γ−ρuΔx)ϕE+(2Γ+ρuΔx)ϕW(16)
3. Gauss-Seidel疊代法
原問題的代數方程
4 Γ ϕ i = ( ρ u δ x + 2 Γ ) ϕ i − 1 + ( 2 Γ − ρ u δ x ) ϕ i + 1 4\Gamma \phi_{i}=(\rho u \delta x+2\Gamma)\phi_{i-1}+(2\Gamma-\rho u \delta x)\phi_{i+1} 4Γϕi=(ρuδx+2Γ)ϕi−1+(2Γ−ρuδx)ϕi+1
令 α = ρ u Γ \alpha = \frac{\rho u}{\Gamma} α=Γρu則
4 ϕ i = ( α δ x + 2 ) ϕ i − 1 + ( 2 − α δ x ) ϕ i + 1 (17) 4 \phi_{i}=(\alpha\delta x+2)\phi_{i-1}+(2-\alpha \delta x)\phi_{i+1}\tag{17} 4ϕi=(αδx+2)ϕi−1+(2−αδx)ϕi+1(17)
假設有 N N N個格點，則該方程得Guass-Seidel疊代格式
{ ϕ 1 k + 1 = 1 4 [ ( α δ x + 2 ) ϕ 0 + ( 2 − α δ x ) ϕ 2 k ] ϕ 2 k + 1 = 1 4 [ ( α δ x + 2 ) ϕ 1 k + 1 + ( 2 − α δ x ) ϕ 3 k ] ⋮ ϕ N k + 1 = 1 4 [ ( α δ x + 2 ) ϕ N − 1 k + 1 + ( 2 − α δ x ) ϕ N + 1 ] \left\{ \begin{aligned} \phi_{1}^{k+1} & = &\frac{1}{4}[(\alpha\delta x+2)\phi_{0}+(2-\alpha \delta x)\phi_{2}^{k}] \\ \phi_{2}^{k+1} & = & \frac{1}{4}[(\alpha\delta x+2)\phi_{1}^{k+1}+(2-\alpha \delta x)\phi_{3}^{k}] \\ \vdots\\ \phi_{N}^{k+1} & = & \frac{1}{4}[(\alpha\delta x+2)\phi_{N-1}^{k+1}+(2-\alpha \delta x)\phi_{N+1}] \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧ϕ1k+1ϕ2k+1⋮ϕNk+1===41[(αδx+2)ϕ0+(2−αδx)ϕ2k]41[(αδx+2)ϕ1k+1+(2−αδx)ϕ3k]41[(αδx+2)ϕN−1k+1+(2−αδx)ϕN+1]
寫為矩陣形式：
( − 4 2 − α δ x 2 + α δ x ⋱ 2 − α δ x 2 + α δ x − 4 ) ( ϕ 1 ϕ 2 ϕ 3 ⋮ ϕ N ) ≐ ( − ( 2 + α δ x ) ϕ 0 0 0 ⋮ − ( 2 − α δ x ) ϕ N + 1 ) (18) \begin{pmatrix} -4 & 2-\alpha \delta x & {} & {} & {} \\ 2+\alpha \delta x & {} & {} & {} & {} \\ {} & {} & \ddots & {} & {} \\ {} & {} & {} & {} & 2-\alpha \delta x \\ {} & {} & {} & 2+\alpha \delta x & -4 \\ \end{pmatrix} \begin{pmatrix} {{\phi }_{1}} \\ {{\phi }_{2}} \\ {{\phi }_{3}} \\ \vdots \\ {{\phi }_{N}} \\ \end{pmatrix} \doteq \begin{pmatrix} -(2+\alpha \delta x){{\phi }_{0}} \\ 0 \\ 0 \\ \vdots \\ -(2-\alpha \delta x){{\phi }_{N+1}} \\ \end{pmatrix}\tag{18} ⎝⎜⎜⎜⎜⎛−42+αδx2−αδx⋱2+αδx2−αδx−4⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛ϕ1ϕ2ϕ3⋮ϕN⎠⎟⎟⎟⎟⎟⎞≐⎝⎜⎜⎜⎜⎜⎛−(2+αδx)ϕ000⋮−(2−αδx)ϕN+1⎠⎟⎟⎟⎟⎟⎞(18)
解得方程在 α = ( − 5 , − 1 , 0 , 1 , 5 ) \alpha=(-5,-1,0,1,5) α=(−5,−1,0,1,5)時的曲線
4. TDMA方法
對于任意一個三對角陣
A i T i = B i T i + 1 + C i T i − 1 + D i (19) A_{i}T_{i}=B_{i}T_{i+1}+C_{i}T_{i-1}+D_{i}\tag{19} AiTi=BiTi+1+CiTi−1+Di(19)
我們都可以将其轉化至
T i − 1 = P i − 1 T i + Q i − 1 (20) T_{i-1}=P_{i-1}T_{i}+Q_{i-1}\tag{20} Ti−1=Pi−1Ti+Qi−1(20)
聯立(19)、(20)可得
T i − C i P i T i = B i T i + 1 + D i + C i Q i − 1 (21) T_{i}-C_{i}P_{i}T_{i}=B_{i}T_{i+1}+D_{i}+C_{i}Q_{i-1}\tag{21} Ti−CiPiTi=BiTi+1+Di+CiQi−1(21)
整理得
T i = B i A i − C i P i − 1 T i + 1 + D i + C i Q i − 1 A i − C i P i − 1 (22) T_{i}=\frac{B_{i}}{A_{i}-C_{i}P_{i-1}}T_{i+1}+\frac{D_{i}+C_{i}Q_{i-1}}{A_{i}-C_{i}P_{i-1}}\tag{22} Ti=Ai−CiPi−1BiTi+1+Ai−CiPi−1Di+CiQi−1(22)
觀察(20)、(22)式，可以發現
P i = B i A i − C i P i − 1 ; Q i = D i + C i Q i − 1 A i − C i P i − 1 (23) P_{i}=\frac{B_{i}}{A_{i}-C_{i}P_{i-1}}; Q_{i}=\frac{D_{i}+C_{i}Q_{i-1}}{A_{i}-C_{i}P_{i-1}}\tag{23} Pi=Ai−CiPi−1Bi;Qi=Ai−CiPi−1Di+CiQi−1(23)
結合代數形式(16)及矩陣形式(18)、(19)與(23)可以得到有 N N N個格點的各個系數向量
A = ( 4 4 4 4 4 4 4 4 4 ⋯ 4 ) T B = ( 2 − α Δ x 2 − α Δ x 2 − α Δ x ⋯ 0 ) T C = ( 0 2 + α Δ x 2 − α Δ x 2 − α Δ x ⋯ ) T D = ( ( 2 + Δ x ) ϕ 0 0 0 ⋯ ( 2 − Δ x ) ϕ N + 1 ) T \begin{aligned} A & = & (4~~~4~~~4~~~4~~~4~~~4~~~4~~~4~~~4~~~~~\cdots~~~~4)^{T} \\ B & = & (2-\alpha \Delta x~~2-\alpha \Delta x~~2-\alpha \Delta x~~\cdots~~0)^{T} \\ C & = & (0~~2+\alpha \Delta x~~2-\alpha \Delta x~~2-\alpha \Delta x~~\cdots)^{T} \\ D & = & ((2+\Delta x)\phi_{0}~~0~~0~~\cdots~~(2-\Delta x)\phi_{N+1})^{T} \end{aligned} ABCD====(4 4 4 4 4 4 4 4 4 ⋯ 4)T(2−αΔx 2−αΔx 2−αΔx ⋯ 0)T(0 2+αΔx 2−αΔx 2−αΔx ⋯)T((2+Δx)ϕ0 0 0 ⋯ (2−Δx)ϕN+1)T
解得方程在 α = ( − 5 , − 1 , 0 , 1 , 5 ) \alpha=(-5,-1,0,1,5) α=(−5,−1,0,1,5)時的曲線
5. 網格獨立化考核
網格數量對數值計算結果有着重要的影響。當網格足夠細密以至于再進一步加密網格已對數值計算結果基本上沒有影響時所得到的數值解稱為網格獨立解。
為了避免浪費計算資源，在精度足夠高的同時，減少網格數，這裡使用TDMA方法計算在不同網格劃分情況下 x x x=0.7處的值，進行網格獨立化考核。
顯然，當格點數大于80時，再增加格點數時 ϕ \phi ϕ(0.7)的值其變化已經不明顯了，故我們認為格點數達到80時，此時的解為網格獨立解。