伞蜥是一种生活在澳大利亚东北部沿海和荒漠边缘以及巴布亚新几内亚的小型蜥蜴。澳大利亚北部种群成年雄性个体体长0.75~1 m，平均体重400 g^[1]，生活在南纬15°附近的内陆沙漠和沿海区域之间。该地区气候炎热，旱季日气温在24~35 ℃，日均气温变化较大，降水量少^[1-2]。伞蜥的最主要特征为巨大的双椭圆结构伞状颈部褶皱(图 1)，褶皱宽度约610 mm，而颈部宽度仅有约120 mm，头和颈的最大平均长度约230 mm^[3]。褶皱具有比表面积大、厚度小、可调节性强的特点，主要作用为调节体温^{[1, 4-5]}。伞蜥经常直立行走，在光照较强时往往处于遮光环境下休息，以减少吸热^[6]。

伞蜥在持续升降体温的过程中，所获得/散失的热量等于自身产热量减去代谢散热量以及与外界环境的总换热量^{[1, 5-6]}。伞蜥的表面积与体积之比远大于其他蜥蜴，所以自身产热对伞蜥的得热/散热量影响较小。在干旱环境下，代谢散热主要来自水分代谢，且可以忽略^{[1, 4]}。而与外界的热量交换主要是与周围环境的自然对流和来自太阳的辐射。由于伞蜥在高光照下遮蔽，而在夜间无光环境下活跃^[8]，很多情况仅需考虑自然对流，因此伞蜥与外界仅涉及自然对流换热的散热方式对伞蜥散热研究至关重要。

目前对伞蜥与外界的换热研究主要是通过动物学实验从体温变化的角度进行。Christian和Bedford^[1]分别在干燥和潮湿2种条件下研究伞蜥的换热影响因素，得到黑暗条件下伞蜥在不同环境温度时的体温变化规律：伞蜥全天平均体温(32.6±2) ℃为最适体温，而在此条件下，伞蜥体温始终高于气温，主要依赖向外界散热，自身影响很小。Christian等^[4]还研究了在足够干旱的条件下，自身水代谢的散热量和进食产生的得热量分别下降86%和78%，可以忽略这种影响，而只研究与外界换热。Griffith和Christian^[2]通过实验发现伞蜥可以很好地面对水源、食物大幅增加或者减少的外界环境，在旱季通过与周围环境换热调节体温，而不依赖于食物摄入。Montandon等^[9]通过显微研究和模拟的方法，给出伞蜥颈部褶皱的细致模型(图 2(a))，即每个叶上都有3条凸起的脊，2个脊之间都有1条凸起的折痕，并通过模拟给出了脊的具体尺寸数值。Haas^[10]给出了更精确的伞蜥颈部褶皱增长模型和数学验证(图 2(c))。Castelli等^[5]通过实验给出了干旱条件高辐射状态下伞蜥体温和巢穴温度随环境温度和辐射强度变化的拟合关系式。Yun^[11]通过制作伞蜥的实体模型(图 2(b))，说明伞蜥颈部褶皱与身体的夹角不同对伞蜥高速奔跑时阻力的变化情况，得到60°~85°夹角褶皱带来的阻力逐渐增加，45°~60°阻力趋于稳定。

现有研究通过测量伞蜥的体表温度变化虽然可以得到直观数据，但不能反映伞蜥颈部褶皱在体温调节中的作用，更难以解释伞蜥褶皱大小、形状和位置等对散热的具体作用。而对伞蜥的模拟研究也仅涉及到伞蜥褶皱形态变化和动力学研究，都不能说明褶皱在伞蜥体温调节过程中的影响，所以难以分析伞蜥的形态结构对散热的影响。因此，有必要构建伞蜥模型进行数值模拟，将数值模拟方法应用于爬行动物的行为散热研究中，研究褶皱及伞蜥整体的自然对流换热，得到定量结果。伞蜥具有独特的生理学结构，生活环境使其具有很强的研究价值，通过数值模拟方法研究伞蜥颈部褶皱在亿万年进化中形成的大小形状对自然对流散热的影响，可能为爬行动物散热方式的研究提供一种传热学方法。

