滨海海底边界层对底栖微体化石仿真模拟的影响

Effect of boundary layer on simulation of benthic microfossils in coastal seabed

摘要

计算流体力学CFD(computational fluid dynamics)近年来在古生物学研究中已有大量应用。该方法在研究生物化石的个体形态、器官功能,以及其与生存环境之间的关系方面具有重要意义。目前古生物领域的CFD仿真对象多为厘米级别体型的生物,而对于毫米级别大小的微体生物的仿真则相对较少。区别于厘米级别的宏体生物,底栖型微体生物的生活环境局限于海床上方垂直高度较小的区域,受海底黏性边界层低流速区域的影响更为明显,因而在设置这类化石仿真的水流环境时应当考虑对于边界层流域水体流速以及海床表面地理环境的还原。本文基于随机表面生成的方法,构建了不平坦的滨海海床表面模型,并模拟了水流在海底边界层附近受海床表面地形影响的流动状态;基于获得的仿真流速数据,我们仿真了寒武纪微型底栖刺细胞动物化石Quadrapyrgite在起伏的海床表面环境中四个不同位置的受水流阻力情况。仿真结果显示,近海床表面处的水体流动可形成明显的低流速区域,且低流速区域的厚度随着入口速度增大而变薄。在起伏表面的凹凸区域中,凸起地带的迎水面形成的低流速区域较薄,流速变化较快;而在凸起地带背水面以及凹陷地带形成的低流速区域较厚,流速较缓且易形成涡旋。流速的差异体现在不同位置的Quadrapyrgite的受阻力大小不同上,其差值可达数倍到数十倍。这种量级的受阻力差值可以在毫米级到厘米级尺度上塑造、影响底栖型生物的分布甚至其自身的摄食行为。本文的研究为CFD方法模拟底栖微型化石的生存环境提供了更加深入的思路和方法。

关键词

底栖微体化石 / 计算流体力学 / 海底边界层 / 傅里叶变换 / 随机表面产生方法

Key words

benthic microfossils / computational fluid dynamics / bottom boundary layer / Fourier transform / random surface generation

0 引言

计算流体力学(computational fluid dynamics,下文简称CFD)是传统的流体力学与数值计算科学相结合的一门交叉学科,其核心思路是利用有限元(finite element method,简称FEM)、有限体积(finite volume method,简称FVM)等方法,对求解域进行离散,并求解流体力学的控制方程以模拟其中的液体或气体等流体流动。随着计算机硬件与软件技术水平的发展,CFD方法已成为工程中成熟的研究技术之一,在许多学科领域已有广泛的应用,如建筑、材料、汽车以及航空航天等。同时,CFD技术也大量应用于现代生物的仿生学研究,如鱼类^[1-2]、水母^[3⇓-5]等。

基于生物化石形态学重建的CFD仿真更是近年来古生物领域中新兴且重要的研究方向之一,尤其是用于验证生物化石的功能学与个体生态学的猜想^[6]。例如:Waters 等^[6]对棘皮动物海蕾有无鞭毛驱动的呼吸状态进行了仿真,表明在鞭毛驱动下的内部水流流速维持在外界水流流速的一半时能达到最佳换氧效率;Rahman等^[7]对棘皮动物海扁果类(cinctans)中不同体型的6个种进行了三维重建与流体分析,表明具有更大背部凸起的海扁果相较于扁平的类型易于将食物随水流引导至口部,而后者相较于前者则更易于降低阻力保持稳定性,体现了生物在进化中对于稳定性与摄食能力的权衡;Song等^[8]对队列状排列的三叶虫化石进行了重建与仿真,证明了三叶虫的队列行进行为更可能是为了减少个体所受的阻力,有利于其顺利完成迁徙活动。上述仿真研究主要从水动力学角度分析了宏体动物的群居或聚集行为以及个体形态学对其适应环境能力的影响,在这类仿真中,建立物理场模型时需着重考虑生物周围的流场环境^[9],因而其在设定仿真物理场的流速参数时,多参考现代海洋中测定的流速数据(如采用一组或多组固流速参数进行定量的模拟分析),但并未考虑实际海床沉积层表面的形态对于底栖生物所处的水流环境的影响。而研究表明,滨海区域(图1A^[10]),脉动较为频繁的海流以及波浪作用(图1B^[10])使得海床表面沉积质颗粒再悬浮于水中并被搬运,并在海床表面产生可阻碍水流的起伏与波痕,进而影响海床表面附近的微观流速分布。这种地形的作用导致的流动变化应当对底栖微体生物的生存环境有着不可忽视的影响。

