引言

在工程应用中，流体和固体结构的相互作用是一个普遍但又复杂的问题.其中，涡激振动（vortex-induced vibration，VIV）与涡激旋转（vortex-induced rotation，VIR）十分重要.涡激振动关注固定物体的弹性振动，考虑弹性形变，例如油气管道在深海中的振动，高压电线、桥梁在风中的振动等^[1，2，3].涡激旋转关注刚体在流体中的旋转响应，一般不考虑物体本身的形变.空气中旋转落地的水果或冰雹、流体中自由运动的物体等均可以简化为涡激旋转；陈卓等^[4]的研究还发现，在流场中自由摆动方柱的涡激旋转能强化微流体通道内热量与物质的交换.对涡激旋转的研究能够更加有利于解释如从飞行器中分离的物品、海洋中旋转能量收集装置桨片^[5]等物体的旋转特性，促进对旋转特性的控制.

针对涡激旋转的研究目前主要集中于实验与数值仿真.早期研究中，研究者主要聚焦于刚体涡激旋转表现出的自旋转（autorotation）现象，最早由Maxwell和Niven^[6]展开.Lugt^[7]通过对Navier-Stokes 方程的数值求解，研究了平板与椭圆柱的自旋转现象，指出自转的主要原因是漩涡脱落与刚体旋转的同步化，因此相关现象也被称为涡激旋转.Skews^[8]通过实验发现，多边形柱的涡激旋转主要发生在边数小于或等于8的情况.特别地，针对方柱的研究，Zaki^[9]、Ryu与Iaccarino ^[10]通过数值仿真观察到了在不同雷诺数下方柱运动的模态.Ryu与Iaccarino^[10]指出，在小雷诺数（Re=45~150）情况下，方柱的涡激旋转存在六个基本模态：静止模态（initial rest position，Re=45），小幅振荡模态（small-amplitude oscillation regime，Re=50，60，70），π/2极限振荡模态（π/2-limit oscillation regime，Re=80，90），随机旋转模态（random rotation regime，Re=100，110），π极限振荡模态（π-limit oscillation regime，Re=120，130），自转模态（autorotation regime，Re=140，150）.[注：为了进行区分，本文中涡激振动对应的弹性振动（vibration）翻译为振动，而涡激旋转模态中在平衡位置周围的抖动（oscillation）翻译为振荡.]这些数值仿真结果，为理论分析方柱涡激旋转现象的原因提供了分类依据，也能指导在工业中涡激旋转效应在能量收集、姿态控制等方面的应用.

基于以上模态分类，少量特殊情况下方柱涡激旋转的理论解释已经提出，但集中于稳定性与驱动力的定性分析.Robertson等人^[11]针对长方形方柱的涡激旋转，通过刚体运动准定常理论（quasi-steady theory）的方法，提出了旋转稳定性的判断条件. Ryu与Iaccarino^[10]早期的研究指出，小幅振荡模态中周期性的振荡主要是由所受力矩与旋转角度之间不同步而驱动的.后来，针对方柱涡激旋转在Re=150下的自转模态，Ryu^[12]通过对于迎风面、逆风面分类研究的方法，计算了定常流场情况下的压强导致的力矩，指出大角度旋转带来的“瞬时”力矩突然换向是维持自转模态旋转的主要原因.同时，在这一模态下，迎流面对力矩的贡献与瞬时总力矩方向相反，阻碍方柱的转动；背流面对力矩的贡献最大且与瞬时力矩同向，促进方柱转动；总力矩方向与旋转方向基本相同.这些研究均只深入到对力矩这一宏观量的定性研究，而没有对力矩来源——流场状态进行解释，难以刻画涡激旋转中的流体力学现象，从机理上解释造成力矩变化的原因.

