在斯托克斯振荡流下具有时间分数阶加速度的牛顿流体自由对流稳定性数值模拟

Numerical Simulation of Free Convection Stability of Newtonian Fluids with Time-Fractional Acceleration under Stokes Oscillatory Flow

摘要: 为揭示时间分数阶动力学对牛顿流体自由对流行为的调控机制，本研究将时间分数阶导数纳入流体自由对流控制方程，结合相平面图分析，通过数值模拟与Mathematica参数化算法(每对控制参数生成10⁴个样本点)，系统探究分数阶特性对流体稳定性及关键物理量的影响。结果表明：时间分数阶导数的非局部特性与记忆效应可显著增强系统不稳定性，催生以不稳定振荡为特征的复杂对流模式；牛顿流体的壁面剪切应力与努塞尔数随分数阶参数($\alpha$ )减小而降低，随幂律指数增大而升高；所建算法与分析方法实现了稳定/不稳定区域的精准界定。本研究完善了复杂流体系统的理论框架，深化了记忆效应与反常扩散对流体稳定性及混沌转变的认知，为涉及牛顿流体的工程与工业应用提供了控制优化的理论支撑与新视角。

Abstract: To reveal the regulatory mechanism of time-fractional dynamics on the natural convection behavior of Newtonian fluids, this study incorporates the time-fractional derivative into the governing equations of fluid free convection. Combined with phase plane analysis, numerical simulations and a Mathematica-based parameterized algorithm (generating 10⁴ sample points for each pair of control parameters) are employed to systematically investigate the effects of fractional characteristics on fluid stability and key physical quantities. The results show that the nonlocal property and memory effect of the time-fractional derivative can significantly enhance system instability, leading to complex convection patterns characterized by unstable oscillations; the wall shear stress and Nusselt number of Newtonian fluids decrease with the reduction of the fractional parameter ($\alpha$ ) and increase with the rise of the power-law index; the established algorithm and analytical method realize the accurate delineation of stable/unstable regions. This study improves the theoretical framework of complex fluid systems, deepens the understanding of the impacts of memory effect and anomalous diffusion on fluid stability and chaos transition, and provides theoretical support and new perspectives for the control and optimization in various engineering and industrial applications involving Newtonian fluids.

1. 引言

流体的流动与传热普遍存在“历史依赖”特性，从聚合物熔体粘弹性响应、多孔介质渗流滞后到非牛顿流体长期力学行为，这类时间累积效应相关现象难以通过传统整数阶微积分“局部作用”假设精准刻画。分数阶微积分凭借非局部性与记忆特性，突破整数阶模型局限，成为描述流体时间效应的有力工具。

时间效应是流体动力学核心特性，反映运动状态对前期历程的依赖，贯穿各类传递过程：粘弹性流体中表现为应力松弛与蠕变，多孔介质中体现为压力传导滞后，分数阶流体系统中则通过记忆效应调控动量与热量输运。其研究兼具理论与工程价值：理论上，Kida等人[1]通过纳维–斯托克斯方程发现不可压缩粘性流体从稳定到混沌的演化路径，Badrinath等人[2]亦通过该方程探索了不稳定流体流动的混沌路径；Akyildiz和Bellout [3]在瑞利–贝纳德热对流系统中，分析了温度相关粘度对混沌行为的影响。工程上，Farajzadeh和Tohidi [4]证实无序平流可增强螺旋盘管内幂律流体混合传热，Vajravelu等人[5] [6]则分别探究了热物理性质对幂律流体薄膜流动、热扩散对幂律流体流动的影响。