在流体力学的概念中,流体在较高雷诺数的湍流状态下流过无滑移的壁面时,由于无滑移壁面的黏性力作用,壁面上方会形成具有流速梯度的薄层,称为边界层(或湍流边界层,boundary layer);而边界层中靠近壁面处的流域主要由黏性力主导,其中相对上方的区域流速较低,该流域被称为黏性子层(viscous sublayer)(图1C)。海洋环境中,在风力、洋流以及潮汐力等作用下,海床沉积物上方也会形成与工程中流体力学边界层有相似特点的边界层^[11]。在滨海区域(图1A^[10]),近海床表面的上方存在1 cm左右厚度的波流边界层(wave-current boundary layer,为便于理解与叙述,后文称之为海底边界层),层中的流动受到震荡波流(图1B^[10])以及壁面黏性力的阻碍而被减速(图1D^[11])。对于生活在海底边界层中的毫米至厘米级别的小型底栖生物来说,其生态活动与该区域的环境密切相关,如水流速度、温度、氧气含量、沉积物浓度等。早先学者们针对海底边界层的流速特性已进行过一些定量的研究以及实地的测定,如:Caldwell和Chriss^[12]利用水下热探头实地测定了Oregon大陆架200 m深的海床表面上方的流速分布情况,并通过线性拟合得到海床表面形成的边界层中黏性子层的厚度约为6 mm;Wengrove等^[13]测定了潮下带区域约1.5 m深的海底边界层流速分布情况^[8]。此外,海底边界层的厚度也会受到其他因素的影响,如区域内生物的高度、富集程度、海床表面的粗糙程度以及上方洋流的速度等^[7,14]。

在古生物学领域,Liu等^[15]基于滨浅海环境的流速参数对寒武纪微体的底栖型刺细胞水母化石进行了流固耦合仿真,模拟了两种具有不同围鞘形态的化石在水中的受力情况,表明两辐射对称的六方锥石(Hexaconularia)适应大陆架水流环境的能力强于四辐射对称的四方塔型壳(Quadrapyrgites),并通过对比二者不同的CFD效应尝试解译寒武纪幸运期-第二期之间微体底栖型水母的生态演替关系。该研究在展示了CFD技术在微体生物化石领域应用的可行性同时,也初步提到海底边界层(黏性子层)对于底栖型微体生物的生存环境的可能影响。由于寒武纪微体化石绝大多数生活在浪基面之上的滨海环境,这一区域海床受到潮汐、陆源物质输入等因素的影响,表面地理环境相比浪基面之下浅海环境更加复杂,现有的较为微观的定量仿真虽可从形态学与水动力学的角度分析这类生物适应环境的能力,却难以还原其原本生活的海洋与地质环境,这给利用CFD方法进一步研究底栖生物的生态活动及其受所处的宏观流动环境的作用造成了局限。而受限于计算资源以及软硬件技术,目前古生物领域将实际海床表面以及边界层环境与底栖微体生物相结合的宏观流体仿真研究相对则存在一定空白,但毫无疑问,模拟底栖微体生物实际生活的地理以及水流环境(即海底边界层流动环境)并开展相应仿真研究,不仅可以提高仿真本身的还原度,亦对研究者们更为清晰地认识底栖微体生物对环境的适应能力有着重要意义。

基于Liu等工作中的部分研究方法和讨论^[15],本文利用COMSOL软件及其官方资料库中提供的随机表面生成方法,构建了具有凹凸起伏的海床表面,并用湍流模型模拟了洋流以不同速度流过海床表面时形成的边界层流态,主要对其中的流速进行了可视化;其次,基于这部分仿真得到的流速数据,本文又根据随机生成的海床表面模型构建了类似真实海床表面的具有起伏的不规则仿真域,并以寒武纪微体底栖刺细胞动物四辐射对称的四方塔型壳(Quadrapyrgites)化石作为对象开展了流体仿真,分析了其在凹凸区域中受阻力情况的差异。仿真结果说明:海底边界层在迎水面的上坡区域相对较薄,流速变化较快,流速相对较高;而在背水面的下坡以及凹陷区域边界层较厚,流速较缓,且易形成涡旋。这一结论同样体现在Quadrapyrgites的受力分析情况中,在相同的水流环境下,固着在迎水上坡区域生长的Quadrapyrgite所受的阻力显著高于固着在凹陷或背水区域生长的Quadrapyrgite。这一结论可作为已有微体化石仿真结论的补充,对于认识和推测微体生物生存的环境有着重要意义。此外,本文首次实现了将底栖微体化石与模拟真实海床环境的仿真域相结合,为微体化石的仿真研究提供了参数设置上的参考,也为后续更加精确地开展宏观的古生物流体仿真研究提供了技术依据。