本研究的主要目的是针对方柱小幅振荡模态，以流场分析为主要方法提出一理论模型，作为通过理论定量解释这一振荡现象的初步尝试.这一模型充分结合钝体绕流中已有的理论模型，考虑方柱带来的“扰流效应”与黏性流体带来的“边界层效应”，结合自由流线理论与边界层理论进行计算.我们将该理论模型称为自由流线-边界层理论模型（Free Streamline-Boundary Layer Model，缩写为FSBL Model）.

本文的组织结构如下：第一节主要介绍模型的基本假设.第二节对远离方柱、黏性体现不明显的流场依据自由流线模型进行计算.第三节对靠近方柱、充分体现流体黏性作用的流场，依据边界层模型进行计算.第四节结合浸没边界法的仿真数据，进行假设验证、模型计算与结果分析.最后在第五节对研究结果进行总结.

1 模型描述及基本假设

本文研究的方柱绕流为图1所示的二维问题.均匀来流速度为U_∞，朝x正方向.方柱边长为L，中心固定，且初始时刻一面正对来流.方柱旋转后，取逆时针方向旋转角度为α.由于研究只考虑方柱小幅振荡的情况，故始终有$\left|\alpha \right|\ll \pi /4$; 且方柱正对来流区为同一表面，记为Σ₁.其余各面沿顺时针方向依次记为Σ₂~Σ₄.为了方便后续的研究，沿Σ₂与Σ₄各建立一辅助右手坐标系（x₊，y₊），（x_-，y_-）.

本文研究中将采用如下假设：

（1）假设流体为不可压缩流体.

（2）假设流体的非定常性对于方柱的力矩及一时刻方柱所受力矩时，可以假设流场为定常流场.在此定常假设下，辅助正交坐标系（x₊，y₊），（x_-，y_-）为惯性系.该定常假设的有效性将会在第4小节进行验证与说明.

（3）方柱的小幅振荡可以近似认为是简谐振动，故其所受力矩M近似满足M∝α，记该比例系数为M_α=-M/α.由此，方柱的转动满足方程：

其中，I=1/6ρ_sL⁴为方柱的转动惯量，ρ_s为固体的密度.由式（1）可得，方柱的振荡频率为$f=1/2\pi \sqrt{M_{\alpha }/I}$.该假设的有效性将会根据第4节数值结果进行验证.值得注意的是，简谐振动的假设也与Ryu^[10]指出的小幅振荡驱动力是所受力矩与旋转角度之间异步性的结论一致.

（4）假设远离方柱的流场变化主要是由方柱绕流造成的，流体的黏性可以忽略不计.并且，由于正对流体的仅为Σ₁面，我们可以认为绕流效应主要是由Σ₁面造成的.

（5）假设方柱绕流区内雷诺数较小，没有因为绕流出现的分离泡，在方柱表面主要表现为边界层效应，该假设称为“边界层假设”.该假设的有效性将会根据第4节数值结果中Σ₂，Σ₄周围流场速度与涡量分布进行验证.

（6）假设固体Σ₃面各点所受压强呈线性变化，黏性力可以忽略不计.该假设的有效性将会在第4节通过数值结果进行验证.需要说明:本文考虑的是α很小的运动模态，而Ryu^[12]的研究表明，在α较大时该假设不再有效.

2 自由流线模型及其理论分析

2.1 自由流线物理模型的基本假设

1869年，Kirchhoff^[13，14]针对圆柱绕流的计算提出了自由流线的概念.自由流线的基本假设是钝体绕流的分离区内部为“死水区”，“死水区”内部流体流动基本均匀且可以忽略不计，压力为常数.“死水区”由两个被称为“自由流线”的间断面围成.“自由流线”自与钝体接触起一直沿着钝体表面直至分离点，在分离点与钝体表面相切，满足流函数为0. 其中，Kirchhoff ^[13，14]与Brodetsky^[14，15]假设自由流线上速度为常数并始终等于来流的速度，针对圆柱绕流计算所得分离角θ_s与压差阻力系数c_Dp同实验所得有较大差异.而 Roshko^[14，16]，Eppler^[14，17]，Woods^[14，18]针对圆柱的研究注意到，假设钝体后方的压力为一合适的负值，并且允许速度沿自由流线进行变化，所得分离角θ_s与压差阻力系数c_Dp更贴近实验所得，能更好模拟出无黏绕流流场的压强分布.以Roshko^[14，16]的模型为例，对关于来流对称性的钝体（如圆柱），该模型假设自由流线在分离点具有相同的速度kU_∞.分离后，速度沿自由流线先匀速变化至与来流平行（“分离区”），再减速至无穷远速度U_∞（“减速区”）.