综上所述，本文使用数值模拟方法研究伞蜥在仅与外界发生自然对流换热的前提下，其颈部褶皱的位置(即褶皱面与躯体之间的夹角θ)以及褶皱面积对其换热影响，创造性地选用实验研究中伞蜥所处最适环境温度下的特定时刻，通过具有代表性的t_∞，t_w，构建还原模型进而计算平均表面努塞尔数Nu的具体数值。同时改变模型中的单一变量，模拟出换热结构发生改变带来的换热影响，这表明伞蜥的换热结构在长期进化中的方向可能与提高散热能力密切相关。从引入数值模拟的方法到改变内部单一影响因素进行深层次的定量求解，从传热角度定量分析在伞蜥的自然进化中，其对外界环境温度的适应性。

1 模型及计算方法 1.1 模型构建

由于伞蜥不同亚种个体差异较大，本文选择澳大利亚北部伞蜥作为具体研究对象，并选择成年雄性平均个体尺寸作为建模参考。在真实伞蜥的结构和尺寸的基础上，结合他人研究采用的模型，得出简化计算模型，如图 3所示。换热较强的部分对计算结果影响大，需要对这些部分做细化，尽可能接近真实伞蜥。褶皱作为最主要研究对象，为构建有说服力的颈部褶皱换热模型，使用文献[9]中的模型主体(文献[9-10]均通过实验对模型还原的合理性进行了细致说明)，同时参考文献[7]的建模，综合得到。其由相互重叠的平面椭圆柱以及7条突出高度接近但不等的细棱，7条棱将伞面八等分，其中中间棱较高较宽。伞蜥的头部简化为圆柱Ⅰ，直径d=12 cm，高6 cm。圆柱Ⅱ为伞蜥的颈部，高18 cm，与头部等宽。圆柱Ⅲ为伞蜥的胸腹部，直径为24 cm，高20 cm。考虑到伞蜥胸腹部对换热影响同样较大，根据文献[1]的实验研究，应该保留躯体主要部分在与真实伞蜥胸腹部体积相同的情况下简化为直径较大的同心圆柱。圆柱Ⅳ则将伞蜥胸腹部以下的部位简化为圆柱，忽略了细长尾部后半段和四肢，这些部位对换热影响很小，温度实验研究并未考虑其有换热影响，认为不影响散热，简化后可以有效节约计算成本。伞蜥模型整体高L=75 cm，模型以圆柱下底面中心为坐标原点。

伞蜥颈部褶皱化简为长轴L=70 cm，2个对称的椭圆面且相互重叠，其交线长58 cm。褶皱宽度依据Montandon等^[9]提供的经验公式计算，得到伞面厚度d₁=1.05 cm，伞面与身体间即与竖直面夹角为θ。伞面上方的竖直棱则根据脊的增长模型计算，中间最高脊按经验公式计算为高1.25 cm、宽2 cm，其余脊均按宽1 cm、高1 cm设置。

1.2 数理模型及计算方法

本文在实验[1]的基础上，选取实验条件为足够干旱、完全避免太阳光照的室内环境，忽略太阳辐射和风速，并假设蜥蜴直立且静止。根据实验结果得到伞蜥体表温度t_b和周围环境温度t_∞的关系，得到伞蜥最适温度区间为(32.6±2) ℃，模拟取室内温度恒定t_∞=29 ℃，伞蜥体表温度恒定t_b=33 ℃。假设室内空间足够大，因此边界条件不会对空间内空气的流动产生影响。模型为高温热源，经计算瑞利数Ra=6.67×10⁵，属于低速层流。本次模拟核心问题是研究伞蜥与周围环境自然对流传热能力，努塞尔数Nu表示对流强烈程度的无量纲数，因此模拟的重点是得到Nu。