1 仿真原理与方法

本文利用COMSOL Multiphysics 5.6版本软件来建立模型并开展仿真。COMSOL Multiphysics软件中提供了建立多种不同几何形体的方法,亦可导入外部的模型文件,用户可快速生成基本几何形体并在此基础上构建仿真域。COMSOL官方的资料库也介绍了利用三角函数生成参数化随机表面的方法(http://cn.comsol.com/model/generation-of-random-surfaces-50281)。本文将这种方法进行简要概括,并根据这种方法来还原海床表面凹凸不平的环境,构建虚拟的海床表面仿真模型。第二部分仿真中采用的Quadrapyrgite模型与Liu等研究^[15]中采用的模型参数基本相同,为确保后续网格划分的稳定性,本文进一步对模型的表面细节进行了修复,并将模型以.stl格式保存,以便后续将其导入软件之中。

1.1 随机海底表面构建的基本原理

自然界中多数周期信号可视为由数个余弦(或正弦)信号叠加而成,这即是工程应用中傅里叶变换的核心思想。因此,一个二维随机表面可用多个不同频率的余弦信号叠加近似表示出,如公式(1)所示。式中A为比例因子,用于调整生成的信号的幅值,即随机表面凹凸部分的高度与深度;N为最大谐波次数,该数值的大小可以决定信号的随机程度;φ_n为第n次谐波的相位值,用一个在 \left[-\frac{\pi }{2},\frac{\pi }{2}\right]上的均匀分布的随机函数u(n)表示。调整上述的参数值即可获得不同样式的随机表面(图2)。同样利用傅里叶变换的原理,将公式(1)扩充一个维度,并将随机函数u(n)变为u(m,n),即可生成一个三维随机表面的横截面,如公式(2)所示。

为了使随机生成的表面更能还原自然界中被风、水等侵蚀的样貌,还需增补一个侵蚀函数,该函数如公式(3)所示。公式中β用来调节侵蚀的程度,该值越大,则随机生成的表面越平滑,反之则越粗糙;式中g(m,n)为二维标准正态分布的随机函数。将公式(2)与公式(3)整合后,最后的生成函数如公式(4)所示。

通过改变式中的A、β、N值,即可以得到不同粗糙程度和起伏高度的随机表面(图3)。在COMSOL中,定义相关的随机函数以及相关的A、β、N变量后,即可利用公式(4)在“参数化曲面”节点中生成不同样式的随机表面。

(1) \begin{array}{c}f\left(x\right)=A\stackrel{N}{\sum _{m=1}}cos(2m\pi x+{\phi }_n)\end{array}

(2) \begin{array}{c}f\left(x\right)=A\stackrel{N}{\sum _{m=1}}\stackrel{N}{\sum _{n=1}}\left[cos\right(2m\pi x+2n\pi y)+u(m,n\left)\right]\end{array}

(3) \begin{array}{c}E(m,n)=\frac{g(m,n)}{(m^2+n^2{)}^{\frac{\beta }{2}}}\end{array}

(4) \begin{array}{l}\mathrm{f}\left(\mathrm{x}\right)=\mathrm{A}\stackrel{N}{\sum _{m=1}}\stackrel{N}{\sum _{n=1}}\left[E\right(m,n\left)\mathrm{c}\mathrm{o}\mathrm{s}\right(2m\mathrm{\pi }x+2n\mathrm{\pi }y)\\ +u(m,n)]\end{array}

1.2 仿真域的构建

基于1.1节中的方法,本文设定A=0.008,N=20,β=1.5,构建了用于仿真的1 m×1 m随机表面,生成的随机表面的最大凹深和凸高分别为-0.04 m和0.02 m。为了确保仿真的可复现性,g(m,n)与u(m,n)两个随机函数在COMSOL中均设定随机种子参数为3,固定该参数的取值可使得每次重新构建模型时生成的随机表面均相同(否则每次重新构建模型时均会生成不同的随机表面)。仿真域的构建如下面步骤所示:首先创建0.8 m×0.4 m×0.4 m(长,宽,高)的长方体;其次将长方体与生成的随机表面的中心区域重合,并作布尔操作去除掉多余的区域,即可获得具有不平坦下表面的海床表面仿真域;最后将仿真域的一端延长至5.6 m,延长区域的长度为仿真域的五倍以上,使水流从延长域流入仿真域时得以充分发展,以形成较为稳定厚度的边界层(图4A)。

基于上述的海床表面仿真域,本文取沿其中心轴的切面截线,并分别将其向两侧拉伸0.005 m,形成一个海床表面仿真域近似的部分切片。为了降低计算开销以及仿真中所用的生物模型与整个仿真域的尺寸差距,该仿真仅选取该切片中具有地形代表性特征的一部分区域作为仿真域(包含了背水面、凹陷、迎水面区域),最终形成的仿真域为长0.2 m、宽 0.01 m、最大高度约 0.025 m的不规则区域(图5A)。Quadrapyrgites模型(图5B,高约0.002 m)被分别置于底面的四个不同的位置(如图5C,下文以Q1、Q2、Q3和Q4代表,从左到右分别对应迎水面区域、浅凹陷区域、深凹陷区域和背水面区域的位置),低于底面的区域同样采用布尔分割方法去除,以形成完整的仿真模型。上述的两部分仿真域中,水的密度和动力黏度设置为1 000 kg/m³和0.001 Pa·s。

1.3 边界条件与网格划分

对于海床表面仿真域,其延长域的一端被设置为入口边界条件,设置入口速度时需使得水流经延长域到达计算域时的流速垂直分布与文献中实际测定近海底边界层的流速数据相同。根据文献中不同的海底边界层测量数据^[12-13],仿真中采用的入口流速分别为0.04 m·s^-1、0.06 m·s^-1、0.08 m·s^-1和0.10 m·s^-1。计算域的一端被设置为零静压的出口条件。整个计算域和延长域的下表面被设置为无滑移壁面(no-slip boundary),用以模拟实际中的海床表面;两个侧面和顶面设为滑移面,用于模拟水域环境(图4A)。根据不同的计算精度需求,仿真域与延长域划分了不同细化程度的网格,以对仿真域进行离散;前者的网格细化程度相对较高,后者细化程度则相对较低(图4A)。整个区域采用自由四面体网格覆盖,并在无滑移边界与流体接触的面添加12层边界层网格,以更精确捕捉边界层中流速变化较快区域中的流动(图4B)。对于该部分仿真,本文开展了网格敏感性测试,用以选择计算量与精度相平衡的网格参数(见后文2.1节)。

对于固着有Quadrapyrgites的仿真域,其边界条件划分与上述海床表面仿真域的相似。该仿真域的入口流速设置采用第一部分海床表面仿真域入口速度为0.04 m·s^-1的仿真结果中与该入口相同位置的流速采样点数据(利用三次样条插值法拟合成与高度z相关的函数)。仿真域的网格细化参数采用“超细化”的预定义设置,最后生成的网格共有3 916 972个单元。

1.4 流场模型与求解器

流场的模型选择依赖流体力学中的无量纲准数雷诺数(Reynolds number,记作Re)大小,其表征着流场中惯性力与黏性力之比。对于海床表面的仿真域,根据雷诺数计算公式^[16],在不同的入口流速下计算得到的当地雷诺数Re_x范围为224 000~560 000。考虑到雷诺数的范围以及近壁面区域流动的精确计算,以平板流动问题中Re_c=5×10⁵作为临界雷诺数,本文选用低雷诺数壁处理方法的SST模型来描述流场中的控制方程,该模型结合了k-ω和k-ε模型的特点^[7],可以适用于较低雷诺数到高雷诺数的流场,亦可以较为精确地求解壁面附近区域的流动。对于固着有Quadrapyrgites的仿真域,本文根据底栖生物生活的主要流场环境资料^[9],并参照Liu等的研究^[15]中的设置,采用层流模型(Re范围约为1至10之间)。最后,所有仿真均采用直接、全耦合求解器计算流场的稳态解。

2 仿真结果

2.1 海床表面仿真的网格敏感性测试

对于海床表面仿真,共采用四组不同细化程度的网格参数的网格进行敏感性测试,网格单元总数范围为527 393~5 164 041。测试时以入口速度为0.06 m·s^-1的模型作为测试模型,在尽可能靠近模拟海床表面的位置选取了149个采样点的流速数据,并比较每组仿真数据的差别。经测试,本文最终采用的网格共有3 049 620个单元,与最细化网格参数设置下的计算结果的平均差值约为2%,精度可以满足仿真需要,因此后文的仿真均以该网格参数进行。

2.2 结果处理与可视化

下文实验的后处理共包含两部分,分别针对两组不同的仿真。其中2.2.1节与2.2.2节基于海床表面仿真结果,分别对边界层表面的流态进行了俯视面与纵切面的可视化,以展示模拟的海底边界层区域的宏观流速分布以及近海床表面的微观流动和范围;2.2.3节实验基于固着有Quadrapyrgites的仿真结果,对过化石模型对称中心纵切面的流速范围进行可视化;同时,分别计算了Q1-Q4位置四个化石所受水流阻力的大小并进行了汇总。

2.2.1 边界层低流速区域的分布范围

为了展示边界层中的流速随地形起伏的变化,图6以俯视视角分别绘制了虚拟海床表面的高度与流速箭头。图中的红蓝颜色表示海床表面的高度,而彩色箭头的大小与颜色用以表示区域中各处的流速。绘制流速箭头时设置显示范围为0~0.01 m·s^-1,以描绘近海床表面低流速区域的流态,因而箭头密集的区域说明该处低流速区域较厚,反之则说明该处低流速区域较薄,水流速度梯度相对较大,变化较快。

仿真结果显示,区域中的箭头数量随着入口流速增加而逐渐减少(图6A-D),在入口流速为0.04 m·s^-1(图6A)时数量最多,0.1 m·s^-1时(图6D)数量最少,说明低流速的区域范围随上方水流速度增大而逐渐减小。此外,每组仿真中箭头主要集中分布在地形中较为低洼区域,说明地形较为突出的区域表面流速较快,而较为凹陷的区域相对流速较慢,即使入口速度增大,也基本保持了这一趋势。

2.2.2 凹凸区域中边界层的厚度以及流动形态

边界层区域中的流动以沿仿真域中心纵向作截面的方法进行可视化(图7A-D,即左列图,分别展示了0.04~0.1 m·s^-1入口流速下的边界层分布)。由于高于海床表面约6 cm以上区域的流动受地形起伏影响较不显著,且靠近海床的位置流速随高度变化更加明显,此处为了便于叙述与观察,四组仿真均以流场中流速低于0.01 m·s^-1的区域的边界层厚度来说明在床表面不同位置的分布的差别。入口速度为0.04 m·s^-1时,流速低于0.01 m·s^-1的区域厚度约为0.5~1.5 cm,而在凸起区域的迎水面厚度则只有不足0.1 cm(图7A)。随着入口速度增大,背水面以及表面的凹陷处的反向回流逐渐明显,在入口速度为0.1 m·s^-1时流速低于0.01 m·s^-1的区域厚度则分别降至为0.3~1.3 cm以及0.02 cm。这可以说明,边界层区域中的总体流速也在随入口速度增加而增大,对应着近海床表面定量流速区域的边界层厚度随之减小,但边界层中总体的流动形态基本未发生明显变化,相对保持一致。

图7中右列图展示了四组仿真中水流在海床表面凹凸区域中边界层内的微观流动情况(左列图红框区域的放大图)。可以观察到,近海床表面的边界层的厚度随地形起伏有所变化:在凸起区域的迎水面流速相对较高;而在凹陷处或凸起区域的背水面流速较低,且形成较大范围的低流速区域。此外,迎水面上方的水流相对平行于海床表面,流动较为平稳;而凹陷区域形成的低流速区域中可以观测到回流与涡旋。

2.2.3 凹凸区域中的生物仿真结果以及受力情况

在固着有Quadrapyrgites的仿真域,由于所选流场模型的差别,仿真域中的流速总体低于海床表面仿真域中的结果,但流态大致相近,水流在壁面上方的流动也具有类似海底边界层中的流速梯度(即靠近海床表面区域的流速低于远离海床表面区域的流速),并形成了部分涡旋和回流区域(图8A)。而不同位置Quadrapyrgites所处的流场环境流速与回流和涡旋的作用均有关。从可视化结果可以看出,Quadrapyrgites在Q1所处的迎水面位置显然流速最大(图8B);在Q4所处的背水面位置流速最低(图8D),其周围的流动呈平行于海床表面的反向回流状态;而Q2、Q3处于凹陷区域中形成的涡旋流域(图8C)中,环境流速相对高于Q4但低于Q1。虽然所处的高度和横向位置略有不同,但Q2和Q3所处的流场环境可视为近似相同。

上述可视化结果同样也体现在Q1-Q4 Quadrapyrgites 所受平行于水流方向阻力的计算中(表1)。阻力大小的绝对值表征流场中物体受到水流冲击力的大小,其正负与参考水流方向有关,与参考水流方向相同则为正值,反之则为负值。由于Q2、Q3和Q4的位置主要为涡旋或回流区,受反向回流的冲击作用,因此所受的阻力体现为负值,但其绝对值仍明显小于Q1所受阻力的值。其中,Q4在仿真域中与Q1的相对高度相差最小,但所受的阻力最小,仅为Q1的2.3%左右;Q2和Q3所受的阻力高于Q4但远低于Q1,分别约为Q1的9%和7%。由于Q2与Q3同处凹陷区域形成的涡旋中,因此二者所受阻力大小的差别相对不显著。

3 讨论

3.1 边界层模拟结果与实际数据的等效性

本文参照实际海底设置的等效仿真主要是针对底栖微体化石生存环境的仿真,因此,四组实验中,水流经延长域流至仿真域入口时的流速可与实际测定的近海底边界层流速数据进行对比,以说明参数设置的合理性。入口流速为0.04 m·s^-1与0.06 m·s^-1下的仿真域入口流速数据与Wengrove 等^[13]测定的近海底边界层流速相对较为吻合,而入口速度为0.1 m·s^-1下的仿真域流速数据与Chriss和Caldwell^[17]测定的海底边界层数据较为吻合,因此,本文设置的流速参数范围能覆盖实际测定的数据,具有合理的现实依据,等效适用于底栖生物生存水流环境的讨论。

3.2 边界层可视化与生物所受阻力

边界层可视化的结果说明,横向的水流流动在海床表面上方形成了厚度不同的低流速区域,即海底边界层区域,其中的流速分布受横向水流的速度与地形环境作用的共同影响,流速低于0.01 m·s^-1的区域随入口速度的增大而减小的同时,也随着海床表面凹凸区域的变化有所减增。由于包括Quadrapyrgites 在内的寒武纪早期底栖微体生物高度大多在几个毫米以内,且其生活环境一般局限于海底表层以及不同深度的底质中,因此都处于海底边界层的范围之内。即使上方横向水流的速度增大,由于边界层中流速分布在垂直方向上的梯度特性,边界层内的相对流速仍低于上方流动,因而可为底栖微体生物的生态活动提供相对稳定的低流速环境。同时,受海床表面地形起伏的影响,边界层中的微观流速分布也存在差异,在迎水面上坡的区域,水流方向几乎平行于坡面,垂直于海床表面方向流的速度梯度较大(即流速在该方向变化较快),流速也相对较高;而在背水坡面与凹陷处区域,低流速区域厚度较大,流速相对较低,且流动较易形成涡旋与回流。因此,即使在洋流环境相近的海域,海底表面的起伏环境也会在海底边界层中产生数倍的流速差。流速差异除对海底边界层区域中的水流形态以及食物颗粒与其他沉积物的掺杂程度产生影响之外,产生的最直接的结果即是影响固着在海床表面生长的底栖微体生物的受阻力大小,进而影响其生态活动的稳定性。

在对不同位置的Quadrapyrgites的仿真中,计算得到Q1-Q4间的阻力差值也体现了这一结论。迎水面位置Q1处所受阻力远大于Q2、Q3和Q4处的,差值最高可达40多倍(与Q4处相比)。导致Quadrapyrgites在Q1处受阻力较大的原因主要是流体流过迎水面位置的上斜坡表面时产生了加速作用;Q4处的背水面位置为流体流动受凸起区域阻碍而产生的低速区,因而此处的生物个体所受阻力最低; Q2和Q3处凹陷区域中存在的涡旋使得该处的流动速度相对不稳定,因此该处的水流产生的冲击力虽高于Q4区域,但仍远低于Q1区域。综合两部分仿真的结论,可以说明海底边界层区域的低流速环境对底栖微体生物提供了一定程度的低阻力保护作用,但微观地形的差异亦使得这种保护能力在不同位置有所差别,可能达到数倍至数十倍的量级。相比之下,可认为海床表面中微观的凹陷环境更有利于Quadrapyrgites这类底栖微体生物开展摄食等活动。

3.3 边界层对底栖型生物生存环境的影响

大陆架的海洋环境是现代海洋生物最重要的栖居地,因而其中的水流与掺杂的沉积物搬运和沉降不仅塑造了海底地形,且对海洋底栖生物的呼吸、摄食、繁衍等活动起着至关重要的作用^[14]。而生活在滨岸环境中的生物更容易受到波动水流的影响,如现生底栖生物有着较为稳固的固着方式(如固着盘)或者依靠主动游泳能力适应较为复杂水流环境。有别于现代具有固着盘的底栖生物,Quadrapyrgites不具有固着盘,也不具有左右移动能力,其呈辐射对称的比较薄而柔软的围鞘使得其无法采用更有利的姿态面对水流,因而其对水流冲击带来的阻力作用更为敏感^[15]。因此,即使Quadrapyrgites毫米级别的体型使得其总体受到的阻力远小于更大型的底栖生物,但较高的流速环境仍可能使得其发生倾伏,因而失去从其周边水流中高效获取食物的机会。而海底边界层中适宜的流速环境有助于沉积物和食物颗粒的沉积,利于提升底栖型生物的摄食效率^[18],也可减小波流对底栖生物的冲击^[9]。此外,仿真的可视化结果也说明,海底边界层流过起伏的海床表面时,在凹陷地形中产生了大量涡旋区域(对应仿真域中Q2和Q3位置)。涡旋的产生提供了相对较为适中的流速与阻力环境,有利于水流中沉积物颗粒的掺杂,也有助于海床表面沉积物颗粒的再悬浮,因此极大地影响了底栖微体生物的食物丰度。由于滨海区域的海流受潮汐力等因素作用而脉动频繁,且其中的海床表面附近的流动受海平面风力大小、温度等因素影响更为显著^[19],因此海床表面的沉积物易受海水流动的搬运^[20]和侵蚀作用,进而形成大量凹凸起伏的波痕地貌。可以推测,因波痕而形成的涡旋环境可以影响营养物质存留,并可进一步影响该区域生活的底栖微体生物种群的规模。

除地形起伏作用以及上方水流的流速的影响外,边界层的流动状态也与海床表面的粗糙程度以及区域内的生物富集程度等因素有关^[7,14,21]。较为粗糙的海床表面可以在增加边界层范围的同时,也会在海底边界层流域内制造更多涡旋^[22]。此外,底栖生物营群居生活可降低群落中单体受到的阻力,实际上提高了群落区域海床表面的粗糙度,并在其中的生物单体之间产生更多涡旋区域,更降低了单体所受的阻力,也影响了营养物质的摄入效率^[21,23]。因此,海底边界层与底栖生物群之间可认为是相互作用并相互影响的,进而改变了底栖生物的生存环境。而较大规模生物群落对于海底边界层流动的宏观影响程度则需要更多计算资源以及进一步的仿真说明,因而这部分内容在本文中暂未涉及。

4 结论

本文利用随机表面生成的方法构建了自然中具有起伏的海床表面,并模拟了洋流流过海床表面时形成的近海底边界层的流动状态以及其中分布于不同位置的生物的受阻力情况。仿真结果说明:

(1)边界层中的流动受上方洋流速度与地形起伏的共同影响。

(2)边界层中的低流速区域在凹陷区域或凸起区域的背水面较厚,在凸起区域的迎水面较薄,因而前者的水流速度相对较低,而后者的水流速度相对较高。这种流速的差异可能对其中底栖微体生物的生存与分布位置产生影响。

(3)边界层的低流速环境可对毫米级别的微体生物提供相应的低阻力保护作用。而这种保护作用受海床表面凹凸起伏的影响,对迎水面位置、凹陷位置和背水面位置的生物个体造成的阻力差可达数倍至数十倍。

已有的古生物仿真研究多将仿真域中模拟海床的底面设为平滑的表面,这种设置方式可在理想情况下研究底栖型生物在不同流速环境中的受力与摄食效率等问题。而本文的研究作为底栖微体化石流体仿真的延伸工作,加入了随机表面生成的方法,首次实现了将底栖微体化石与模拟起伏海床表面环境的仿真域相结合,进一步还原了底栖微体生物原本的生活环境及所处的水动力状况。由于底栖微体生物的生态特征较为复杂,根据其在沉积物中的生活位置可以分为表栖类型和内栖类型,因而不同类型底栖微体生物生存空间有所差别,除在沉积底质表面上以外,还可以出现在表层之下的不同深度,而本文采用的CFD方法尚难以精确还原内栖生物的生活环境。此外,受频繁多变的海洋环境以及实际海床表面沉积质结构的影响,海底边界层中的流动应当远比本文实验中模拟的状态更为复杂(如沉积底质表面的粗糙度以及其中的孔隙可能会对海底边界层的流动状态存在进一步的微观影响),但本文的仿真仍在毫米至厘米级别的尺度上尽可能地还原了海底边界层区域中的流动环境,为古生物与流体仿真相结合的研究带来了新的进展。本文的工作不仅提供了仿真数据的参考,对于古生物领域研究者更深入地认识古生物与其生活的海洋、地质环境之间的联系亦有着重要意义。

衷心感谢王小龙老师对本文审稿工作的大力支持以及中国地质大学(北京)万晓樵教授提出的宝贵审稿意见!

参考文献

[1] BORAZJANI I, SOTIROPOULOS F. Numerical investigation of the hydrodynamics of carangiform swimming in the transitional and inertial flow regimes[J]. Journal of Experimental Biology, 2008, 211(10): 1541-1558.

[2] CURATOLO M, TERESI L. The virtual aquarium: simulations of fish swimming[C]// Proceedings of European COMSOL conference. Stockholm: COMSOL, 2015.

[3] GEMMELL B J, COSTELLO J H, COLIN S P, et al. Passive energy recapture in jellyfish contributes to propulsive advantage over other metazoans[J]. Proceedings of the National Academy of Sciences of the United States of America, 2013, 110(44): 17904-17909.

[4] GEMMELL B J, COLIN S P, COSTELLO J H. Widespread utilization of passive energy recapture in swimming medusae[J]. Journal of Experimental Biology, 2018, 221(1): jeb.168575.

[5] MCHENRY M J, JED J. The ontogenetic scaling of hydrodynamics and swimming performance in jellyfish (Aurelia aurita)[J]. The Journal of Experimental Biology, 2003, 206(22): 4125-4137.

[6] WATERS J A, WHITE L E, SUMRALL C D, et al. A new model of respiration in blastoid (Echinodermata) hydrospires based on computational fluid dynamic simulations of virtual 3D models[J]. Journal of Paleontology, 2017, 91(4): 662-671.

[7] RAHMAN I A, O’SHEA J, LAUTENSCHLAGER S, et al. Potential evolutionary trade-off between feeding and stability in Cambrian cinctan echinoderms[J]. Palaeontology, 2020, 63(5): 689-701.

[8] SONG H C, SONG H J, RAHMAN I A, et al. Computational fluid dynamics confirms drag reduction associated with trilobite queuing behaviour[J]. Palaeontology, 2021, 64(5): 597-608.

[9] LABARBERA M. Feeding currents and particle capture mechanisms in suspension feeding animals[J]. Integrative and Comparative Biology, 1984, 24(1): 71-84.

[10] 许淑梅, 孔家豪, 李三忠, 洋板块地层学的概念、模式、组成及失序变化[J]. 古地理学报, 2022, 24(1): 1-10.

[11] TROWBRIDGE J H, LENTZ S J. The bottom boundary layer[J]. Annual Review of Marine Science, 2018, 10: 397-420.

[12] CALDWELL D R, CHRISS T M. The viscous sublayer at the sea floor[J]. Science, 1979, 205(4411): 1131-1132.

[13] WENGROVE M E, FOSTER D L. Field evidence of the viscous sublayer in a tidally forced developing boundary layer[J]. Geophysical Research Letters, 2014, 41(14): 5084-5090.

[14] GRANT W D, MADSEN O S. The continental-shelf bottom boundary layer[J]. Annual Review of Fluid Mechanics, 1986, 18: 265-305.

[15] LIU P, ZHANG Y H, YANG X G, et al. Hydrodynamic simulations of millimeter-scale Cambrian sedmentary medusozoans[J]. Journal of Geophysical Research: Biogeosciences, 2022, 127(10): e2022JG006854.

[16] RAHMAN I A. Computational fluid dynamics as a tool for testing functional and ecological hypotheses in fossil taxa[J]. Palaeontology, 2017, 60(4): 451-459.

[17] CHRISS T M, CALDWELL D R. Universal similarity and the thickness of the viscous sublayer at the ocean floor[J]. Journal of Geophysical Research: Oceans, 1984, 89(C4): 6403-6414.

[18] PRATT M C. Living where the flow is right: how flow affects feeding in bryozoans[J]. Integrative and Comparative Biology, 2008, 48(6): 808-822.

[19] CSANADY G T. Circulation in the coastal ocean[M]// Advances in geophysics: Vol. 23. Amsterdam: Elsevier, 1981: 101-183.

[20] 周济福, 谢益芹. 近海岸波流边界层及泥沙输运研究[C]// 第十六届全国水动力学学术会议暨第三十二届全国水动力学研讨会论文集. 上海: 中船编印社有限公司, 2021: 7.

[21] GIBSON B M, RAHMAN I A, MALONEY K M, et al. Gregarious suspension feeding in a modular Ediacaran organism[J]. Science Advances, 2019, 5(6): eaaw0260.

[22] DENNY M. Biology and the mechanics of the wave-swept environment[M]// Biology and the mechanics of the wave-swept environment. Princeton: Princeton University Press, 2014: 117-133.

[23] 刘平, 杨晓光, 唐清钦, 寒武纪早期锥管状化石Anabarites trisulcatus群居行为的计算流体力学分析[J]. 微体古生物学报, 2021, 38(3): 234-240.

基金资助

中国科学院战略性先导科技专项(B类)(XDB26000000)

国家自然科学基金项目(41672009)

国家自然科学基金项目(41621003)

国家自然科学基金项目(41720104002)

国家自然科学基金项目(42172016)

国家重点研发计划项目(2023YFF0803601)