标签：frac Delta 数值 rho delta partial 方法 lambda

1-1 数值方法B

非稳态扩散方程

主体公式

\[\frac{\partial (\rho c_p T)}{\partial t}=\frac{\partial}{\partial x}\Big(\lambda\frac{\partial T}{\partial x}\Big)+S(x,t,T) \]

该式左边为时间导数，表示某一点的能量随时间的变化；右边为空间变化，描述热在空间扩散的过程。

其边界条件为：\(T(x,~t=0)\)

其解决思路为：首先从t=0出发，准备初始条件（即上述边界条件）；

当\(t=1\Delta t\)时，对方程离散化，获取线性方程，并利用上述边界条件，得到的解为\(T(x,~t=1\Delta t\)；

对\(t=2\Delta t\)，新的初始条件为\(T(x,~t=1\Delta t)\)，同样进行离散化，并用新边界条件求解。

重复直到\(t=t_{final}\)，此时求出的解为\(T(x,~t=t_{final})\)。

对每一步求解，我们都需要求解一个新的线性系统[A][T]=[B]。

离散化

​​

考虑无源项的公式：

\[\frac{\partial(\rho c_pT)}{\partial t}=\frac{\partial}{\partial x}\Big(\lambda\frac{\partial T}{\partial x}\Big) \]

进行时间和空间上的积分：

\[\rho c_p\int_t^{t+\Delta t}\int_V\frac{\partial T}{\partial t}dVdt =\int_t^{t+\Delta t}\int_V\frac{\partial}{\partial x}(\lambda\frac{\partial T}{\partial x})dVdt \]

\[\rho c_p\int_w^e\left[\int_t^{t+\Delta t}\frac{\partial T}{\partial t}dt\right]Adx =\int_t^{t+\Delta t} \int_w^e \frac{\partial}{\partial x }(\lambda \frac{\partial T}{\partial x})Adxdt \]

\[\rho c_p(T_P^1-T_P^0)A\Delta x =\int_t^{t+\Delta t}\left[\lambda_e\frac{T_E-T_P}{\delta x_e}-\lambda_w\frac{T_P-T_W}{\delta x_w}\right]Adt \]

这里，\(T_P^0\)是t时刻点P处的温度，\(T_P^1\)是\(t+\Delta t\)时刻点P处的温度。

但这一步还不够，为得到温度随时间的变化，设温度\(T(t)\)在时刻t和\(t+\Delta t\)之间呈线性变化。

有如下公式：

\[\int_t^{t+\Delta t}T_Pdt=[fT_P^1+(1-f)T_P^0]\Delta t \]

其中，该式等号左边部分是在时间区间\([t,~\Delta t]\)内，对点P处温度进行积分，得到的是温度的累计效果；

右边部分是用点P处温度在时刻t和\(t+\Delta t\)的加权平均值，乘以时间步长\(\Delta t\)得到的近似值。

在右边部分中，f是时间离散化的权重系数：

• f=0：显式欧拉法/前向差分法，只使用当前时刻温度；
• f=1：隐式欧拉法/后向差分法，只使用下一时刻温度；
• f=0.5：中点法，对当前时刻和下一时刻的温度取平均值。

则一维非稳态热传导公式变为：

\[\rho c_{p}(T_{P}^{1}-T_{P}^{0})\Delta x=\left\{f\left[\lambda_{e}\frac{T_{E}^{1}-T_{P}^{1}}{\delta x_{e}}-\lambda_{w}\frac{T_{P}^{1}-T_{W}^{1}}{\delta x_{w}}\right]+(1-f)\left[\lambda_{e}\frac{T_{E}^{0}-T_{P}^{0}}{\delta x_{e}}-\lambda_{w}\frac{T_{P}^{0}-T_{W}^{0}}{\delta x_{w}}\right]\right\}\Delta t \]

其中：

• 项\(\rho\)​\(c_p\)​\((T_{P}^1\)​\(-T_P^0)\)​\(\Delta x\)是储能项，表示控制节点P在时间段\(\Delta t\)内因为温度变化产生的能量变化，\(\rho c_p \Delta x\)表示节点P控制体积的总热容量；

• 项

  \[f[\lambda_{e}\frac{T_{E}^{1}-T_{P}^{1}}{\delta x_{e}}-\lambda_{w}\frac{T_{P}^{1}-T_{W}^{1}}{\delta x_{w}}] \]

  表示基于未来时刻\(t+\Delta t\)的传导项。
  其物理意义为：-热通量=热导率*温度梯度（即\(\lambda \frac{\Delta T}{\delta x}\)），表示热量从高温区流向低温区的速率。
  其中，\(\lambda_{e}\frac{T_{E}^{1}-T_{P}^{1}}{\delta x_{e}}\)表示热量从节点E传导向节点P，\(-\lambda_{w}\frac{T_{P}^{1}-T_{W}^{1}}{\delta x_{w}}\)表示热量从节点P传导向节点W，这里的负号表示节点P损失热量。
  ​\(\delta x_e\)和\(\delta x_w\)是节点距离，它们决定了温度梯度的大小。

重新排列该方程，最终可得：

\[-f\frac{\lambda_w}{\delta x_w}T_W+[\rho c_p\frac{\Delta x}{\Delta t}+f(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e})]T_P-f\frac{\lambda_e}{\delta x_e}T_E=\rho c_p\frac{\Delta x}{\Delta t}T_P^0+(1-f)[\frac{\lambda_w}{\delta x_w}T_W^0+\frac{\lambda_e}{\delta x_e}T_E^0-(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e})T_P^0] \]