相关研究实现多维度拓展：研究体系从牛顿流体延伸至粘弹性、纳米、多孔介质等多元流体，Khayat [7] [8]揭示了粘弹性流体热对流的混沌特性及流体弹性对混沌转变的调控作用；研究范畴涵盖多物理场耦合与混沌分析，Widmann等人[9]发现对流回路层流与湍流的不同混沌动力学特征，Kolodner等人[10]明确二元流体对流的色散混沌模式，Albers和Sprott [11]指出高维系统通过Neimark-Sacker分岔序列通往混沌，Labonia和Guj [12]则发现环形区域几何构型影响自然对流向混沌的转变。此外，Boujelbene等人[13]结合傅里叶定律研究了幂律流体热力学，Shahane等人[14]关注不确定性对三维自然对流的影响，Weber和Arfken [15]为非线性方法与混沌研究提供了重要理论支撑。​

然而，当前研究仍存在不足：不同分数阶算子适用性差异不明、多因素耦合调控机制不清、工程定量表征方法待完善。基于此，本文结合分数阶微积分，系统分析时间效应对流动稳定性、剪切应力传递及传热效率的影响，为复杂流体系统研究与应用提供参考。

2. 物理模型与数值计算

2.1. 物理模型

本研究以Stokes第二问题为研究对象，重点研究无限大平板上牛顿流体自由对流的稳定与不稳定性问题(图1)。Stokes第二问题描述了无限大平板突然启动或振荡运动时上方粘性流体的非定常流动响应，该问题不仅揭示了流体动量扩散的瞬态特性，也是研究更复杂流体力学问题的基础，在海洋工程、航空航天工程、生物医学工程等领域具有重要应用价值。在无限大平板假设下，N-S方程简化为线性偏微分方程：

$\frac{\partial u}{\partial t}=v\frac{{\partial }^{2}u}{\partial y^2}$

其中，$v$ 为运动粘度，$u\left(y,t\right)$ 为流向速度。

Figure 1. Model of stokes’ second problem

图1. Stokes第二问题模型

我们考虑一个具有时间分数阶加速度的牛顿流体模型，其在重力作用下流经无限大平板。流场限定在$y>0$ 一个区域内，其中y是在垂直于表面的方向上进行测量。无穷远处的流体具有恒定速度，平板的运动在$t=0$ 时发生变化。在时刻$t=0^{+}$ ，平板在其自身平面内受到两种不同类型的速度$u_w=U\cos \omega t$ 和$U\sin \omega t$ 影响，从而产生诱导流，其中$\omega$ 是平板振荡频率。同时，平板温度升高到$T_w$ ，然后保持在同一水平。

2.2. 控制方程

基于复杂流动的记忆性与反常动力学特征，分数阶导数的积分核结构赋予了方程描述历史效应的能力，其物理本质对应流体运动的亚扩散行为，为研究斯托克斯振荡流场下牛顿流体自由对流的稳定性提供了更精准的数学工具。基于此，自由对流影响下瞬态分数阶流动的动量守恒方程和能量守恒方程可表示为如下形式[16]：

动量守恒方程

$\rho {\chi }^{\alpha -1}\frac{{\partial }^{\alpha }u}{\partial t^{\alpha }}=\mu \frac{{\partial }^{2}u}{\partial y^2}-\sigma B_0^{2}u+\rho {\beta }_{f}g\left(T-T_{\infty }\right)$ (1)

能量守恒方程

$\rho c_p{\chi }^{\alpha -1}\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}=k\frac{{\partial }^{2}T}{\partial y^2}+\mu {\left(\frac{\partial u}{\partial y}\right)}^2$ (2)

式中，$u\left(y,t\right)$ 为流体的速度，$t$ 为时间，$\alpha$ 为分数阶参数，$\chi$ 为时间尺度参数，$\mu$ 为流体的动力粘度，$\rho$ 为流体的密度，$y$ 为竖直方向上的坐标，${\beta }_f$ 为流体的热膨胀系数，$g$ 为重力加速度，$T\left(y,t\right)$ 为流体的温度，$k$ 为流体的导热系数，$c_p$ 为流体的定压比热容。