结合实验进行假设，忽略内部导热，将模型设置为无滑移定温壁面，模型计算时取整体的一半。流体物性参数取定性温度t_m=31 ℃下的数值，除Z方向动量方程中的密度项外，其与温度变化关系通过布西涅斯克假设确定

$ \left(\rho_{\infty}-\rho\right)=\rho \beta\left(t-t_{\infty}\right) . $

(1)

式中:

$ \beta=-1 / \rho(\partial \rho / \partial t)_p . $

(2)

该假设需要温差很小，以减小由于热物性和体积热膨胀系数β随温度变化造成的影响。本研究的温差高达4 ℃，为确保布西涅斯克假设有效性，在Ra不变的条件下将模型等比例放大1倍使温差减小到0.8 ℃。此时空气定温t_∞=30.6 ℃，模型定温t_w=31.4 ℃，模型长度L=1.5 m。由于温差极小，热膨胀系数可按理想气体β=1/t_m近似。空气作为该研究中不可压缩的流动介质，忽略能量方程中的黏性耗散与辐射效应，数学模型涉及到的方程可以简化为：

连续性方程

$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 . $

(3)

X方向动量方程

$ \begin{gathered} u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial x}+ \\ \vartheta\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2}\right) . \end{gathered} $

(4)

Y方向动量方程

$ \begin{gathered} u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial y}+ \\ \vartheta\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}\right) \end{gathered} $

(5)

Z方向动量方程

$ \begin{gathered} u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial z}+ \\ \vartheta\left(\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}+\frac{\partial^2 w}{\partial z^2}\right)+g \beta\left(t-t_{\infty}\right) \end{gathered} $

(6)

能量方程

$ u \frac{\partial t}{\partial x}+v \frac{\partial t}{\partial y}+w \frac{\partial t}{\partial z}=\alpha\left(\frac{\partial^2 t}{\partial x^2}+\frac{\partial^2 t}{\partial y^2}+\frac{\partial^2 t}{\partial z^2}\right) . $

(7)

上述数理模型通过有限体积法离散求解，其边界条件如图 4所示。连续性方程收敛标准设置为10^-4、动量方程和能量方程的收敛标准均被设置为10^-5，具体计算设置参见文献[12]。由于研究Ra较大、边界层厚度薄、流动过程剧烈，需要非常密的网格，且模型和计算区域很大，考虑节约计算成本的同时不影响计算精度，采用非结构化网格，从几何体向外逐渐降低网格密度进行加密，每层的网格增长率为1.05，但仍需要千万量级网格量。由于模型不规则，对伞面和棱的边缘以及棱与圆柱之间存在的大量尖角做钝化处理并调整网格斜度至0.4以下，以提高网格生成质量。

1.3 计算模型与网格区域验证

研究前需要对计算模型数值方法进行验证。首先选用与此模拟相似的直立短圆柱的自然对流换热进行验证^[13]，边界条件与此模拟相同，并比较模拟得到的Nu与经验关联式计算值。对比结果如表 1所示，可知二者符合良好，最大相对误差仅在3.5%以下，证明采用该模型计算具有可信度。

因为此模型较大，有必要验证网格的加密程度是否影响到结果，以及计算区域是否足够大不影响计算结果的准确性，并避免计算区域过大而增加不必要的计算，所以需对网格和计算域进行无关性验证。从表 2可知，计算域扩大和网格数量增多时Nu变化均小于2%，说明区域和网格设置均满足要求。

2 结果及讨论

