非稳态导热微分方程 + 无限大平板的非稳态传热问题

标签：partial varTheta dfrac mu 无限大 frac 稳态 导热 mathrm

外层肉和内层肉的球坐标传热微分方程组

\[\begin{cases} \dfrac1a\dfrac{\partial T_1}{\partial t}=\dfrac1{r^2}\dfrac\partial{\partial r}\left(r^2\dfrac{\partial T_1}{\partial r}\right) \\[2ex] \dfrac1a\dfrac{\partial T_3}{\partial t}=\dfrac1{r^2}\dfrac\partial{\partial r}\left(r^2\dfrac{\partial T_3}{\partial r}\right) \end{cases} \]

初值条件

\[t=0:0\le r\le R_3,T_1=T_2=T_3=T_0 \]

边界条件

\[\begin{matrix} \displaystyle t>0,r=0,\frac{\partial T_1}{\partial r}=0 \\[2ex] \displaystyle t>0:r=R_1,h_\infty(T_1-T_2)=-k\frac{\partial T_1}{\partial r} \\[2ex] \displaystyle t>0:r=R_2,h(T_2-T_3)=-k\frac{\partial T_3}{\partial r} \\[2ex] \displaystyle t>0:r=R_3,h(T_3-T_\infty)=-k\frac{\partial T_3}{\partial r} \end{matrix} \]

为了解这类方程，引入过余温度

\[\varTheta=T-T_\infty \]

则外层的方程和条件化为

\[\begin{matrix} \displaystyle\frac1a\frac{\partial\varTheta}{\partial t}=\frac1{r^2}\frac\partial{\partial r}\left(r^2\frac{\partial\varTheta}{\partial r}\right)=\frac{\partial^2\varTheta}{\partial r^2}+\frac2r\frac{\partial\varTheta}{\partial r} \\[2ex] \displaystyle t>0,r=R_3,h\varTheta_3+k\frac{\partial\varTheta_3}{\partial r}=0 \end{matrix} \]

再巧妙地设定一个新变量 \(\varPhi=\varTheta r\)，则有 \(\dfrac{\partial\varPhi}{\partial t}=r\dfrac{\partial\varTheta}{\partial t},\dfrac{\partial\varPhi}{\partial r}=r\dfrac{\partial\varTheta}{\partial r}+\varTheta,\dfrac{\partial^2\varPhi}{\partial r^2}=2\dfrac{\partial\varTheta}{\partial r}+\dfrac{\partial^2\varTheta}{\partial r^2}\)。原方程化为简洁的形式

\[\frac{\partial\varPhi}{\partial t}=a\frac{\partial^2\varPhi}{\partial r^2} \]

分离变量法，寻求方程有如下形式的解 \(\varPhi(r,t)=X(r)\tau(t)\)，从而化为 \(\dfrac1{a\tau}\dfrac{\mathrm d\tau}{\mathrm dt}=\dfrac1X\dfrac{\mathrm d^2X}{\mathrm dr^2}\)。此式左端是 \(t\) 的函数，右端是 \(x\) 的函数。为等式成立，两边只能等于常数。要是这个常数为 \(0\) 也太扯了，正数稍微算一下也不对，肯定是负数，取 \(-a^2\)。则 \(\dfrac1{a\tau}\dfrac{\mathrm d\tau}{\mathrm dt}=\dfrac1X\dfrac{\mathrm d^2X}{\mathrm dr^2}=-a^2\)，有

\[\begin{aligned} \displaystyle\tau&=C_1\mathrm e^{-a^2at} \\ \displaystyle X&=C_2\sin ax+C_3\cos ax \\ \displaystyle\varPhi&=X\tau=e^{-a^2at}(A\sin ax+B\cos ax) \end{aligned} \]

由径向传热的对称性，\(r=0\) 时 \(\varPhi=0\)，易得 \(B=0\)，进而 \(\varPhi=A\mathrm e^{-a^2at}\sin ax\)。但其实这一步是外推的思想，我不是特别确定，下面算错了不怪我。