在动量方程和能量方程(1)~(2)中，涉及分数阶导数的首项与其他项的物理单位不同。为保持量纲一致性，引入时间尺度参数$\chi ={\omega }^{-1}$ ，使得动量方程中的首项系数$\rho$ 变为$\rho {\chi }^{\alpha -1}$ ，能量方程中的首项系数变为$\rho c_p{\chi }^{\alpha -1}$ 。这种调整确保了所有项的单位一致。无量纲化后，该因子会适当消去，从而得到一个统一的无量纲能量方程。

2.3. 初始和边界条件

初始条件：($t=0$ , $y>0$ )

$u\left(y,0\right)=0,T\left(y,0\right)=T_{\infty }$ (3)

边界条件：

在$y=0$ ，$t>0$ ，

$u\left(0,t\right)=u_0\cos \left(\omega t\right)或u_0\sin \left(\omega t\right),T\left(0,t\right)=T_W$ (4)

在$y\to \infty$ ，$t>0$ ，

$u\left(y,t\right)\to 0,T\left(y,t\right)\to T_{\infty }$ (5)

2.4. 无量纲化

上述方程和边界条件通过以下方式进行无量纲化：

$\begin{array}{l}\tau =\omega t,\eta =\sqrt{\frac{\omega \rho }{\mu }}y,U=\frac{u}{u_0},\theta =\frac{T-T_{\infty }}{T_W-T_{\infty }}\\ Gr=\frac{{\beta }_{f}g\Delta T}{\omega u_0},M=\frac{\sigma B_o^2}{{\rho }_f\omega },Ec=\frac{u_o^2}{c_p\Delta T_W},Pr=\frac{v_f}{{\alpha }_f}\end{array}$ (6)

经过式(6)给出无量纲关系，对方程(1)~(5)进行无量纲化，可以得出无量纲控制方程及无量纲定解条件如下：

动量方程和能量方程无量纲形式

$\frac{{\partial }^{\alpha }U}{\partial {\tau }^{\alpha }}=\frac{{\partial }^{2}U}{\partial {\eta }^2}-MU+Gr\theta$ (7)

$\frac{{\partial }^{\alpha }\theta }{\partial {\tau }^{\alpha }}=\frac{1}{Pr}\frac{{\partial }^2\theta }{\partial {\eta }^2}+Ec{\left(\frac{\partial U}{\partial \eta }\right)}^2$ (8)

初始条件和边界条件

$\begin{array}{lll}U\left(\eta ,0\right)=0;\hfill & \theta \left(\eta ,0\right)=0\hfill & \eta \ge 0,\hfill \\ U\left(0,\tau \right)=\cos \left(\tau \right)/\sin \left(\tau \right);\hfill & \theta \left(0,\tau \right)=1\hfill & \tau \ge 1,\hfill \\ U\left(\eta ,\tau \right)\to 0;\hfill & \theta \left(\eta ,\tau \right)\to 0\hfill & \eta \to \infty .\hfill \end{array}$ (9)

以无量纲形式表示的关键量，如壁面处的剪切应力和热传递为：

$\begin{array}{l}{C}_f={\tau }_W\left(t\right)={\frac{\partial U}{\partial \eta }|}_{\eta =0}\\ N{u}_x{\left(R{e}_x\right)}^{-1/2}={-\frac{\partial \theta }{\partial \eta }|}_{\eta =0}\end{array}$ (10)

3. 数值方法

为了在区域$0\le \eta \le {\eta }_{\max }$ ，$0\le \tau \le {\tau }_{\max }$ 且${\tau }_{\max }\le 100\pi$ 上求解问题(7)至(9)，当$\alpha =1$ 时，我们使用线方法的数值方法。依次采用四阶有限差分法和五阶隐式向后差分的方法对方程组(7)至(8)在$0\le \eta \le {\eta }_{\max }$ 和$0\le \tau \le {\tau }_{\max }$ 上进行离散化。在$0<\alpha <1$ 的情况下，我们利用$2-\alpha$ 阶隐式差分格式来逼近Caputo意义下的分数阶导数。由于应用了有限差分法，带有初始条件和边界条件(9)的偏微分方程组(7)至(8)被转化为一个常微分方程组。然后，利用隐式向后差分的方法求解了一个常微分系统。问题的数值解由长度为$m$ 的两个向量$\overrightarrow{W}\left(\tau \right)$ 和$\overrightarrow{\Theta }\left(\tau \right)$ 表示表示，其中$m$ 是网格点${\eta }_1,{\eta }_2,\cdots ,{\eta }_m$ 的数量，且${\eta }_1=0$ 和${\eta }_m={\eta }_{\max }$ 。我们定义微分矩阵$M_1$ 和$M_2$ ，使得函数及其一阶和二阶导数近似如下：