自由流线参数方程及速度的求解主要通过物理平面（z）、复势平面（F）、归一化复速度倒数平面[ζ=U/W=U（dz/dF），U为归一化速度]与实轴上半空间$\{t\in \mathbb{C}\mid \mathcal{R}（t）⩾0\}$的保角变换来进行，$\mathcal{R}（\cdot ）$表示取实部.在计算正对来流平板时，在Roshko的模型中，分离区在归一化复速度倒数平面（ζ）中的轨迹为半圆.

依据第1章的假设，远离方柱的流体黏性可以忽略不计，因此可以认为“自由流线”区域满足无黏的要求.假设“自由流线”与斜板分离后，在改变方向的同时减速，且在归一化复速度倒数平面（ζ）中的轨迹为一倾斜的椭圆，椭圆两长轴端点为分离点，为了便于后期计算，假设椭圆半短轴长为1. 由于归一化参数U的引入，这一关于半短轴的假设是不失一般性的.这一假设相对于Roshko的假设，其速度变化更加平滑.

2.2 自由流线轨迹及速度理论模型

自由流线在物理平面（z）、复势平面（F）、归一化复速度倒数平面（ζ）中的轨迹如图2所示，为了方便将其保角变换至实轴上半空间，我们还设置一辅助平面（$\tilde{\zeta }$），将椭圆保角变换为圆.

首先对归一化复速度倒数平面（ζ）中轨迹进行详细说明: 由于自由流线与障碍物交点处（C）速度为0，C在ζ平面上表现为无穷远.而无穷远处（I）速度为U_∞=k_∞U，分离点处（S₊，S_-）速度大小分别为k₊U，k_-U.因此在ζ平面上I坐标为1/k_∞，分离点的坐标分别为S₊=（1/k₊）ie^iα，S_-=-（1/k_-）ie^iα. 由于椭圆半短轴长度及长轴方向已经给定，实际椭圆只有2个自由度.考虑到物理意义，研究取k₊/k_∞，k_-/k_∞作为自由参数，用于表示两分离点与无穷远处速度的比值.该自由参数可以通过数值仿真数据或实验测得.

各个复平面之间的变换由以下各式给出：

其中k_∞，λ，β，m，χ，C为曲线参数，综合式（4）、式（5），可以化简得到

式中正、负号分别代表着上方与下方自由流线.综合式（1）、式（2）、式（5），可得物理空间自由流线轨迹与复势的关系：

式中z₀为零势能点在z平面中的坐标.复速度与复势的关系：

根据自由流线定义，自由流线应取$F\in {\mathbb{R}}^{+}$.

下面对曲线参数（k_∞，λ，β，m，χ，C）进行计算.在计算时，首先假设k₊，k_-已知，求出其他量，特别是k_∞与k₊，k_-的关系，此后可通过二分法，在给定k₊/k_∞，k_-/k_∞时，确定k₊，k_-，k_∞的值.

参数k_∞，λ，β，m，χ的确定依赖于归一化复速度倒数平面（ζ）与圆平面在端点与零点处的对应关系.为了方便，记$\frac{1}{k}=\frac{1/k_{+}+1/k_{-}}{2}$.可得：

其中，式（12）为χ满足的方程，其解为式（14）中的一个，正负号的确定根据其是否满足式（12）.

参数C的确定可以依赖二分法，其确定标准应当为使得两分离点间的距离为斜板的长度，其中分离点处复势的确定方式为