即下述线性方程：

\[a_W\boldsymbol{T}_W+a_P\boldsymbol{T}_P+a_E\boldsymbol{T}_E=b \]

其中，\(a_w=-f\frac{\lambda_w}{\delta x_w}\)，\(a_p=\rho c_p\frac{\Delta x}{\Delta t}+f(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e})\)，\(a_e=-f\frac{\lambda_e}{\delta x_e}\)；

\(b=\rho c_p\frac{\Delta x}{\Delta t}T_P^0+(1-f)[\frac{\lambda_w}{\delta x_w}T_W^0+\frac{\lambda_e}{\delta x_e}T_E^0-(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e})T_P^0]\)。

下面依次讨论三种不同的权重系数f下，方程的形式。

1. f=0：

   \[\left(\rho c_p\frac{\Delta x}{\Delta t}\right)T_P=\frac{\lambda_w}{\delta x_w}T_W^0+\frac{\lambda_e}{\delta x_e}T_E^0+[\rho c_p\frac{\Delta x}{\Delta t}-\left(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e}\right)]T_P^0 \]

   重新排列上述公式：

   \[T_P+\left(\frac{\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e}}{\rho c_p\frac{\Delta x}{\Delta t}}-1\right)T_P^0=\frac{\frac{\lambda_w}{\delta x_w}T_W^0+\frac{\lambda_e}{\delta x_e}T_E^0}{\rho c_p\frac{\Delta x}{\Delta t}} \]

   对照\(a_W T_W+a_P T_P+a_E T_E=b\)，其系数如下：
   ​\(a_P=\rho c_P \frac{\Delta x}{\Delta t}\)，\(a_E=a_W=0\)，
   ​\(b=\frac{\lambda_w}{\delta x_w}T_W^0+\frac{\lambda_e}{\delta x_e}T_E^0-(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e}-\rho c_p\frac{\Delta x}{\Delta t})T_P^0\)
   为确保解的稳定性，当前节点当前时间的温度\(T_P^0\)的系数：
   ​\(\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e}-\rho c_p\frac{\Delta x}{\Delta t}\)需要为负数，确保温度在每个时间步内无不合理增长。
   该系数的意义为：如果系数为正，中心点的温度\(T_P\)会不受控地增长，使得解发散。

   可得：\(\Delta t<\frac{\rho c_p(\Delta x)^2}{2\lambda}\)时解稳定，此时\(\lambda_{w}=\lambda_{e},~\delta x_{w}=\delta x_{e}= \Delta x\)。
   该方法为一阶精度，误差随着\(\Delta t\)线性增长；解耦，每个节点温度可以独立求解，但需严格控制\(\Delta t\)步长。

2. f=0.5：

   \[\begin{aligned}&a_{P}=\rho c_{p} \frac{\Delta x}{\Delta t}+\frac{1}{2} ( \frac{\lambda_{w}}{\delta x_{w}}+\frac{\lambda_{e}}{\delta x_{e}}),\\&a_{E}=\frac{1}{2}\frac{\lambda_{e}}{\delta x_{e}},\\&a_{W}=-\frac{1}{2}\frac{\lambda_{w}}{\delta x_{w}},\\&b=\frac{1}{2}[ \frac{\lambda_{w}}{\delta x_{w}}T_{W}^{0}+\frac{\lambda_{e}}{\delta x_{e}}T_{E}^{0}-\left(\frac{\lambda_{w}}{\delta x_{w}}+\frac{\lambda_{e}}{\delta x_{e}}-2\rho c_{p}\frac{\Delta x}{\Delta t}\right)T_{P}^{0}]\end{aligned} \]

   解稳定的条件同f=0时。
   该方法为二阶精度，\(\Delta t\)小的时候该方法最精确；非解耦，需要线性系统的解。

3. f=1：

   \[\begin{aligned} &a_{P}=\rho c_{p}\frac{\Delta x}{\Delta t}+\frac{\lambda_{w}}{\delta x_{w}}+\frac{\lambda_{e}}{\delta x_{e}}, \\ &a_{E}=-\frac{\lambda_{e}}{\delta x_{e}}, \\ &a_{W}=-\frac{\lambda_{W}}{\delta x_{W}}, \\ &b=\rho c_p\frac{\Delta x}{\Delta t}T_P^0 \end{aligned} \]

   解稳定的条件：\(-\rho c_p\frac{\Delta x}{\Delta t}<0\)，对任何\(\Delta t\)都成立。
   该方法为一阶精度，需要线性系统的解。

1-1.2 数值方法B 续

‍

‍

标签：frac,Delta,数值,rho,delta,partial,方法,lambda
From： https://www.cnblogs.com/yukibvvd/p/18513271/11-numerical-method-b-zkkftc