$\begin{array}{l}U\left({\eta }_i,\tau \right)\approx W_i\left(\tau \right),\frac{\partial U}{\partial \eta }\approx M_1\cdot \overrightarrow{W}\left(\tau \right),\frac{{\partial }^{2}U}{\partial {\eta }^2}\approx M_2\cdot \overrightarrow{W}\left(\tau \right)\\ \theta \left({\eta }_i,\tau \right)\approx {\Theta }_i\left(\tau \right),\frac{\partial \theta }{\partial \eta }\approx M_1\cdot \overrightarrow{\Theta }\left(\tau \right),\frac{{\partial }^2\theta }{\partial {\eta }^2}\approx M_2\cdot \overrightarrow{\Theta }\left(\tau \right)\end{array}$ (11)

通过将数值解(11)代入方程(7)至(8)中，在$\alpha =1$ 的情况下我们使用线方法将方程(7)至(8)转变为常微分方程组：

$\begin{array}{l}\frac{\text{d}{W}_i}{\text{d}\tau }=M_2\cdot W_i-M{W}_i+Gr{\Theta }_i,W_i\left(0\right)=\delta {\eta }_i;\\ \frac{\text{d}{\Theta }_i}{\text{d}\tau }=\frac{1}{Pr}{\left(M_2\cdot \Theta \right)}_i+Ec{\left(M_1\cdot W\right)}_i^2,{\Theta }_i\left(0\right)=0;\\ i=2,\cdots ,m-1.\end{array}$ (12)

这里引入一个小的常数$\delta \le {10}^{-12}$ ，以避免在$n<1$ 的情况下的奇异性。为了排除边界上的代数关系，我们对边界条件进行微分，并以以下形式添加初始条件：

$\begin{array}{l}\frac{\text{d}{W}_1}{\text{d}\tau }=-\sin \left(\tau \right),\frac{\text{d}{W}_m}{\text{d}\tau }=0,W_1\left(0\right)=1,W_m\left(0\right)=0,\\ \frac{\text{d}{\Theta }_1}{\text{d}\tau }=0,\frac{\text{d}{\Theta }_m}{\text{d}\tau }=0,{\Theta }_1\left(0\right)=1,{\Theta }_m\left(0\right)=0.\end{array}$ (13)

在$0<\alpha <1$ 的情况下，我们有：