首先，将提出的数理模型按照上文所示方法应用到计算流体力学应用软件中进行数值模拟，并得到一系列结果。为方便结果讨论，对关键参数进行描述。如图 5所示，伞蜥颈部褶皱长度即褶皱模型长轴L₀，伞蜥头部直径即模型中圆柱直径d，二者比值L₀/d表示伞蜥颈部褶皱的大小。L₀/d褶皱与竖直面夹角θ是本研究的主要变量。角α表示伞面边缘与各棱之间的夹角，以下研究α=0°、45°、90°、135°和180°处在伞面上的线段(记作L_α)，其中起点设置为伞面与圆柱的相交处(x=0)，而终点在伞面外侧，定义无量纲数m=X/L_d。上述线段中，α=180°对应线段，沿z轴正方向平移1和5 mm后记作L₁和L₂，沿z轴负方向平移13.5 mm(伞面下表面下方3 mm)后记作L₃。局部表面努塞尔数为Nu_α，整个模型的平均表面努塞尔数为Nu。

2.1 温度云图及流线图

由于模型复杂，为更好地表示褶皱以及突起的脊的换热作用，选取平行于伞面的平面，高于伞面3 mm和低于伞面3 mm，可以表示出伞面上的温度场和流线分布。而整体模型则用XOZ和YOZ截面表示。如图 6所示，圆柱和伞面附近的空气被加热，由于密度差在伞面附近形成一股不断向上流动的气流，同时高温区域向低温区域传递热量。圆柱Ⅱ和Ⅲ周围形成的气流沿壁面向上流动，受到伞面阻挡后继续沿着伞面下表面流动，当流动到边缘时流向上表面。侧面进入的空气在接触到伞面后沿伞面流动并在边缘改变流动方向，在伞上表面形成空腔，由于负压增大出现倾斜向上的表面环流，反映在YOZ截面上则为椭圆环状绕流。而在图 6(e)中表现为这股气流温度较周围空气高，所以在L_0°和L_180°正上方边缘出现沿切线的温度梯度，且这段换热较强烈。在图 6(a)中则为在边缘上方的温度分布出现了突起，在伞上表面由于有棱的存在，棱周围温度高于其他同一高度处的温度，形成温度梯度，在一定程度上增强了换热。

2.2 特定位置处的温度和速度分布

为更好地表示褶皱表面附近的换热和空气流动，通过将L_180°沿Z轴方向平移得到线段L₁、L₂和L₃上的温度和速度分布(图 7(a)和7(b))以及将L_α向上平移3 mm的线段记作L′_α，得到对应的温度和速度分布(图 7(c)和7(d))。

沿伞面方向和圆柱面的距离x为横坐标。模型与环境之间的传热受到3方面的作用：圆柱、伞面和从下表面流动至上表面的气流。L₂在(0, 0.1)内主要受到圆柱Ⅰ的作用，其温度和速度分布接近单圆柱表面自然对流，即温度随x增加而降低，温度梯度降低，而壁面速度较低，在圆柱的速度边界层内速度迅速增加，之后随x增大而减少。在(0.1, 0.25)内不再受到圆柱的影响，且距离边缘较远，仅受伞面作用，然而由于L₂与伞面距离恒定，温度不变，速度逐渐降低。(0.25, 0.36)内受到下表面高温气流影响温度升高，在x=0.36 m处到达伞面边缘，此时由于下表面气流速度较大，所以在x=0.36 m处速度剧烈上升。L₃反映下表面的空气流动情况，由于向上流动的气流速度较大，被伞面阻挡后沿伞面流动，边缘流动速度依旧高于上表面，因此会造成L₁和L₂在边缘附近速度增大。L₁由于离伞面较近，沿线温度始终高于气流温度，故温度在边缘才会下降。α=0°和180°在对称面上，而α=45°、90°和135°则在3条棱上，且为各棱的中线。总体温度和速度沿线变化趋势与之前相同，但需注意在棱上的线段，由于棱同样是高温热源，附近温度较高，线上平均温度较高，且高于热气流温度，边缘温度不会升高但速度增大。中间棱高度更高，速度更大。

2.3 Nu_α分布