式中$\mathcal{R}[\cdot ]$表示取实部.由于C是函数g中的参数，进而也是z（F）的参数，在这里将其记为z（F; C）.C的取值应满足：

式中I{ · }表示取虚部.由此，关于自由流线轨迹及其速度的具体求解步骤为：

（1）根据数值仿真或实验结果测量得到k₊/k_∞，k_-/k_∞，进而由式（9）-式（11）、式（13）-式（14）递归确定参数k_∞，λ，β，m，χ.

（2）结合式（3）、式（6）、式（7）、式（15）、式（16），使用二分法递归确定参数.

（3）根据式（6）-式（7）确定自由流线轨迹及其速度.为了方便后续推导，令Σ₁表面的速度模长为$W_1^{FS}$（x₁，y₁）.将所求得的上下自由流线轨迹分别以图1所示辅助坐标系（x_±，y_±）表示.令上下自由流线（“分离区”）在z_±坐标系中轨迹为y_±=$y_{\pm }^{FS}$（x_±），速度在（x_±，y_±）坐标系中方向分别为$u_{\pm }^{FS}$（x_±），$v_{\pm }^{FS}$（x_±），记其模长为$W_{\pm }^{FS}\left(x_{\pm }\right)=\sqrt{u_{\pm }^{FS}{\left(x_{\pm }\right)}^2+v_{\pm }^{FS}{\left(x_{\pm }\right)}^2}$.

3 边界层模型及其理论分析

对于靠近Σ₂，Σ₄面的流场，由于流体的黏性作用，Σ₂，Σ₄表面会形成一边界层，可借助边界层理论对该表面流场进行计算.1908年，Blasius^[19]通过数量级分析将Navier-Stokes方程简化，利用自相似原理求得了平板下的流场速度.Falkneb与Skan^[20]则在就楔形劈尖的边界层求解时，通过对远离边界的流场进行Bernoulli假设，求解了存在压强梯度的边界层方程.在本文研究模型中，以自由流线速度为足够远处速度，通过数量级分析将Navier-Stokes方程简化，利用自相似性求得流场速度.

3.1 边界层方程的简化

为了方便表述，本章的研究采用z_±坐标系，记s=±，记速度沿x_s，y_s方向的分量分别为u_s，v_s.

首先分析二维黏性流体的方程，有:

以上三个方程中，v为运动黏度，ρ_f为流体密度.由于定常假设，时间偏导项$\partial u_s/\partial t$与$\partial v_s/\partial t$为0，在下方分析中略去.下面对以上三个方程进行数量级分析.

假设边界层“足够远”处速度量级为U，V，特征板长度为l，该坐标的边界层高度为δ，且$\delta \ll l$，那么由式（17），可得：

$V\sim \frac{U}{l}\delta$

对于式（18），各项数量级依次为

$\frac{U^2}{l};\frac{VU}{\delta }=\frac{U^2}{l};\frac{1}{{\rho }_f}\frac{\delta p}{\delta x};\frac{\nu U}{l^2};\frac{\nu U}{{\delta }^2}$

考虑到等式左边数量级为U²/l，等式右边由于$\delta \ll l$，νU/l²可以省略，则

$\frac{U^2}{l}\sim \frac{1}{{\rho }_f}\frac{\delta p}{\delta x}\sim \frac{\nu U}{{\delta }^2}$

可得

$\delta \sim \sqrt{\frac{\nu l}{U}}=\frac{l}{\sqrt{{\mathrm{R}\mathrm{e}}_l}}$

在本文模型中，取足够远速度，即为$u_{\pm }^{FS}$，为了方便后面的计算，定义边界层厚度

${\delta }_s\left(x_s\right)=\sqrt{\frac{\nu x_s}{u_{fs}^s}}=\frac{x_s}{\sqrt{{\mathrm{R}\mathrm{e}}_{x_s,s}}},s=\pm$

对于式（18），同样进行数量级分析：

