为方便计算，对区域进行离散化处理，采用隐式差分格式进行离散计算。隐式差分格式如图：
每层材料内部
对第 j j j层材料:
其中， λ j \lambda_j λj表示第 j j j层的热扩散率， c j c_j cj表示第 j j j层的比热容， ρ j \rho_j ρj表示第 j j j层的材料密度。
∂ T j ∂ t = λ j c j ρ j ∂ 2 T j ∂ x 2 \dfrac{\partial T_j}{\partial t} =\dfrac{\lambda_j}{c_j \rho_j}\dfrac{\partial^2 T_j}{\partial x^2} ∂t∂Tj=cjρjλj∂x2∂2Tj
对方程左端采用向后差分格式，方程右端采用中心差分格式，有：
T j ( x i , t k + 1 ) − T j ( x i , t k ) Δ t = λ j c j ρ j T j ( x i + 1 , t k + 1 ) − 2 T j ( x i , t k + 1 ) + T j ( x i − 1 , t k + 1 ) ( Δ x ) 2 \dfrac{T_j(x_i,t_{k+1})-T_j(x_i,t_{k})}{\Delta t} =\dfrac{\lambda_j}{c_j \rho_j} \dfrac{T_j(x_{i+1},t_{k+1})-2T_j(x_i,t_{k+1}) + T_j(x_{i-1},t_{k+1})}{(\Delta x)^2} ΔtTj(xi,tk+1)−Tj(xi,tk)=cjρjλj(Δx)2Tj(xi+1,tk+1)−2Tj(xi,tk+1)+Tj(xi−1,tk+1)
整理得:
− T j ( x i , t k ) = λ j Δ t c j ρ j ( Δ x ) 2 F ( X ) − T j ( x i , t k + 1 ) -T_j(x_i,t_{k}) = \dfrac{\lambda_j \Delta t}{c_j \rho_j (\Delta x)^2} F(X) - T_j(x_i,t_{k+1}) −Tj(xi,tk)=cjρj(Δx)2λjΔtF(X)−Tj(xi,tk+1)
F ( X ) = T j ( x i + 1 , t k + 1 ) − 2 T j ( x i , t k + 1 ) + T j ( x i − 1 , t k + 1 ) F(X) = T_j(x_{i+1},t_{k+1})-2T_j(x_i,t_{k+1}) + T_j(x_{i-1},t_{k+1}) F(X)=Tj(xi+1,tk+1)−2Tj(xi,tk+1)+Tj(xi−1,tk+1)
令 F ( X ) = r j F(X) = r_j F(X)=rj
整理得:
T j ( x i , t k ) = − r j T j ( x i + 1 , t k + 1 ) − r j T j ( x i − 1 , t k + 1 ) + ( 1 + 2 r j ) T j ( x i , t k + 1 ) T_j(x_i,t_k) = -r_jT_j(x_{i+1},t_{k+1})-r_jT_j(x_{i-1},t_{k+1})+(1+2r_j)T_j(x_i,t_{k+1}) Tj(xi,tk)=−rjTj(xi+1,tk+1)−rjTj(xi−1,tk+1)+(1+2rj)Tj(xi,tk+1)
左边界条件
T j ( x , 0 ) = 37 ° T_j(x,0) = 37\degree Tj(x,0)=37°
此前求得的左边界条件：
− k 1 ∂ T 1 ∂ x + k o u t T 1 ∣ x = 0 = k o u t T o u t -k_1 \dfrac{\partial T_1}{\partial x} + k_{out}T_1|_{x=0} = k_{out}T_{out} −k1∂x∂T1+koutT1∣x=0=koutTout
对其进行离散化处理：
− k 1 T ( 2 , k ) − T ( 0 , k ) 2 Δ x = k o u t ( T o u t − T ( 1 , k ) ) -k_1 \dfrac{T_(2,k) - T(0,k)}{2\Delta x} = k_{out}(T_{out} - T(1,k)) −k12ΔxT(2,k)−T(0,k)=kout(Tout−T(1,k))
整理可得：
T ( 0 , k ) = 2 Δ x k o u t T o u t k 1 − 2 Δ x k o u t T ( 1 , k ) k 1 + T ( 2 , k ) T(0,k) = \dfrac{2\Delta xk_{out}T_{out}}{k_1} - \dfrac{2\Delta xk_{out}T(1,k)}{k_1} + T(2,k) T(0,k)=k12ΔxkoutTout−k12ΔxkoutT(1,k)+T(2,k)
根据每层材料内部，有
T j ( 1 , k − 1 ) = − r 1 T j ( 2 , k ) − r 1 T j ( 0 , k ) + ( 1 + 2 r 1 ) T j ( 1 , k ) T_j(1,k-1) = -r_1T_j(2,k)-r_1T_j(0,k)+(1+2r_1)T_j(1,k) Tj(1,k−1)=−r1Tj(2,k)−r1Tj(0,k)+(1+2r1)Tj(1,k)
联立得:
T j ( 1 , k − 1 ) = − r j T j ( 2 , k ) − r j ( 2 Δ x k o u t T o u t k 1 − 2 Δ x k o u t T ( 1 , k ) k 1 + T ( 2 , k ) ) + ( 1 + 2 r j ) T j ( 1 , k ) T_j(1,k-1) = -r_jT_j(2,k)-r_j(\dfrac{2\Delta xk_{out}T_{out}}{k_1} - \dfrac{2\Delta xk_{out}T(1,k)}{k_1} + T(2,k))+(1+2r_j)T_j(1,k) Tj(1,k−1)=−rjTj(2,k)−rj(k12ΔxkoutTout−k12ΔxkoutT(1,k)+T(2,k))+(1+2rj)Tj(1,k)
T j ( 1 , k − 1 ) + r 1 2 Δ x k o u t T o u t k 1 = − 2 r 1 T j ( 2 , k ) + ( 1 + 2 r 1 + r 1 2 Δ x k o u t k 1 ) T j ( 1 , k ) T_j(1,k-1)+r_1\dfrac{2\Delta xk_{out}T_{out}}{k_1} = -2r_1T_j(2,k)+(1+2r_1+r_1\dfrac{2\Delta xk_{out}}{k_1})T_j(1,k) Tj(1,k−1)+r1k12ΔxkoutTout=−2r1Tj(2,k)+(1+2r1+r1k12Δxkout)Tj(1,k)
右边界条件
此前求得的右边界条件：
k 4 ∂ T 1 ∂ x + k o u t T 4 ∣ x = x s k i n = k s k i n T s k i n k_4 \dfrac{\partial T_1}{\partial x} + k_{out}T_4|_{x=x_{skin}} = k_{skin}T_{skin} k4∂x∂T1+koutT4∣x=xskin=kskinTskin
对其进行离散化处理：
k 4 T ( N + 1 , k ) − T ( N − 1 , k ) 2 Δ x = k s k i n ( T s k i n − T ( N , k ) ) k_4 \dfrac{T_(N+1,k) - T(N-1,k)}{2\Delta x} = k_{skin}(T_{skin} - T(N,k)) k42ΔxT(N+1,k)−T(N−1,k)=kskin(Tskin−T(N,k))
整理可得：
T ( N − 1 , k ) = T ( N + 1 , k ) − 2 Δ x k s k i n T s k i n k 4 + 2 Δ x k s k i n T ( N , k ) k 4 T(N-1,k) = T(N+1,k) -\dfrac{2\Delta xk_{skin}T_{skin}}{k_4} + \dfrac{2\Delta xk_{skin}T(N,k)}{k_4} T(N−1,k)=T(N+1,k)−k42ΔxkskinTskin+k42ΔxkskinT(N,k)
根据每层材料内部，有
T j ( N , k − 1 ) = − r 1 T j ( N + 1 , k ) − r 1 T j ( N − 1 , k ) + ( 1 + 2 r 1 ) T j ( N , k ) T_j(N,k-1) = -r_1T_j(N+1,k)-r_1T_j(N-1,k)+(1+2r_1)T_j(N,k) Tj(N,k−1)=−r1Tj(N+1,k)−r1Tj(N−1,k)+(1+2r1)Tj(N,k)
联立得:
T j ( N , k − 1 ) = − r j T j ( N + 1 , k ) − r j ( T ( N + 1 , k ) − 2 Δ x k s k i n T s k i n k 4 + 2 Δ x k s k i n T ( N , k ) k 4 ) + ( 1 + 2 r ) T j ( N , k ) T_j(N,k-1) = -r_jT_j(N+1,k)-r_j(T(N+1,k) -\dfrac{2\Delta xk_{skin}T_{skin}}{k_4} + \dfrac{2\Delta xk_{skin}T(N,k)}{k_4})+(1+2r)T_j(N,k) Tj(N,k−1)=−rjTj(N+1,k)−rj(T(N+1,k)−k42ΔxkskinTskin+k42ΔxkskinT(N,k))+(1+2r)Tj(N,k)
T j ( N , k − 1 ) − 2 Δ x k s k i n T s k i n k 4 = − 2 r j T j ( N + 1 , k ) + ( 1 + 2 r j − r j 2 Δ x k s k i n k 4 ) T j ( N , k ) T_j(N,k-1)-\dfrac{2\Delta xk_{skin}T_{skin}}{k_4} = -2r_jT_j(N+1,k)+(1+2r_j-r_j\dfrac{2\Delta xk_{skin}}{k_4})T_j(N,k) Tj(N,k−1)−k42ΔxkskinTskin=−2rjTj(N+1,k)+(1+2rj−rjk42Δxkskin)Tj(N,k)
N = L N=L N=L
本文参考自【数学建模】2018年国赛A题详解（一）
更多建模相关信息看Github仓库(包括真题练习)