$\begin{array}{l}{\frac{{\partial }^{\alpha }{W}_i}{\partial {\tau }^{\alpha }}|}_{\tau ={\tau }_k}=\frac{\Delta {\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left(W_i^{k+1}-W_i^k+{\displaystyle \sum _{j=1}^k\left(W_i^{k-j+1}-W_i^{k-j}\right)c_j^{\alpha }}\right)\\ {\frac{{\partial }^{\alpha }{\Theta }_i}{\partial {\tau }^{\alpha }}|}_{\tau ={\tau }_k}=\frac{\Delta {\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left({\Theta }_i^{k+1}-{\Theta }_i^k+{\displaystyle \sum _{j=1}^k\left({\Theta }_i^{k-j+1}-{\Theta }_i^{k-j}\right)c_j^{\alpha }}\right)\end{array}$ (14)

其中$\Delta \tau$ 是积分步长，${\tau }_k=k\Delta \tau$ ，$k=0,1,2,\cdots ,N$ ，$c_j^{\alpha }=(j{+}^{1-\alpha }-j^{1-\alpha }$ ，$j=1,2,3,\cdots ,k$ ，$W_i^k=W_i\left({\tau }_k\right)$ ，${\Theta }_i^k={\Theta }_i\left({\tau }_k\right)$ 等等。为了求解常微分方程组，我们使用了在Mathematica 14版中实现的标准求解器NDSolve。在$0<\alpha <1$ 的情况下，我们在每一步都使用方程(7)至(8)的线性近似形式：

$\begin{array}{l}\frac{{\partial }^{\alpha }{W}_i}{\partial {\tau }^{\alpha }}=M_2\cdot W_i^{k+1}-M\cdot W_i^{k+1}+Gr{\Theta }_i^{k+1},\\ \frac{{\partial }^{\alpha }{\Theta }_i}{\partial {\tau }^{\alpha }}=\frac{1}{Pr}{\left(M_2\cdot \Theta \right)}_i^{k+1}+Ec{\left[{\left(M_1\cdot W\right)}^2\right]}_i^{k+1},\\ i=2,\cdots ,m-1;k=0,1,2,\cdots ,N-1.\end{array}$ (15)

在这种情况下，我们直接使用边界和初始条件(9)：

$\begin{array}{l}{W}_1^{k+1}=\cos {\tau }_k和\sin {\tau }_k,W_m^{k+1}=0,k=0,1,\cdots ,N-1,\\ {\Theta }_1^{k+1}=1,{\Theta }_m^{k+1}=0,k=0,1,\cdots ,N-1,\\ W_i^o=0,{\Theta }_i^o=0,i=1,2,\cdots ,m.\end{array}$ (16)

由于模型(7)至(8)具有拟抛物型，因此，对于具有四阶微分矩阵的差分格式，以及在$2-\alpha$ 阶算法(9)的情况下，我们有以下估计[17]：

$\max \left(\left|W_i\left({\tau }_j\right)-U\left({\eta }_i,{\tau }_j\right)\right|,\Vert {\Theta }_i\left({\tau }_j\right)-\theta \left({\eta }_i,{\tau }_j\right)\Vert \right)\le C\left(\Delta {\tau }^a+h^b\right)$ (17)

其中$C,a,b$ 是一些常数。理论上应该是$a=2-\alpha$ ，$b=4$ 。

4. 追踪稳定/不稳定解的算法

因为热对流有稳定和不稳定模式，为此我们需要在参数空间(Ec、Gr、M)中区分这些模式——见图2。由于在边界上有边界条件(9)的偏微分方程组，且(7)至(8)不能用标准线性近似来估计稳定性。为了确定

Figure 2. Stable and unstable regions under different boundary conditions

图2. 不同边界条件下的稳定和不稳定区域

稳定和不稳定解，我们使用了基于Mathematica 14中When Event的数值方法。通过假设稳定解在时间上具有有限范数$\Theta$ ，而不稳定解的范数$\Theta$ 随时间增加。设范数的极限为$\Vert \overrightarrow{\Theta }\left(\tau \right)\Vert \le {\Theta }_{\max }$ 。求解$\Vert \overrightarrow{\Theta }\left(\tau \right)\Vert ={\Theta }_{\max }$ 方程，我们可以定义一个函数$\tau \left(\alpha ,Ec,Gr,M\right)$ ，然后用它来可视化稳定和不稳定区域——见图3。稳定区域定义为具有$\tau \left(\alpha ,Ec,Gr,M\right)\ge t_m$ 的点的集合，不稳定区域是具有$\Vert \overrightarrow{\Theta }\left(\tau \right)\Vert \ge {\Theta }_{\max }$ 的点的集合，以计算得到的点$\tau \left(\alpha ,Ec,Gr,M\right)$ 的集合作为输出。

图2中绘制了稳定(红色)和不稳定(紫色)区域作为格拉晓夫数(Gr)和埃克特数(Ec)的函数。稳定区域表示流体流动保持可预测和规则的条件，而不稳定区域表示流动变得混沌且对初始条件敏感的条件，在稳定区域中导致更规则的行为，在不稳定区域中增加复杂性。这些区域之间的过渡由绿色、蓝色阴影标记，表示不同程度的稳定性。

比较速度是$\text{cos}\left(\tau \right)$ 或$\text{sin}\left(\tau \right)$ 的函数时的稳定和不稳定区域，我们观察到稳定和不稳定区域的总体分布在两种情况下仍然相似。然而，这些区域的具体边界和范围可能由于周期性驱动力的差异而略有变化。分数阶导数始终影响两种速度函数的稳定性，在不稳定区域中增加复杂性和混沌行为，同时在稳定区域中有助于更可预测的动力学。

5. 结果与讨论

在本节中，我们根据数值计算结果生成了相平面图、表面摩擦和努塞尔数分布的图形。本节将通过分析在余弦和正弦边界条件下的相平面图，进一步分析分数阶导数对牛顿流体自由对流混沌行为的影响。诸如表面摩擦和努塞尔数等物理量对于理解具有不同分数阶参数牛顿流体的剪切应力和传热特性至关重要。图形结果的详细解释如下：

5.1. 相平面图

图3中的相平面图描绘了牛顿流体在边界条件$U\left(0,\tau \right)=\text{cos}\left(\tau \right)$ 和$U\left(0,\tau \right)=\text{sin}\left(\tau \right)$ 下速度和温度分布之间的关系，反应了分数阶导数对自由对流稳定性的影响。针对分数阶参数$\alpha$ 的四个不同值($\alpha =1$ 、$\alpha =0.9$ 、$\alpha =0.8$ 和$\alpha =0.7$ )绘制了相平面图，提供了关于不同的$\alpha$ 值如何影响牛顿流体的稳定性和动力学的直观认识。

Figure 3. Effect of fractional-order parameters on the phase portraits of Newtonian fluids

图3. 分数阶参数对牛顿流体相平面图的影响

在边界条件下的牛顿流体的相平面图显示出随着$\alpha$ 的减小从周期性行为向混沌行为的转变。在$\alpha =1$ 时，轨迹更加规则，表现出稳定的周期性行为。随着$\alpha$ 的减小，系统开始显示出更加不稳定和混沌的轨迹。$\alpha =0.9$ 、$\alpha =0.8$ 和$\alpha =0.7$ 时的相平面图显示出逐渐复杂的模式，反映了系统对扰动的不同敏感性以及分数动力学的影响。分数阶导数引入了一种非局部效应和记忆效应，显著地改变了流动和热行为，增强了牛顿流体在自由对流下的不稳定性。

比较正弦和余弦两种边界条件下的相平面图存在比较大的差异，这是因为非线性或混沌系统中，初始条件的极小差异会被不断放大，导致系统长期行为的不可预测性。

5.2. 壁面剪切应力

图4呈现了牛顿流体在两种边界速度条件$U\left(0,\tau \right)=\text{cos}\left(\tau \right)$ 和$U\left(0,\tau \right)=\text{sin}\left(\tau \right)$ 下，壁面剪切应力(以摩擦系数C_f表征)随时间的演化规律；同时通过相平面图，直观展示了不同分数阶参数($\alpha =1,0.8,0.6,0.4$ )对C_f行为模式的调控作用。

Figure 4. Effect of fractional-order parameters on wall shear stress

图4. 分数阶参数对壁面剪切应力的影响

从量化规律来看，壁面剪切应力C_f随分数阶参数$\alpha$ 的减小呈显著下降趋势；更关键的是，$\alpha$ 的降低会驱动系统的动力学行为发生质的转变：当$\alpha =1$ (纯牛顿流体)时，C_f的相平面轨迹呈现出高度规则的周期性特征，对应稳定的流动状态；而随着$\alpha$ 减小(如$\alpha =0.8,0.6,0.4$ )，轨迹模式逐渐复杂化，系统开始表现出不稳定的混沌特征。这一转变的物理根源在于：较低的$\alpha$ 会增强流体的“记忆效应”与粘弹性——记忆效应使流体的应力响应依赖于前期的运动状态，而粘弹性则会抑制流体的动量传递效率，二者共同削弱了流体与壁面之间的动量交换强度，从而降低了壁面处的剪切作用(即表面摩擦)。本质上，分数阶参数$\alpha$ 可视为系统偏离纯牛顿流体行为的“非牛顿程度指标”：$\alpha$ 越小，流体的非牛顿特性越显著，壁面剪切应力的衰减幅度也越大。

进一步对比两种边界条件的影响：$U\left(0,\tau \right)=\cos \left(\tau \right)$ 与$U\left(0,\tau \right)=\sin \tau$ 相比，结果导致更高的C_f值。这是因为$\text{cos}\left(\tau \right)$ 具有更高的初始值，并且在周期开始时变化更突然，导致更高的剪切力，从而导致更高的表面摩擦。相比之下，$\text{sin}\left(\tau \right)$ 从零开始并且变化更缓慢，导致较低的剪切力和表面摩擦。

5.3. 热传递速率

在图5所呈现的数值模拟结果中，针对牛顿流体的热传递过程，其努塞尔数(Nu)的变化呈现出显著的分数阶参数依赖性：随着分数阶参数的减小，努塞尔数呈现出同步降低的趋势。这一现象的物理机制可解释为：较低的分数阶参数会增强流体的“记忆效应”，同时提升其粘弹性特征——粘弹性的增强会抑制流体内部的动量与热量输运效率，进而削弱对流传热的强度；而记忆效应则会使流体的热响应更依赖于前期状态，进一步降低了瞬时传热能力，最终共同导致努塞尔数的下降。由此可见，分数阶参数本质上改变了系统的热力学行为模式，是调控整体热传递速率的关键参数之一。

Figure 5. Effect of fractional-order parameters on heat transfer rate

图5. 分数阶参数对热传递速率的影响

进一步对比两种边界条件下的努塞尔数变化可以发现：当边界速度条件取为$U\left(0,\tau \right)=\text{cos}\left(\tau \right)$ 时，对应的努塞尔数明显高于$U\left(0,\tau \right)=\text{sin}\left(\tau \right)$ 的情况。这一差异源于两种三角函数的初始特性与变化趋势：$\text{cos}\left(\tau \right)$ 在初始时刻($\tau =0$ )的取值为1，具备更高的初始速度幅值，且在周期起始阶段的变化率(导数绝对值)更大，这种“高初始值 + 快变化”的特征会驱动流体形成更强烈的流动扰动，从而强化了对流传热的速率，最终对应更高的努塞尔数；而$\text{sin}\left(\tau \right)$ 在初始时刻的取值为0，且初期变化相对平缓，流动扰动的强度较弱，对流传热的驱动作用不足，因此对应的努塞尔数也处于较低水平。

此外，从图5的曲线波动特征还可观察到：对于本研究中的恒定粘度牛顿流体，其努塞尔数随时间的变化呈现出明显的振荡特性，且这种振荡的幅度与频率会随着流动条件(如分数阶参数、边界速度)的改变而发生显著变化，说明牛顿流体的热传递过程对流动参数的变化较为敏感。

6. 结论

本研究考察了由时间分数阶动力学控制的牛顿流体自由对流中的稳定性。通过将分数阶导数纳入数学模型，我们发现历史依赖性是调控牛顿流体流动稳定性的关键特性，其对流动行为的时间累积效应显著影响了系统的动力学响应。主要发现如下：

1) 分数阶导数引入的非局部与记忆效应，增强了系统不稳定性。

2) 壁面剪切应力和努塞尔数随着分数阶参数$\alpha$ 的减小而减小。

3) 较低的分数阶参数$\alpha$ 值会引入更强的记忆效应和增加的粘弹性，降低对流传热效率。

4) 边界条件的初始差异会被非线性系统放大，导致相平面图呈现显著的行为分异。