$\begin{array}{c}\frac{UV}{l}=\frac{U^2\delta }{l^2};\frac{V^2}{\delta }=\frac{U^2\delta }{l};\frac{1}{{\rho }_f}\frac{\delta p}{\delta y}\\ \frac{\nu V}{l^2}=\frac{U^2{\delta }^3}{l^4};\frac{\nu V}{{\delta }^2}=\frac{U^2\delta }{l^2}\end{array}$

可得

$\frac{1}{{\rho }_f}\frac{\delta p}{\delta y}\sim \frac{U^2\delta }{l^2}$

则$\frac{1}{{\rho }_f}\frac{\delta p}{\delta y}\ll \frac{1}{{\rho }_f}\frac{\delta p}{\delta x}$，压强在y方向的偏导可忽略不计，即

同时，对于边界层的压强p_s，有

由于忽略了流场黏性，有Bernoulli方程：

综合式（18）和式（22），有

记

$\begin{array}{c}{A}_s\left(x_s\right)=\frac{W_s^{FS}\left(x_s\right)}{{\left[u_s^{FS}\left(x_s\right)\right]}^2/x_s}\frac{\mathrm{d}{W}_s^{FS}\left(x_s\right)}{\mathrm{d}{x}_s}\\ B_s\left(x_s\right)=\frac{x_s}{u_s^{FS}\left(x_s\right)}\frac{\mathrm{d}{u}_s^{FS}\left(x_s\right)}{\mathrm{d}{x}_s}\end{array}$

在数值计算中，A_s，B_s参数中微分的计算可以通过差分实现.式（23）可简化为：

同时有边界条件：

3.2 根据自相似性化简边界层方程

考虑边界层的流函数，有

$u_s=\frac{\partial {\psi }_s}{\partial y_s},v_s=-\frac{\partial {\psi }_s}{\partial x_s}$

则式（24）可写为

边界条件式（25）-式（27）可化简为

不失一般性地假设ψ_s（x_s=0，y_s=0）=0，那么由式（30）可得

${\psi }_s\left(x_s,y_s=0\right)=0$

对于实际速度，我们尝试求解在x_s方向上具有形式自相似性的速度，即存在函数Φ_s使：

对式（31）在y_s方向积分，可以得到：

${\psi }_s\left(x_s,y_s\right)=f_s\left(\frac{y_s}{{\delta }_s\left(x_s\right)}\right){\delta }_s\left(x_s\right)u_s^{FS}\left(x_s\right)$

