【《流固耦合数值计算方法综述》12000字】

--流固耦合数值计算方法综述目录TOC\o"1-3"\h\u22142流固耦合数值计算方法综述 1179171.1流固耦合计算流程 1111791.2计算流体动力学基本理论 5236351.1.1流体动力学控制方程 5273151.1.2湍流模型 6200261.3计算结构动力学基本理论 10267831.4计算域离散与动网格技术 12142601.4.1计算区域的离散 12266501.4.2动网格方法 13216711.5数据交换方法 14209161.6方柱下游柔性板流致振动研究 1578391.6.1计算模型 1518201.6.2计算域离散化及边界条件 16198621.6.3结果分析 18从工程应用角度来看，流体力学和固体力学的发展是为了解决工程中遇到的实际问题。流固耦合研究是目前工程应用上的热点问题，也是实际工程设计中对计算流体力学和计算固体力学的有效应用。随着计算资源的不断增加，大量理论工作进入了各种流固耦合应用的求解框架。本章介绍了研究所采用的流固耦合数值计算方法，及其实现基础和理论，并且模拟了方柱下游柔性板流致振动现象，用于验证该方法的有效性。1.1流固耦合计算流程流固耦合问题根据耦合形式的不同可以分为两类。第一类问题的特征是流体区域和固体区域全部重叠或部分重叠在一起，难以明显地分开，需要针对具体的物理现象建立本构方程，如渗流问题。第二类问题的特征是耦合作用仅仅发生在两相交界表面处，在方程上的耦合是通过交界处边界条件的平衡和协调关系引入的，如图1.1。整个仿真区域Ω包含两个子区域，即流场ΩF及其边界ΓF和固场ΩS及其边界ΓS。显然，Ω=ΩF+ΩS。子区域之间存在一个公共边界Γ，定义为Γ=ΩF∩ΩS，并且ΓΓF，ΓΓS，称为流固耦合交界面。本文研究的水平轴潮流能水轮机流固耦合问题属于第二类。图1.1表面流固耦合问题几何区域定义Figure1.1Definitionofgeometricregionforsurfacefluid–solidinteractionproblem求解这第二类流固耦合问题的数值方法可大致归为两种：整体求解法ADDINNE.Ref.{54C15CA1-61D8-491E-39B}[129-132]和分区求解法ADDINNE.Ref.{2EAD7603-D24C-4BB4-A318-5E4C7DFCECAF}[133-135]。整体求解法是通过改写流体方程和固体方程，将其放在同一方程组中，采用统一算法同时进行求解，流体和固体交界面的状态信息隐含在求解过程中ADDINNE.Ref.{C869C7D3-774B-48F7-B1B8-73108C0765BF}[136]。这种方法对于多学科交叉问题能够取得较好的计算精度，但是需要更多的计算资源，适用于求解梁、板、壳等一些简单结构ADDINNE.Ref.{D2B253D4-7AAE-4FA4-BC79-39E3AAD8B7F7}[137-140]。分区求解法是将流体方程和固体方程交替进行求解，每一个物理场计算得到的交界面信息代入另一个物理场继续计算ADDINNE.Ref.{3A39C91E-B6C5-4A17-B505-5D6A951A089C}[141,142]。与整体求解法相比，分区求解法中物理场计算之间相对独立，每一时刻计算方程数量少，耦合作用体现于每次迭代时交界面的信息传递ADDINNE.Ref.{4DA5E08A-64BF-4C86-92E7-34A0923395DD}[143,144]。分区求解法的优势在于既能考虑到流场和固场离散化方法的差异，又能充分利用现有成熟的计算流体力学和计算固体力学程序，灵活选择不同的求解器，因此广泛应用于大型复杂的工程问题ADDINNE.Ref.{985D6F7C-FA88-4470-B6BF-DB9C4ECCAE01}[145,146]。潮流能水轮机结构复杂，所用网格数量多，模拟时间长，因此选择采用分区求解法。从数据传递方向出发，分区求解法进一步可划分为单向流固耦合和双向流固耦合。单向流固耦合是指耦合交界面上的数据传递是单向的，通常将流场计算得到的交界面荷载传递到固体结构分析，但没有固场计算结果返回给流场的过程。也就是说，只考虑流体分析对固体分析的影响，而结构变形对流体影响小至可以忽略不计。例如，旋转机械的结构强度分析。另外，研究已知运动轨迹的刚体对流场的影响，在某种程度上也可以看作是一种单向流固耦合问题。例如快速启闭阀门对流体运动的瞬时影响分析。由于固体的运动已知，且看作为刚体，变形忽略不计，所以此类问题通常可以通过自定义函数（UserDefinedFunctions，简称UDF）的方式单独在CFD求解器中完成。双向流固耦合是指耦合交界面上的数据传递是双向的，也就是，每个物理场在每次迭代计算初都接收并更新另一物理场在上次迭代完成后传递而来的交界面物理量，并且在迭代收敛后将新的物理量传递给另一物理场用作下次迭代。由于求解有限体积法离散方程和有限元法离散方程所采用的时间积分不同，为了避免流场计算和固场计算存在时间上的滞后，提高数值计算的精度和稳定性，在一个时间步内，耦合交界面需要进行多次数据交换，各物理场进行多次迭代计算，直到流场、固场和交界面信息均满足精度要求后，才可以同时进入下一个时间步的计算。对大多数耦合作用现象，如果只考虑静态结构性能，采用单向耦合分析足以。但如果要考虑振动等动力学特性，双向耦合分析必不可少。因为，双向耦合分析主要是为了解决振动和大变形问题而进行的。例如潮流能水轮机在流场中的静态结构分析属于单向耦合分析，而要是考虑叶片在流场中的振动情况，就需要用到双向耦合方法进行分析。图1.2给出了本文研究所采用的双向流固耦合计算基本流程。在每次耦合子迭代循环中先计算流场再计算固场。单向流固耦合计算时，如果只考虑流体对固体的影响，则取消更新耦合边界位移这一步骤；如果只考虑固体对流体的影响，可以取消更新耦合边界荷载这一步骤，也可以如之前所述，单独在CFD求解器内完成计算。图1.2双向流固耦合计算基本流程Figure1.2Basicflowdiagramoftwo-wayfluid–solidinteractioncalculation从计算流程图可以看出，本文的流固耦合方法是一个典型的Dirichlet–Neumann分区求解方法。其特征为流体方程是在流固耦合交界面的位移确定后再求解，而固体方程是在交界面荷载确定后再求解。例如，在第n个时间步的第1个耦合子迭代循环开始时，流场依据上一个时间步得到的耦合交界面位移，dnΓF,1=df（1.1）在固场更新了耦合交界面荷载分布之后，fnΓS,1=fd（1.2）此时，第n个时间步的第1次迭代完成，需要判断耦合交界面处的收敛情况。式（1.2）表明，该方法可以将交界面位移dΓ视为流体和固体相互作用的主要变量，而界面荷载r（1.3）式中：r为位移残差；εFSI1.2计算流体动力学基本理论流固耦合数值分析作为同时考虑流体和固体的一种研究方法，不可避免地需要涉及到流场数值计算和固场数值计算。本节和下一节分别阐述计算流体力学和计算固体力学的理论基础。计算流体动力学是结合了流体动力学理论、数值分析计算方法和计算机理论技术发展起来的一门学科，是近年来流体力学研究领域中的重要技术之一。它利用计算平台对流体流动、换热、燃烧以及其他复杂的物理现象进行模拟，求解其中的质量传递、能量传递、动量传递以及化学反应问题，可以获得流场内部压力、速度和温度等物理参数的空间分布情况与随时间变化情况。1.1.1流体动力学控制方程自然界中流体流动遵循质量守恒方程（连续性方程）、动量守恒方程（Navier–Stokes方程，简称N–S方程）以及能量守恒方程。本文研究对象水平轴潮流能水轮机在工作状态下，热量交换很小，可以忽略不计。周围水体可视为不可压缩牛顿流体，且不考虑流体密度和粘性的变化，于是，连续性方程可写为：∂（1.4）N–S方程为：∂（1.5）式中：ui表示i方向速度；i和j取值范围（1，2，3），分别代表笛卡尔坐标系中X、Y、Z这3个坐标轴方向；t表示物理时间；ρ表示流体密度；p表示压力；ν表示运动粘度。1.1.2湍流模型本文研究中，计算域内的流动大多处于湍流范围，直接求解方程（1.4）和（1.5），对计算资源的需求巨大，当前无法满足。需要对其进行近似和简化处理，分别采用了雷诺时均方法中的SSTk−（1）SSTk−在雷诺时均方法中，流动的瞬时速度和瞬时压力被分解为时均分量和脉动分量：u（1.6）p=（1.7）式中：ui表示i方向瞬时速度；ui表示i方向速度的时均分量；ui'将（1.6）和（1.7）代入（1.4）和（1.5）中，可以得到不可压缩湍流运动的雷诺时均方程：∂（1.8）∂（1.9）与方程（1.4）和（1.5）相比，雷诺时均方程多出了−ρui'uj'有关项，−ρui'uj−（1.10）k=（1.11）式中：δij是“Kroneckerdelta”符号（当i=j时，δ区别于流体动力粘度μ，湍动粘度μt取决于流动状态，是空间坐标的函数。为了求解湍动粘度μt，Launder和SpaldingADDINNE.Ref.{B0171D66-26E4-44DC-A511-FA1C0DAE7BAF}[148]引入湍动能耗散率εWilcoxADDINNE.Ref.{8589589F-7EA8-4C0A-98D5-42821A046A58}[149]引入了比耗散率ω，形成了k−ω两方程模型，可以较好地预测边界层流动、低雷诺数流动以及旋转流动，但缺点是计算结果对剪切层外的k和ω敏感性高ADDINNE.Ref.{C599D4F5-71FE-4190-88AE-6111F9C94DED}[150]。SSTk−ωADDINNE.Ref.{FC7216E5-7247-4001-9D82-03F110F55F04}[151]模型是一种混合模型，粘性力占主导地位的边界层内部区域，使用低雷诺数k−ω模型，逐渐过渡到边界层外部完全湍流区域使用高雷诺数k−ε模型。相比于k−ε模型，可以更精确地预测边壁附近的流动。相比于标准k−ω湍动能k和比耗散率ω分别由以下输运方程获得：∂（1.12）∂（1.13）式中：Jk和Jω分别表示k和ω的有效扩散率；Gk是因平均速度梯度造成k的产生项；Gω是ω的产生项；Hk和Hω分别表示k和ω因湍流造成的耗散；Sk和Q（1.14）F（1.15）Ψ（1.16）Q（1.17）式中：y表示到壁面的法向距离；Qω∗湍动粘度μt表达为kμ（1.18）F（1.19）Ψ（1.20）式中：Π是应变率幅值；α∗潮流能水轮机在工作时，叶片周围会出现流动分离现象，且主要基于翼型绕流，靠来流对叶片作用产生力矩使水轮机旋转。预测边壁附近流动的准确性对水轮机水动力性能结果非常重要。相较于其他雷诺时均模型，本文第四章和第五章的模拟采用了更合适的SSTk−（2）大涡模拟方法在湍流中，会产生不同尺度的涡。质量、动量和能量主要是通过大尺度涡传递。大尺度涡更依赖于具体求解问题，受流动涉及到的边界条件和几何形状影响大。小尺度涡受其影响小，更趋于各向同性，运动具有共性。大涡模拟方法的基本原理就是通过细密的网格筛选，直接计算能够被网格捕获的大尺度涡，而不能被网格捕获的小尺度涡通过亚格子模型处理。区别于直接数值模拟，大涡模拟只需要直接计算大尺度涡，从而可以使用相对更稀的网格和更大的时间步长。区分能够直接计算的大涡和需要建模的小涡就是网格过滤的过程：φ（1.21）式中：φ表示点x处的一个变量；φ表示过滤后的变量；V表示网格单元的体积；v表示网格单元尺度；Ω因此，大涡模拟为了能够直接计算尽可能多尺度的能量谱，需要足够精细的网格。图1.3描绘了某个节点处湍动能量谱与尺度的关系。需要强调的是，这在计算域内各个节点处是有区别的。图1.3湍动能量谱与尺度的对数关系Figure1.3Energyspectrumagainsttheinvertedlengthscale(logscales)对于不可压缩流动，过滤后的连续性方程和N–S方程分别为：∂（1.22）∂（1.23）式中：τij为亚格子尺度应力张量。过滤后的N–S方程（1.23）与雷诺时均方程（1.9）在形式上类似，区别在于这里的变量仍为瞬时值，而非时均分量，湍流应力的表达式也不同。亚格子应力τij属于未知量，为此，学者们提出了多种亚格子尺度模型ADDINNE.Ref.{EA09BA42-ECD1-48C7-B7F4-EFE2450FA657}[152-156]τ（1.24）μ（1.25）S（1.26）L（1.27）S（1.28）g（1.29）式中：τkk表示亚格子应力的各向同性部分，未被建模，但被加入到过滤后的静压项；Sij表示解析速度场的应变率张量；LS表示亚网格尺度混合长度；κ为冯·卡门常数；为了更好地捕捉计算域中的瞬时流动及其对叶片振动的影响，本文第五章的模拟采用了大涡模拟WALE模型。1.3计算结构动力学基本理论结构动力学基本控制方程为：M（1.30）式中：M表示结构质量矩阵；C表示结构阻尼矩阵；K表示结构刚度矩阵；d，d和d分别表示节点加速度向量，节点速度向量和节点位移向量；f上式可作为静、动力学平衡通式。在任意t时刻，上述方程可以视为考虑了惯性力项Md用于多自由度体系的非线性动力响应分析的逐步积分法有许多种，包括无条件稳定或条件稳定的。所谓无条件稳定是指无论选取的时间步长为多少，迭代求解过程都是稳定收敛的；而条件稳定是指求解过程的稳定性依赖于时间步长的选取。求解动力学平衡方程又分显式和隐式两种方法。两者的区别在于考虑动力学方程平衡状态的时刻。显式算法是完全基于t时刻的可用值，得到t+∆t时刻的瞬态解。中心差分时间积分法是应力分析应用中最常用的显式算法，优点是，在忽略阻尼项和使用主对角质量矩阵的情况下，不需要在每一步时间内求解方程组，只需要进行向量矩阵运算。这在数值计算中显得非常高效。然而，它仅在条件下稳定，稳定极限约等于弹性波穿过模型中最小单元尺寸的时间。由于这些限制，结构动力学中的显式方法主要用于大模型和短周期问题，例如碰撞研究。隐式算法为了增加时间步长上限，不仅基于t时刻的值，而且需要用到t+∆t时刻的相同变量值，从而获得t+∆t时刻的瞬态解。因此，必须求解非线性方程组。在结构问题中，隐式积分算法的时间步长通常比简单显式算法的稳定极限大一到两个数量级ADDINNE.Ref.{3DA02F60-DA7D-48F0-AAE1-3CD22CEDBD81}[157,158]。Newmark时间积分法ADDINNE.Ref.{E481F686-2FAF-479D-83EA-E8F797E68656}[159]是最常用的数值求解有限元半离散运动方程的直接时间积分方法之一。本文研究中的潮流能水轮机材料采用的是铝合金T6082-T6，属于线弹性材料。在线性结构动力学系统中，矩阵是常数，向量是时间的函数。式（1.30）可以改写为：M（1.31）式中：dt，dt和dt分别表示t时刻的节点加速度向量，节点速度向量和节点位移向量；fd（1.32）d（1.33）式中：a和b均为Newmark积分系数。结合（1.31）（1.32）和（1.33），t+∆t时刻的未知变量可以通过ta0（1.34）d（1.35）d（1.36）式中：各个积分常数分别取值为a0=1a∆t2，a1=b对于有限元半离散运动方程，选择合适的时间积分方法最重要的因素是精度、稳定性和耗散。在条件稳定的时间积分算法中，稳定性受选择的时间步长影响；而在无条件稳定的时间积分算法中，时间步长的选择可以不考虑稳定性。在Newmark时间积分法中，数值算法的耗散可以由积分系数b控制：b（1.37）a（1.38）当Newmark积分系数满足上述条件时，该时间积分法可以无条件稳定ADDINNE.Ref.{FA68B600-62F1-48FF-BB8F-C7A02ADBD0FC}[160]。1.4计算域离散与动网格技术流固耦合问题的研究既要考虑流场对结构的作用，又要考虑结构对流场的影响。前者通常体现在荷载的传递，后者则主要表现为结构变形的反馈。这就涉及到流场和固场的计算域中耦合交界面在数值计算过程中的运动变化。1.4.1计算区域的离散通常在进行数值计算之前，首先要离散化计算区域，即将空间上连续的计算区域划分成许多个子区域，并确定每个区域的节点，从而生成网格。然后，将偏微分格式的控制方程转化为各个节点上的代数方程组。求解这些代数方程组得到节点位置上的值，再根据节点值确定计算域内其他位置上的值。由于变量在节点之间的分布假设及推导离散方程的方法不同，形成了有限差分法（FiniteDifferenceMethod，简称FDM）、有限元法（FiniteElementMethod，简称FEM）和有限体积法（FiniteVolumeMethod，简称FVM）等不同类型的离散化方法。有限差分法是数值模拟偏微分方程最早采用的方法。其基本思路是将计算区域划分为差分网格，用有限个网格节点代替连续的求解域，利用Taylor级数展开等方法将控制方程中的导数项用函数值的差商代替，从而将连续函数的微分方程离散为网格节点上定义的差分方程组，每个方程中包含了该节点及附近节点上的待求解函数值。这是一种将微分问题转化为代数问题的近似数值解法。该方法的优点是其建立在经典的数学逼近理论的基础上，容易理解和接受，形式简单，使用方便；缺点在于仅当网格极其细密时，离散方程才满足积分守恒。有限元法借鉴了有限差分方法的离散处理内核，将计算区域划分为有限个互不重叠的微小单元，在每个单元内选择一些合适的节点构造插值函数，借助极值原理，将问题的控制方程转化为所有单元上的有限元方程。再将局部单元总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程组就得到各节点上待求的函数值。与有限差分法只考虑网格节点上因变量的数值相比，有限元法还需要假定因变量值在网格节点之间的变化规律，即构造插值函数，并将其作为近似解。有限元法的优点是对于求解区域的单元划分没有特别限制，因此适合处理具有复杂边界的计算域；缺点是计算量大，内存需求高，求解速度慢，因此常用于网格数量相对较少的固场计算。有限体积法将计算区域划分为有限个互不重叠的控制体，使得每个网格节点周围都有一个控制体，将待求解的微分方程对每一个控制体积分，得到一组离散方程。该方法的未知数是网格点上的因变量。每个控制体都能满足因变量的积分守恒，也就意味着整个计算区域满足积分守恒。该方法可视为有限差分法和有限元法的中间物。有限体积法只寻求网格节点上的因变量数值，但在寻求控制体积的积分时，必须假定因变量值在网格节点之间的分布，得到离散方程后，便可忘掉插值函数。必要时，可以对微分方程中不同的项采取不同的插值函数。有限体积法因其良好的守恒性和网格适应性，计算效率高等优点，广泛应用于计算流体力学领域。本文研究中，流体控制方程离散采用有限体积法，固体控制方程离散采用有限元法，两者可以通过交换耦合面的数值信息有效结合，从而实现流固耦合数值计算。1.4.2动网格方法部分研究人员认为流固耦合数值计算可以分为三个子场，流场和固场这两个真实物理场再加上动网格这一虚拟场ADDINNE.Ref.{4E98AF37-D096-4BE8-9949-1EA6F912C66F}[161]。固场求解是基于拉格朗日坐标系，着眼于质点。不设限制的情况下，每个网格节点都具有转动和平动6个自由度。网格节点的运动位移是固场中待求解计算的基本物理量之一，也正因为计算域内部节点间的相对运动，才导致了内应力的产生。流场求解大多基于欧拉坐标系，着眼于空间点，网格节点与计算空间保持固定。待求解的基本物理量通常为：速度、压力、温度和组分，并不计算网格节点位移。考虑计算区域中边界的变形或运动，则需要通过特殊手段来解决。为避免网格退化，甚至出现负体积，本文研究采用了动网格技术。本小节主要对基于Fluent平台上的动网格技术进行扼要阐述，详细理论公式可参考文献[161]，此处不再赘述。Fluent主要提供了三种方式控制网格节点位置的更新ADDINNE.Ref.{FEC47A90-87C1-47E2-8D66-0FF7DB0E858A}[162]。一是网格光顺（Smoothing）法。基本原理是将计算区域中的网格运动利用光顺方程近似为线弹性体运动，不改变网格的拓扑结构，保持节点之间的连接关系不变，不存在节点的生成或消除。适用于小幅度的网格运动问题，计算量小。二是动态铺层（Layering）法。基本原理是检测运动边界相邻网格的高度是否达到所设定的阈值决定网格分裂还是合并。适用于处理变形幅度大的边界运动，仅限二维问题中的四边形，三维问题中的三棱柱或六面体网格。三是网格重构（Remeshing）法。基本原理是在边界运动过程中，不断检测计算区域中的网格质量，并标记低质量网格，在动网格更新的过程中，对标记的网格重新划分，直至满足用户预设的扭曲率和尺寸要求。可以处理大幅度的网格运动，计算量大。适用于二维问题中的三角形网格、三维问题中的四面体网格以及边界层网格。通常与网格光顺法联合使用。考虑固体的变形和位移对流场的影响最直接的体现就是流场边界形状的变化，因此，流场网格在计算过程中需要更新。网格更新是流固耦合计算中关键环节之一，要求既能准确描述运动边界的运动变形状态，保持运动边界的连续，又能维持计算域网格单元合理的形状。在本文研究中，相对于网格尺寸，叶片变形大，并且为了适应复杂的模型几何外形，避免负体积网格的同时，保证运动更新后的网格质量，联合使用了光顺法和网格重构法。1.5数据交换方法在分区求解时，流固耦合交界面数据交换不同于同一物理场内交界面的数值共享，特点是该过程发生在不同的求解器之间，不同的离散化格式之间。以往研究人员在优化流固耦合分区求解法时，应用了不同插值方法处理交界面网格不匹配问题，例如，投射法ADDINNE.Ref.{9720A3D7-9F65-4DF1-BEDD-0E4711E581D8}[163-165]，最近邻插值ADDINNE.Ref.{93A81F4A-5674-4493-AFBD-73D034F0F89F}[166]，线性插值，径向基函数插值ADDINNE.Ref.{BCB508CE-75C7-41F5-8CE5-4A00D435C1D4}[167-170]。在传递广延性质（如位移）和强度性质（如荷载）时，通常又分别采用保形插值和守恒插值。这是根据区域在离散化后，所需传递的值是否具有“部分加和性”所决定的。另外，由于非匹配网格，发送方网格上的物理量，必须转换为基于目标网格类型的数据才能传递。尽管流场和固场相交界处轮廓是一致的，但各自单元类型和节点位置截然不同。对此非匹配网格，必须在不同网格进行邻域搜索和插值，从而实现数据传递。为了减少网格不匹配所消耗的计算资源和带来的误差，本文研究中的耦合交界面网格单元编号及其所包含的节点编号、坐标均保持一致。因此，流场计算得到的节点压力信息无须邻域搜索和插值即可精确一一映射至固场，而固场计算得到的节点位移同样精确反馈至流场。Fluent在动网格更新时与Ansys一样，不需要更改交耦合界面上的网格拓扑，从而保证数据交换的稳定性。Fluent和Ansys作为商业通用软件，可以输入或输出多种文件格式，这些文件格式可分为二进制文件和文本文件。其中inp文件是两者均可以输入输出的文件格式，可以通过文本编辑器清楚地看到其数据结构，方便读写各类信息，便于处理。此外，选择只输出所需要交换的流固耦合交界面上的节点坐标和/或荷载，从而降低文件大小、内存消耗，提高速度。尽管Fluent输出的inp文件与Ansys输出的inp文件在数据结构上不完全相同。但Fluent向Ansys传递数据时，只需以*CLOAD,OP=NEW为关键词检索inp文件中的节点荷载信息段，其中包含了节点编号、荷载方向、荷载大小，截取并保存，以备调用。同理，Ansys向Fluent传递数据时，只需检索inp文件中的节点编号和坐标，保存为可调用的UDF文件。1.6方柱下游柔性板流致振动研究模拟钝体尾流场内的流致振动被广泛应用于验证流固耦合数值计算方法的有效性ADDINNE.Ref.{2D738503-F75B-4835-BBD9-8E45C8737E9D}[171-179]。本节选择对方柱下游柔性板的流固耦合现象进行模拟。流体绕过方柱时会形成交替发放的漩涡，漩涡与方柱分离后，产生的脉动流体力激发柔性板的连续振动，柔性板的振动反过来又会改变方柱头部漩涡的发放状态。该算例最先是由Wall等人ADDINNE.Ref.{72807CED-81CF-4EDD-B2D7-FEFE5B8EBE67}[171,172]提出，作为二维模型，几何和物理性质均有明确的定义。各个区域都很容易建模，柔性板变形大，与流场相互作用明显，特别适合于研究耦合策略和算法。其准确性被不同研究组ADDINNE.Ref.{4E2C6879-63A4-4416-BFD9-635C863B68FA}[173,174,178,180]证实，已经成为检验耦合方法的基准算例之一。1.6.1计算模型本章算例模型和边界条件设置与文献[171，174，178]保持一致。在不可压缩流体中淹没了一个固定位置的方柱刚体。柔性板一端连接在方柱的下游侧面，未变形状态时处于水平位置，与来流方向平行。方柱尺寸为1cm×1cm×1cm，柔性板尺寸为4cm×1cm×0.06cm，计算区域如图1.4所示。流体的密度为1.18kg/m3，动力粘度为1.82×10-5sPa·s。柔性板的密度为100kg/m3，弹性模量为0.25MPa，泊松比为0.35。在柔性板自由端设置监测点A。文献[171]给出的是二维计算结果，采用完全匹配型流固耦合交界面，应用多项式界面预测方法和动态松弛因子不动点迭代算法。文献[174]同样进行的二维耦合计算，选择线性插值方法处理非匹配流固耦合交界面，应用Newton-Raphson方法作为耦合迭代算法。文献[178]给出了三维计算结果，选择SmartBucket插值方法处理非匹配流固耦合交界面，通过商业软件Ansysworkbench中的SystemCoupling模块连接Fluent和Ansys，其中，在SystemCoupling模块中通常采用的是固定松弛因子不动点迭代算法。本章采用完全匹配型流固耦合交界面，通过本文开发的接口程序连接Fluent和Ansys，对该模型进行三维耦合计算，计算过程中将上一次迭代时固场求解得到的交界面位移不做任何处理直接传递给流场求解。图1.4柔性板模型布局（单位：cm）Figure1.4Layoutoftheflexibleplate(unit:cm)1.6.2计算域离散化及边界条件尽管该算例的计算域划分六面体网格非常容易，考虑到柔性板会有较大的变形，特别是自由端附近网格容易发生畸变，出现负体积网格，导致计算失败。因此，流场网格更新需要采用光顺法配合网格重构法，流场采用四面体网格进行离散，如图1.5。固场网格同样了采用纯四面体网格划分。流场和固场网格单元数量分别为171万和24万，节点数量分别为32万和5万。流固耦合交界面网格的单元数量为4万，节点数量为2万。由于方柱设定为刚体，固场计算仅考虑柔性板部分。图1.5（c）展示了变形后的柔性板附近网格，可以看出该流固耦合计算方法能够使得网格在计算过程中保持良好的质量。（a）XZ平面网格（b）柔性板附近网格（c）变形后柔性板附近网格图1.5柔性板模型流场网格Figure1.5Fluidmeshoftheflexibleplate.(a)MeshinXZplane.(b)Mesharoundtheflexibleplate.(c)Mesharoundthedeformedflexibleplate.在流场中，进口设置为0.513m/s速度的均匀来流；出口设置为压力为零的出口边界；Z方向上下的两个平面设置为自由滑移壁面；Y方向前后的两个平面设置为可变形对称边界，变形方式约束在所处平面内；方柱表面设置为无滑移壁面；柔性板表面设置为可变形无滑移壁面，变形方式由UDF决定。在固场中，柔性板连接方柱的一端设置为固定约束；与流体接触的三个面设置为流固耦合交界面；设置大变形计算；时间积分为默认Newmark方法。计算时间步长设置为1×10-3s，流速和位移的迭代收敛标准为1×10-5。1.6.3结果分析文献[171]和文献[174]给出了初始振动过程中自由端A点在Z方向上的位移dAZ随时间变化规律，在图1.6（a）中与本文结果进行了比较。可以看出，在0s到1.5s这段时间内，本文得到的dAZ曲线波动幅值比文献[171]和文献[174]更大。这是由于计算过程中固场向流场传递交界面位移的处理方法不同导致的。文献[171]同时应用了多项式界面预测方法和动态松弛因子不动点迭