Nu_α反映褶皱表面每一点的温度变化，由表面局部温度梯度决定且沿扇面表面连续变化，受表面温度梯度和流动速度影响。该参数可以让我们进一步了解流体的流动传热特性。为更直观地说明模型表面的换热情况，取在L_α上Nu_θ的沿线变化。

如图 8所示，α=0°和180°在扇面的对称面上，换热较内表面剧烈，其中L_0°上Nu达到27.825，是整体换热的2.1倍。从圆柱表面起，圆柱和棱的交接处温度等于模型温度，温度梯度从0开始，随x的增加，换热剧烈，温度梯度大，Nu_α先急剧上升，到距离圆柱约0.15 m处达到最大后，圆柱影响越来越弱，周围温度梯度和速度都减小，Nu_α减小，但是在边缘后受到热气流作用，换热变剧烈，Nu_α继续增加。但在棱上的空气，靠近圆柱时受到相邻棱和圆柱的共同加热，温度梯度较小，Nu_α低于在对称面上的，同样受到热气流作用后换热加强，Nu_α增加。通过Nu_45°和Nu_135°的曲线可以发现椭圆虽然不关于α=90°对称，但不对称影响很小。

2.4 Nu分布及θ和L₀/d对其影响

Nu反映整个几何体上的平均传热速率，直接反映模型的传热能力，通过积分模型表面上的Nu_α来获得，是本次模拟最核心结果。由图 8知伞面边缘以及α=0°和180°处，Nu较大，因此伞面上的Nu大于整体Nu，褶皱的散热能力更强。根据模拟，分别求解模型伞结构的Nu=16.236，而除伞外Nu=14.784。此外褶皱表面积是其他身体部位的2.36倍，故褶皱换热比重很大。考虑伞蜥站立时颈部褶皱与身体呈一定夹角，以及同种群伞蜥褶皱长轴与头的直径的比值也接近，改变θ和L₀/d可以得到不同的Nu，因此可以研究伞蜥颈部褶皱对自然对流传热能力的影响。

图 9为Nu和θ变化关系。模拟时首先在θ=45°~85°范围内每5°取一个工况点进行数值计算，得到一个初步的Nu-θ变化规律：Nu随θ增大而增大，到达65°附近的极大值后随θ增大而减小。由于改变θ后模型的特征长度不变，根据Nu计算式Nu=hl/k，传热速率有且仅有唯一极大值，即伞蜥在不同θ下有最大的传热速率。随后，为更精确地得到最大传热速率对应的角度，在60°~70°之间，将θ调整为63°和67°重新做2组模拟，得到极大值点在63°~67°之间。同理再模拟64°和66°，将得到的结果进行拟合，可以得到Nu-θ的拟合曲线。同时因为模型形状大小均未发生变化，由Nu_90°和Nu_180°随θ的变化得到Nu_α与X轴围成的面积，可以反映整体模型Nu-θ的变化规律。如图 9(a)和9(b)所示，Nu_90°和Nu_180°曲线与X轴围成面积均在θ=65°附近最大，反映Nu在θ=65°附近有最大值，此时传热速率最大，伞蜥与周围环境自然对流散热作用最强。

3个亚种的伞蜥生活环境的温度条件不同，其体型虽然相当但是颈部褶皱大小差别很大，所以环境温度可能对伞蜥褶皱的面积有重要影响。改变扇面大小，分析在不同的伞面面积下自然对流传热速率的变化。固定褶皱以外的部分不变，等比例改变伞面尺寸，使扇面长轴L₀每次变化10%，模拟过程和收敛标准与之前一致，得到在不同L₀/d下Nu的数值。由于改变伞面面积会直接改变特征长度，根据Nu/l=h/k，k为流体热导率保持不变，所以不同L₀/d对应Nu/l的变化反映了对流传热系数的变化，同样通过不同L₀/d下Nu_90°随x的变化反映Nu的变化。不同模型与周围空气的自然对流换热系数有且仅有一个最大值，在误差允许范围内，认为6.3≤L₀/d≤6.6时对流换热系数最大，此时模型与外界自然对流换热能力最强。如图 10(a)所示，棱长增加不会改变Nu_90°的变化趋势，但是传热速率较大的区间长度几乎不变，仅所取温度截面上棱的平行线与扇面上棱之间的距离增加，且随棱长增加边缘处的Nu降低，故随x增加Nu_90°单调递减。随着褶皱面积的增大，Nu变化趋势同理，所以Nu逐渐减小且变化幅度较小，但是褶皱表面积增大速度远大于体积增大速度，所以特征长度减小速率更快，因此Nu/l-L₀/d有且仅有唯一极大值。