此时的外边界条件已经化为 \(t>0,r=R_3,\left(h-\dfrac k{R^3}\right)\varPhi+k\dfrac{\partial\varPhi}{\partial r}=0\)，代入有 \(\left(h-\dfrac k{R^3}\right)\sin ar+ka\cos ar=0\)，即 \(\tan aR^3=\dfrac{ka}{k/R_3-h}=\dfrac{aR_3}{1-\textit{Bi}}\)，式中 \(\textit{Bi}=hR_3/k\) 为毕渥数。上式特征方程即可确定满足边界条件的所有 \(a\) 特征值 \(a_i\)，解为它们的线性组合。若用 \(\mu_n=a_nR_3\) 代入上式，得

\[\varPhi(r,t)=\sum_{n=1}^\infty A_n\mathrm e^{-a_n^2at}\sin\frac{\mu_n}{R_3}r \]

由初始条件，\(t=0\) 时，\(\varPhi(r,t)=r\varTheta_0=r(T_0-T_\infty)=\sum_{n=1}^\infty A_n\sin\dfrac{\mu_n}{R_3}r\)。一通积分变换猛如虎！

\[(T_0-T_\infty)\int_0^{R_3}r\sin\frac{\mu_mr}{R_3}\mathrm dr=\sum_{n=1}^\infty A_n\int_0^{R_3}\sin\frac{\mu_mr}{R_3}\frac{\mu_nr}{R_3}\mathrm dr=\frac{R_3(\mu_n\sin\mu_m\cos\mu_n-\mu_m\sin\mu_n\cos\mu_m)}{\mu_m^2-\mu_n^2} \]

系数 \(A\) 由以下方程给出

\[A_n=\cfrac{(T_0-T_\infty)\int_0^{R_3}r\sin\cfrac{\mu_nr}{R_3}\mathrm dr}{\int_0^{R_3}\sin\cfrac{\mu_nr}{R_3}\sin\cfrac{\mu_nr}{R_3}\mathrm dr}=\frac{2R_3(T_0-T_\infty)(\sin\mu_n-\mu_n\cos\mu_n)}{\mu_n(\mu_n-\sin\mu_n\cos\mu_n)} \]

一通代入求出来了外层的分布，其实我也不知道对不对，但我估计弹幕里也不会有人闲得看到这吧哈哈哈

\[T=T_\infty+\sum_{n=1}^\infty\frac{2(T_0-T_\infty)(\sin\mu_n-\mu_n\cos\mu_n)}{(\mu_n-\sin\mu_n\cos\mu_n)}\cfrac{\sin\cfrac{\mu_nr}{R_3}}{\mu_n\cfrac r{R_3}}\mathrm e^{-a^2at} \]

一口老血吐出来，内层还要再依赖于复杂对流边界条件算进去，我决定放弃，改用数值解法。下面的图都是用数值解法求出来的。

---

将毛肚看作厚度为 \(2L\) 的无限大平板。初始均匀温度 \(T_0\)。在 \(t=0\) 时刻，平板两个侧面与温度 \(T_\infty\) 的流体接触，流体与平板间的对流换热系数为 \(h\)。易得导热微分方程

\[\frac{\partial T}{\partial t}=a\frac{\partial^2T}{\partial x^2} \]

初始条件

\[t=0:0\le x\le L,T=T_0 \]

边界条件

\[\begin{matrix} \displaystyle t>0:x=0,\frac{\partial T}{\partial x}=0 \\[2ex] \displaystyle t>0:x=L,h(T-T_\infty)=-k\frac{\partial T}{\partial x} \end{matrix} \]

引入无量纲过余温度 \(\varTheta=\dfrac{T-T_\infty}{T_0-T_\infty}\)，将上述条件齐次化

\[\begin{matrix} \displaystyle \frac{\partial\varTheta}{\partial t}=a\frac{\partial^2\varTheta}{\partial x^2} \\[2ex] \displaystyle t=0:0\le x\le L,\varTheta=1 \\[2ex] \displaystyle t>0:x=0,\frac{\partial\varTheta}{\partial x}=0 \\[2ex] \displaystyle t>0:x=L,h\varTheta+k\frac{\partial\varTheta}{\partial x}=0 \end{matrix} \]