式中，$f_s^{\mathrm{\text{'}}}\left({\eta }_s\right)={\mathrm{\Phi }}_s\left({\eta }_s\right)$，记${\eta }_s=\frac{y_s}{{\delta }_s\left(x_s\right)}$.将式（32）代入式（28），可得

同时由式（25）-式（27）可得边界条件：

计算中需要在每一个位置x_s对式（33）、式（34）求解.求解时可以采用试射法，基本思路是假设f″_s（0）=a为一给定数，数值求解方程并计算$f_s^{\mathrm{\text{'}}}\left[\frac{y_s^{FS}\left(x_s\right)}{{\delta }_s\left(x_s\right)}\right]$，通过递归使$f_s^{\mathrm{\text{'}}}\left[\frac{y_s^{FS}\left(x_s\right)}{{\delta }_s\left(x_s\right)}\right]$为1.

3.3 力矩的计算及基本假设

求得边界层速度分布后，可根据定常假设计算方柱所受力矩M.该力矩由压强力矩M_p与黏性力矩M_μ组成.

在压强力矩M_p的计算中，首先分析Σ₂与Σ₄面所受力矩，其基本思路是由边界层公式进行计算.由式（21）、式（22）得：

其中，n=2，4并且分别对应s=+，-.实际计算中取p_∞.对于Σ₁面，由于自由流线流过方柱表面，所以

对于Σ₃面，可以对其表面压强作线性拟合.有力矩

对于黏性力矩的求解，只考虑Σ₂与Σ₄面，有：

式中μ为黏度，积分可得力矩，有

式（38）、式（39）中，n=2，4分别与s=+，-相对应.

4 数值算例

本章展示在针对小幅振荡的自由流线-边界层理论模型与数值计算的结果对比.

4.1 数值分析方法

研究采用浸没边界法对均匀来流中方形柱体涡激旋转的动力学问题进行数值模拟，分别对流体和刚性固体的控制方程进行数值求解.在控制方程中，浸没边界法使用附加力项来表示流体和固体之间的相互作用.数值计算方法详情参见文献^[21，22].

在流体运动上，研究计算域长为32L，宽为16L，且设置计算域边界速度均为U_∞沿x正方向.计算时，流体连续性方程与Navier-Stokes方程均需进行无量纲化，其对应特征量如表1所示，无量纲化后所得Navier-Stokes方程如式（40）所示.式中，雷诺数Re=ρ_fU_∞L/μ，F_ib表示浸没边界法所施加的力.

在下文进行的所有计算中，本文第2、3章所列流场方程与所有流场物理量也均采用表1特征量进行无量纲化.为简化，所有无量纲化之后的物理量均采用原记号标记.

所得连续方程与Navier-Stokes方程可以由欧拉网格形式进行划分，依据Peskin^[23]所提出的方法，将流体速度和力源项定义在网格节点上，将流体压强定义在网格中心.时间离散采用全隐式时间推进法，空间离散采用Crank Nicholson 格式来离散对流项和扩散项.

对于固体部分，浸没边界法只关注固体边界，即浸没边界的位移与速度.力矩、角加速度、角速度和角度的关系同样进行无量纲化，其中转动惯量特征量为ρ_sL⁴，时间特征量为L/U_∞，角度特征量为1.在下文进行的刚体运动计算中，本文第2、3节所列流场方程与所有流场物理量也均采用这些特征量进行无量纲化.为了简化，所有无量纲化之后的物理量均采用原记号标记.离散方法上，固体被视作刚体，采用显式欧拉法离散刚体的角度、角速度、角加速度之间的关系式.

根据浸没边界法，式（38）中流体与固体的相互作用力可以表示为：

式中，α、β为常系数，U_ib（s，t）是由界面附近流体速度插值得到的浸没边界处的流体速度，$\tilde{U}$（s，t）为浸没边界的运动速度.其中，浸没边界处的流体速度可以由界面附近流场速度插值得到：

式中，V表示边界附近流体所占的区域，u（x，t）表示流体速度，X（s，t）表示固体边界节点的位置，x表示流体欧拉网格的位置坐标，δ是Dirac Delta 函数.

根据以上方法，针对不同ρ_s/ρ_f与Re进行了研究.我们发现，在雷诺数Re=70与密度比ρ_s/ρ_f=8时存在明显且稳定的小幅振荡情况.下面的章节中，我们将以该情况为例对理论模型进行验证.

4.2 假设有效性检验

浸没边界法所得t=0~600时间区间内M-t与α-t曲线如图3所示.

下面首先分析各假设的有效性.对于定常假设，我们对边界层进行数量级分析，式（18）中，令δU为U变化量.公式的数量级为（δU）²/t，而时间偏导项$\partial u_s/\partial t$数量级为fδU，由数值仿真数据可得，$f\sim 0.1\ll 2U/L=2\sim \delta U/l$，故可以将非定常项忽略不计，定常假设成立.

对于简谐振动假设，下面考察小幅振荡期间力矩M与旋转角度α的线性.t=550~600时间窗口内的M-α散点图如图4所示.对M与α数据进行线性回归可得R²=1.0000，其线性极强，简谐振动假设成立.

由于定常假设，可以取任意时刻的方柱旋转角度α，将流场视为定常进行计算.下文以t=600为例研究流场情况.

首先验证方柱表面流场情况，分析方柱绕流区是否存在局部分离泡.图5展示了方柱绕流区流线与涡量大小|w|分布情况.由图5可见，绕流区内流线主要呈现为层流特征，涡量主要集中于方柱拐角处，没有明显独立小涡出现.图6展示了Σ₂与Σ₄表面邻域的速度剖面图.由图6可见，两表面的速度分布呈现明显的边界层特征，没有局部涡的形成，边界层假设成立.

对于Σ₃面的线性压强假设，将浸没边界法所得压强结果以压差系数$c_p=\left(p-p_{\mathrm{\infty }}\right)/\frac{1}{2}\rho /U_{\mathrm{\infty }}^2$表示.取t=600时刻，各个面的c_p随该面距离坐标变化如图7所示，由图7可以看出，除端点处外，Σ₃面的c_p曲线较其他面c_p变化具有较强的线性，且端点处变化较其他曲线变化可以忽略不计，故线性假设基本成立.

4.3 模型与数值结果对比

对于自由流线-边界层模型，根据第2、3章的方程编制MATLAB程序进行计算.由浸没边界法得到t=600时刻数值仿真数据，α=-0.00329215，由上、下方角落邻域内压强计算得k₊/k_∞=1.346886，k_-/k_∞=1.387032.带入MATLAB程序中，所得力矩结果如表2所示. 由浸没边界法得到该时刻力矩为M=0.003921，M_α=1.1910.在t=550~600区段振荡频率为f=0.1332.

可以看出，两方法所得M_α与f基本相当.从表2可以得出，维持方柱小幅度振荡的力主要来源于Σ₁面的压强力矩及Σ₂与Σ₄面的压强力矩差，即迎流面（Σ₁）压强的分布不均匀与上下面（Σ₂与Σ₄）因转动与自由流线速度差异导致的压强不对称性.值得注意的是，与自转模态不同，小幅振荡模态下迎流面（Σ₁）的力矩贡献与瞬时力矩方向相同，并且为主要组成力矩之一；而背流面（Σ₃）对于力矩的贡献与瞬时力矩方向相反，并且可以忽略不计.这也使得此时的力矩M主要表现为恢复力矩，其方向与α方向相反，方柱呈现振荡形态.

下面比较模型所得的压强结果与浸没边界法所得Σ₂与Σ₄面的压强结果，以分析可能的误差原因.结果以上方所定义压差系数c_p表示.

由图8可见，计算所得压差系数总体偏大，趋势与数值仿真结果基本相仿.其原因一方面是浸没边界法所得压强，并非边界处压强，而为靠近边界点的压强此处无法进行完全比较.另一方面，所取k₊/k_∞与k_-/k_∞将会极大影响c_p计算精度，而浸没边界法所得数据可能存在偏差.同时，定常假设导致压强中缺乏非定常变化项以及对迎流面黏性力矩无法测算，使得压强计算存在系统性偏差.

5 结论

本文综合钝体绕流的自由流线与边界层理论，提出了流固耦合中方形柱体涡激旋转小幅振荡模态下的自由流线-边界层理论解释模型，并且通过浸没边界法进行数值仿真，对比验证了该理论的有效性.研究总结如下：

（1）从数值仿真结果可以得出，方柱在小幅振荡模态中可以近似为简谐振动，并且由自由流线-边界层理论推算所得力矩和频率与数值仿真基本吻合.

（2）从力矩来源上讲，方柱涡激振荡维持力矩主要来源于迎流面压强的分布的不均匀以及上下两面压强分布的不对称性.

（3）小幅振荡模态与自转模态下压强的分布不同，所导致模态的主要驱动力来源也不同，这也是这一模态下表现为振荡的主要原因.

尽管在压强的计算上，由于忽略了非定常性与迎流面与背流面表面黏性，自由流线-边界层理论相较于数值仿真存在一定的差异，但是本研究将研究点聚焦于靠近方柱周围流场的局部化描述，给涡激旋转与部分流固耦合问题的理论解释提出了一个新思路.

参考文献