通过之前的模拟发现，改变θ以及改变L₀/d，模型与外界环境之间都存在最大的对流传热速率，意味着在θ=64°~66°和6.3≤L₀/d≤6.6这样2种情况下伞蜥都存在最佳散热情况。为说明这2个最大值的意义，与文献中澳大利亚北部伞蜥种群的成年雄性个体统计数据进行对比，得到结果如表 3所示。

由表 3可知，种群之间个体虽大小不同，但是伞蜥种群的实际θ和L₀/d都在一个正态范围内，可以认为对于文献[1, 4]中所研究的个体数量400的伞蜥种群而言，其成年雄性个体褶皱的平均大小和自然状态下的平均角度是对于这个种群而言最正常的数值。而模拟的结果与同类别种群的平均值比较接近。对于θ，模拟出的最佳换热对应的角度与实际的平均值相差10%，且低于大部分个体的数值。认为可能有以下原因：其一，伞蜥的褶皱角度自身具有调节能力，所以个体不一定均在自然状态下；其二，虽然前文已充分说明忽略部分身体结构和忽略自身产热和代谢的必要性和合理性，但是不可避免产生偏差。同样，对于L₀/d，由于模拟自身精度仅能确定一个范围，这个范围可以认为与平均值接近。对于2组模拟都可以在误差允许范围内认为伞蜥颈部褶皱θ和L₀/d的平均值是伞蜥散热效率最高的情况，所以伞蜥平常的褶皱角度和伞蜥褶皱的大小可能提高伞蜥与外界环境的自然对流散热效率。因此可以认为在亿万年进化中，伞蜥褶皱的面积和张开角度可能都受到自然对流的较大影响，由于其生活环境接近赤道，全年温度较高，伞蜥的褶皱体现了对高温环境的适应性。同时自然对流可能是伞蜥进化方向的重要影响因素之一。

3 结论

本文基于数值模拟的方法，对伞蜥的独特生物学结构颈部褶皱在不同θ和L₀/d下自然对流散热进行模拟研究，并根据模拟结果拟合得到对流换热强度随θ和L₀/d的变化关系。结果表明：

1) 模型在L_0°和L_180°上换热最为剧烈，L_0°上Nu达到27.825，是整体换热的1.8倍。同时边缘热气流造成换热加强，使得伞面Nu达到16.236，是整体的1.045倍。因此形成了更强的换热。

2) 伞蜥在只有自然对流条件下，张角在65°附近时得到最大Nu=15.52，此时具有最大的换热能力，接近伞蜥在自然状态下的褶皱角度(65°~75°)，与所研究种群平均值差10.5%。同时伞蜥在仅有自然对流条件下，褶皱长轴与头的直径比在6.3~6.6之间有最大的对流传热系数，与所研究种群平均值差小于5%。该结论说明伞蜥的自然进化是有利于向更高升/降温速率的方向进行，从而提高伞蜥对自然环境的适应能力。

3) 通过2组模拟，可以认为在伞蜥的进化过程中，伞蜥褶皱的大小和与身体角度均朝有利于增大伞蜥与环境的自然对流换热速率的方向进行。自然对流换热可能是决定伞蜥进化方向的重要原因之一。