分离变量法，寻求方程有如下形式的解 \(\varTheta(x,t)=X(x)\tau(t)\)，从而化为 \(\dfrac1{a\tau}\dfrac{\mathrm d\tau}{\mathrm dt}=\dfrac1X\dfrac{\mathrm d^2X}{\mathrm dx^2}\)。此式左端是 \(t\) 的函数，右端是 \(x\) 的函数。为等式成立，两边只能等于常数。要是这个常数为 \(0\) 也太扯了，正数稍微算一下也不对，肯定是负数，取 \(-a^2\)。则 \(\dfrac1{a\tau}\dfrac{\mathrm d\tau}{\mathrm dt}=\dfrac1X\dfrac{\mathrm d^2X}{\mathrm dx^2}=-a^2\)，有

\[\begin{aligned} \tau&=C_1\mathrm e^{-a^2at} \\ X&=C_2\sin ax+C_3\cos ax \end{aligned} \]

特解为 \(\varTheta=\tau X=C_1\mathrm e^{-a^2at}(C_2\sin ax+C_3\cos ax)=\mathrm e^{-a^2at}(B\sin ax+A\cos ax)\)，此时有

\[\frac{\partial\varTheta}{\partial x}=\mathrm e^{-a^2at}(Ba\cos ax-Aa\sin ax) \]

当 \(x=0\) 时，\(\dfrac{\partial\varTheta}{\partial x}=\mathrm e^{-a^2at}Ba\)。根据零点的边界条件，\(\dfrac{\partial\varTheta}{\partial x}=0\)，显然有 \(B=0\)，进而

\[\varTheta=\mathrm e^{-a^2at}A\cos ax \]

再利用对流边界条件 \(h\mathrm e^{-a^2at}A\cos aL+k\mathrm e^{-a^2at}(-Aa\sin aL)=0\)，得 \(\cot(aL)=\dfrac{ka}h=\dfrac{aL}{\textit{Bi}}\)，式中 \(\textit{Bi}\) 为毕渥数。上式特征方程即可确定满足边界条件的所有 \(a\) 特征值 \(a_i\)，解为它们的线性组合

\[\varTheta(x,t)=\sum_{n=1}^\infty A_n\mathrm e^{-a_n^2at}\cos a_nx \]

由初始条件 \(t=0\) 时，\(\varTheta(x,t)=\sum_{n=1}^\infty A_n\cos a_nx=1\)，显然 \(\sum_{n=1}^\infty A_n\cos a_nx\) 是 \(1\) 对于 \(\cos a_nx\) 展开的傅里叶级数，因此

\[A_n=\frac{\int_0^L\cos a_nx\mathrm dx}{\int_0^L\cos^2a_nx\mathrm dx}=\frac{2\sin a_nL}{a_nL+\sin a_nL\cos a_nL}=\cfrac{\cfrac1{a_n}\sin a_nL}{\cfrac L2+\cfrac1{4a_n}\sin2a_nL} \]

至此得到了方程的解的形式为

\[\varTheta=\frac{T-T_\infty}{T_0-T_\infty}=2\sum_{n=1}^\infty\frac{\sin a_nL}{a_nL+\sin a_nL\cos a_nL}\mathrm e^{-a_n^2at\cos a_nx} \]

如果毛肚和水换热很快，那问题倒是简单了，特征方程化为 \(\cot(aL)=\dfrac{aL}{\textit{Bi}}=0\)，相应地，解为 \(a_nL=\dfrac{2n-1}{2}\pi\)。代入可得无限大平板的温度分布为

\[\varTheta=\sum_{n=1}^\infty\frac{4(-1)^{n+1}}{(2n-1)\pi}\cos\hspace{-0.25em}\left(\frac{2n-1}2\frac{\pi x}L\right)\mathrm e^{-\left(\frac{2n-1}2\pi\right)^2\frac{at}{L^2}} \]

但我们不会那么赖皮地用第一类边界条件，还是老老实实按第三类边界条件算吧。代入毛肚的参数 \(c_p=3204~\text{J/kg}\cdot\text K\)，密度 \(\rho=1170~\text{kg/m}^3\)，导热系数 \(k=0.5~\text{W/m}\cdot\text K\)，故 \(a=\dfrac{k}{c_p\rho}=1.3\times10^{-7}~\text{m}^2\text{/s}\)。设毛肚厚度 \(2L=1.8~\text{mm}\)，水沸腾对流换热系数取 \(h=4000~\text{W/m}^2\text K\)。代入求解并画出表面和内部升温的曲线。

标签：partial,varTheta,dfrac,mu,无限大,frac,稳态,导热,mathrm
From： https://www.cnblogs.com/laoshan-plus/p/18365003