«上一篇
文章快速检索     高级检索
下一篇»
  哈尔滨工程大学学报  2019, Vol. 40 Issue (4): 683-688  DOI: 10.11990/jheu.201709125
0

引用本文  

李裕龙, 廖洪烈, 胡湛, 等. 非匹配网格的三维流固耦合问题[J]. 哈尔滨工程大学学报, 2019, 40(4): 683-688. DOI: 10.11990/jheu.201709125.
LI Yulong, LIAO Honglie, HU Zhan, et al. 3D partitioned fluid-structure analysis based on non-matching meshes[J]. Journal of Harbin Engineering University, 2019, 40(4): 683-688. DOI: 10.11990/jheu.201709125.

基金项目

国家重点研发计划(2016YFC0402601);国家自然科学基金项目(51609269)

通信作者

欧素英, E-mail:ousuying@mail.sysu.edu.cn

作者简介

李裕龙, 男, 副研究员;
廖洪烈, 男, 工程师;
欧素英, 女,讲师,博士

文章历史

收稿日期:2017-11-23
网络出版日期:2018-11-09
非匹配网格的三维流固耦合问题
李裕龙 1, 廖洪烈 2, 胡湛 3, 欧素英 1     
1. 中山大学 海洋工程与技术学院, 广东 广州 119077;
2. 广州船舶及海洋工程设计研究院, 广东 广州 119077;
3. 中山大学 海洋科学学院, 广东 广州 119077
摘要:为研究流固耦合问题中的非匹配网格问题的准确可靠的数值计算方法,本文提出了一种基于非匹配网格的三维Common-Refinement方法,将其应用于离散后不可压缩流体与非线性超弹性体的非重叠子区域之间的接触面。采用Petrov-Galerkin有限元法离散不可压缩流体,并对大变形弹性结构体使用连续Galerkin有限元法进行离散。同时使用任意的Lagrangian-Eulerian(ALE)方法处理流固网格的大幅变形,并且采用全解耦的隐式分区方法去分别求解流固两相。为了满足两者之间液体和弹性耦合界面间牵引力的平衡条件,研究了共同细化方法的空间插值的准确性和可靠性。根据一系列的网格划分方案,通过改流体和固体网格之间的网格匹配系数系统地评估Common-Refinement方法的准确性和精度。将本方法应用于三维标准的圆柱体-弹性板问题,并与文献中标准解进行了对比。求解结果表明了这种方法在流固耦合问题中具有足够的准确性和可靠性。
关键词流固耦合    非匹配网格    Common-Refinement方法    区域分解    有限元    
3D partitioned fluid-structure analysis based on non-matching meshes
LI Yulong 1, LIAO Honglie 2, HU Zhan 3, OU Suying 1     
1. College of Marine Engineering and Technology, Sun Yat-sen University;
2. Guangzhou Marine Engineering Corporation Guangzhou 119077, China;
3. College of Marine Science, Sun Yat-sen University, Guangzhou 119077, China
Abstract: To investigate an accurate and reliable numerical calculation method to solve the non-matching mesh problem in fluid-structure coupling, this paper presents a 3D common-refinement method for non-matching meshes between discrete non-overlapping sub-domains of incompressible fluid and nonlinear elastic structure.The incompressible fluid flow was discretized by using a stabilized Petrov-Galerkin finite element method (FEM), and the large deformation structural formulation was discretized by using a continuous Galerkin FEM.An arbitrary Lagrangian-Eulerian formulation was used to process large deformation of fluid-structure mesh, and the fully decoupled implicit partition procedure was applied to solutions of the fluid and solid phases.To satisfy traction equilibrium condition along the fluid-elastic interface with non-matching meshes, this study investigates the accuracy and reliability of the spatial interpolation in the common refinement method (CRM).According to a series of mesh division schemes, this study systematically assessed the accuracy and precision of CRM by varying grid mismatch between the fluid and solid meshes, thereby demonstrating the second-order accuracy of CRM data transmission through uniform refinements of fluid and solid meshes along the interface.This method was further applied to a 3D benchmark problem of cylinder-elastic plate and compared with the standard solution in the literature.Results show that this method is accurate and reliable in solving the fluid-structure coupling problem.
Keywords: fluid-structure coupling    non-matching mesh    common refinement    domain decomposition    finite element method    

许多涉及多个物理场相互作用的科学和工程模拟通常需要精确和守恒的方法,通过非匹配离散网格来传输物理数据。这些问题包括电磁学、接触力学,共轭传热和流体-结构相互作用。特别在流固耦合应用方面,需要流体和固体中各自网格离散要能够精确地捕捉相互物理作用,其中涉及多尺度和复杂的多模态耦合动力学。这种非匹配网格在流固耦合应用中非常常见。非守恒的数据传输和局部不精确的插值和投影可能导致不稳定性预测和较差的估计,特别是当激励频率接近弹性结构的固有频率时。因此,准确的流固耦合模拟需要对非匹配网格的耦合界面边界条件进行精确而守恒的处理。

对于基于Eulerian-Lagrangian流固耦合系统,其数值模拟存在2种主要方法,即整体模拟和分区域模拟[1-3]。在整体模拟方法中,流体和结构方程使用完全耦合的方式一起联立求解。整体求解方式稳定度较高,即使存在非常强的附加质量效应时也如此,但是这种方法缺少灵活性和可塑性。另一方面,分区域算法则是按照一定的次序分别求解流体和结构2个子区域的控制方程[4-5],进而可以促进现有的流体和固体数值求解器之间进行合适的选择搭配,进行空间和时间的离散。

在一个典型的三维流固耦合问题分区域数值模拟中,耦合界面上的面网格拓扑一般并不相同。这意味着其非结构网格的数据结构也并不相同,由于离散化要求,它们的几何坐标也可能不一致。这种非匹配网格和相关的数据传输问题也存在于其他情形中,如自适应网格划分和多重网格法[6-7]

本文的主要目标有2个方面,其一是基于三维Common-Refinement方法对非匹配网格的误差进行系统的量化和分析,目的是减少匹配网格的必要性,进而设计一个更具有一般性的针对三维流固耦合交界曲面上的非匹配网格节点的投影方案。早期的工作[8-9]建立在一维和二维假定中,且流体假设为可压缩。对于弹性板问题,在可压缩流动系统的附加质量对时间长度具有依赖,而不可压缩系统的附加质量随着时间间隔趋于为零渐近为一个恒定值。本文结合不同的非匹配网格处理耦合界面,对非匹配网格的误差进行了系统的分析。

第2个目标是探讨了非匹配网格在三维流固耦合模拟中的应用。本文采用连续Galerkin(CG)有限元离散化方法对非线性超弹性体结构模型进行离散,采用Petrov-Galerkin有限元离散求解不可压缩流动。由于使用了区域分解策略,本文首先各自生成了流体区域和结构域的三维非结构网格,且两者网格的数据结构互相独立。流体作用在结构表面上的力作为表面牵引力矢量来处理,相对的结构体的位移给流体域的变形提供了边界条件。将三维Common-Refinement方法应用于三维标准钝体-弹性板问题,得到了满意的结果。

1 流体-结构控制方程

首先给出了完整的流体-结构系统的控制方程。对不可压缩流动的控制方程使用任意的拉格朗日-欧拉形式,弹性体结构动力方程则使用拉格朗日形式,以及相应的耦合界面的边界条件。

基于贴体移动边界方法模拟不可压缩流体与弹性结构的相互作用。假设Ω(t)fRd表示在时间t的流体域, 其中d为空间尺度。不可压缩粘性流体在Ω(t)f中的运动由N-S方程给出:

$ \begin{array}{*{20}{c}} {{\rho ^f}\frac{{\partial {{\mathit{\boldsymbol{\bar u}}}^f}}}{{\partial t}}\left| {_{\hat x}} \right. + {\rho ^f}\left( {{{\mathit{\boldsymbol{\bar u}}}^f} - {{\mathit{\boldsymbol{\bar w}}}^f}} \right) \cdot \nabla {{\mathit{\boldsymbol{\bar u}}}^f} = \nabla \cdot {{\bar \sigma }^f} + }\\ {\nabla \cdot {\mathit{\boldsymbol{\sigma }}^{{\rm{sgs}}}} + {\mathit{\boldsymbol{b}}^f},{\rm{on}}\;{\mathit{\Omega }^{\rm{f}}}\left( {\rm{t}} \right)} \end{array} $ (1)
$ \nabla \cdot {{\mathit{\boldsymbol{\bar u}}}^f} = 0,{\rm{on}}\;{\mathit{\Omega }^f}\left( t \right) $ (2)

式中:${\mathit{\boldsymbol{\bar u}}^f} = {\mathit{\boldsymbol{\bar u}}^f}\left( {{{\mathit{\boldsymbol{\bar x}}}^f}, t} \right) $$ {\mathit{\boldsymbol{w}}^f} = {\mathit{\boldsymbol{w}}^f}\left( {{\mathit{\boldsymbol{x}}^f}, t} \right)$分别表示流体和网格的速度; 空间点位置表示为$ {\mathit{\boldsymbol{x}}^f} \in {\mathit{\Omega }^f}\left( t \right);{\rho ^f}$是流体的密度; bf表示作用在流体上的体积力; σsgs来自于大涡模拟滤波; σf是牛顿流体的柯西应力张量,其表达式为:

$ {{\mathit{\boldsymbol{\bar \sigma }}}^f} = - {{\bar p}^f}I + {\mu ^f}\left( {\nabla {{\mathit{\boldsymbol{\bar u}}}^f} + {{\left( {\nabla {{\mathit{\boldsymbol{\bar u}}}^f}} \right)}^{\rm{T}}}} \right) $ (3)

式中:$ {\overline \rho ^f}$表示滤波后的流体应力;μf为流体的动力粘性系数;空间和时间坐标分别表示为xft。滤波后的N-S方程的弱形式可以写作:

$ \begin{array}{*{20}{c}} {\int_{{\mathit{\Omega }^f}\left( t \right)} {{\rho ^f}\left( {{\partial _t}{{\mathit{\boldsymbol{\bar u}}}^f} + \left( {{{\mathit{\boldsymbol{\bar u}}}^f} - {\mathit{\boldsymbol{w}}^f}} \right) \cdot \nabla {{\mathit{\boldsymbol{\bar u}}}^f}} \right)} \cdot {\phi ^f}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\Omega + }}\\ {\int_{{\mathit{\Omega }^f}\left( t \right)} {\left( {{{\mathit{\boldsymbol{\bar \sigma }}}^f} + {\mathit{\boldsymbol{\sigma }}^{{\rm{sgs}}}}} \right)} :\nabla {\phi ^f}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\Omega } = }\\ {\int_{{\mathit{\Omega }^f}\left( t \right)} {{\mathit{\boldsymbol{b}}^f} \cdot {\phi ^f}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\Omega }} + \int_{\mathit{\Gamma }_h^f\left( t \right)} {{h^f} \cdot {\phi ^f}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\Gamma }} } \end{array} $ (4)
$ \int_{{\mathit{\Omega }^f}\left( t \right)} {\nabla \cdot {{\mathit{\boldsymbol{\bar u}}}^f}q\left( x \right){\rm{d}}\mathit{\Omega }} = 0 $ (5)

式中:${\partial _t} $表示时间偏导算子$\partial \left( \cdot \right)/\partial t;{\phi ^f}$q分别是速度和压力的试函数;Γhf(t)表示流域边界;nf为流体域边界的法向量。流体域的更新使用任意拉格朗日欧拉(ALE)方法。

使用虚功原理来表示运动方程和作用在结构上的应力平衡。对于流体-固体系统中结构大幅度应变和动态变形,使用非线性超弹性体方程去描述[10]。假设固体的密度为ρs,固体的每一个点都由各自位置的矢量来描述。假设xsΩ0s表示弹性体的初始位置,同时使用ds(xs, t)∈Ωs(t)表示空间点xs在时刻t的位移。使用虚功原理后,弱形式的控制方程为:

$ \begin{array}{*{20}{c}} {\int_{\Omega _0^s} {\tau _{ij}^s\delta L_{ij}^s{\rm{d}}\mathit{\Omega }} - \int_{\Omega _0^s} {{\rho ^s}b_i^s\delta v_i^s{\rm{d}}\mathit{\Omega }} + \int_{\Omega _0^s} {{\rho ^s}\frac{{\partial u_i^s}}{{\partial t}}\delta v_i^s{\rm{d}}\mathit{\Omega }} - }\\ {\int_{\mathit{\Gamma }_2^s} {t_i^s\delta v_i^s{\eta ^s}{\rm{d}}\mathit{\Gamma }} = 0} \end{array} $ (6)

式中:τijs=Jsσijs是Kirchhoff应力,Js=det(Fs)是雅克比变形梯度张量,σijs是超弹性体的柯西应力;$\delta L_{ij}^s = \partial \delta v_i^s/\partial \phi _j^s$是虚速度梯度;bis是体积力;$\delta v_i^s = \delta v_i^s\left( {{\mathit{\boldsymbol{x}}^s}} \right)$是虚速度场。柯西应力σijs的表达式为:

$ \sigma _{ij}^s = \frac{{{\mu ^s}}}{{{{\left( {{J^s}} \right)}^{5/3}}}}\left( {B_{ij}^s - \frac{1}{3}B_{kk}^s{\delta _{ij}}} \right) + {K^s}\left( {{J^s} - 1} \right){\delta _{ij}} $ (7)

Bijs=FiksFjks是柯西-格林张量;μs是剪切模量。变形梯度张量为:

$ F_{ij}^s = {\delta _{ij}} + \partial d_i^s/\partial {x_j} $ (8)

这里给出三维流固耦合问题中耦合边界条件的简单描述。假设流体域为Ω0f,结构域为Ω0s,耦合边界${\mathit{\Gamma }^{fs}}\left( t \right) = \partial {\mathit{\Omega }^f}\left( t \right) \cap \partial {\mathit{\Omega }^s}\left( t \right)$,在耦合界面Γfs(t)上牵引力和速度必须保持连续。耦合界面边界条件可写作:

$ {{\mathit{\boldsymbol{\bar u}}}^f}\left( {{\mathit{\boldsymbol{\varphi }}^s}\left( {{\mathit{\boldsymbol{x}}^s},\mathit{\boldsymbol{t}}} \right),\mathit{\boldsymbol{t}}} \right) = {\mathit{\boldsymbol{u}}^s}\left( {{\mathit{\boldsymbol{x}}^s},\mathit{\boldsymbol{t}}} \right) $ (9)
$ \int_{{\varphi ^s}\left( {\gamma ,t} \right)} {{\mathit{\boldsymbol{\sigma }}^f}\left( {{\mathit{\boldsymbol{x}}^f},\mathit{\boldsymbol{t}}} \right)} \cdot \mathit{\boldsymbol{n}}{\rm{d}}\mathit{\Gamma }\left( {{\mathit{\boldsymbol{x}}^f}} \right) + \int_\gamma {{\mathit{\boldsymbol{t}}^s}{\rm{d}}\mathit{\Gamma }} = 0 $ (10)

式中:φs是位置矢量;ts是作用在结构边界上的流体牵引力矢量;us为弹性体在时间t的速度;n为耦合边界的外法线矢量。上述条件能够使流体速度完全等于界面的速度,并且弹性体的运动受流体牵引的支配,包括压力和剪切应力对物体表面的影响。

2 离散化后的流固耦合系统

借助于有限元法,流体空间Ωf被离散为互不重叠的有限单元Ωee=1, 2, …, nel, 其中nel为总的单元数。本文使用generalized-α方法去处理时间项[10]。伴随着generalized-α对时间项的处理,使用变分形式的流体控制方程可写作:

$ \begin{array}{*{20}{c}} {\int_{\mathit{\Omega }_{\rm{h}}^f\left( {{t^{n + 1}}} \right)} {{\rho ^f}\left( {{\partial _t}\mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,{\rm{n}} + \alpha _{\rm{m}}^f} + \left( { - w_{\rm{h}}^{f,{\rm{n}} + {\alpha ^f}}} \right) \cdot \nabla \mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,{\rm{n}} + {\alpha ^f}}} \right)} \cdot }\\ {{\phi ^f}{\rm{d}}\mathit{\Omega } + \int_{\mathit{\Omega }_{\rm{h}}^f\left( {{t^{n + 1}}} \right)} {\left( {\mathit{\boldsymbol{\bar \sigma }}_{\rm{h}}^{f,{\rm{n}} + {\alpha ^f}} + \mathit{\boldsymbol{\sigma }}_{\rm{h}}^{{\rm{sgs}},{\rm{n}} + {\alpha ^f}}} \right):\nabla {\phi ^f}{\rm{d}}\mathit{\Omega }} - }\\ {\int_{\mathit{\Omega }_{\rm{h}}^f\left( {{t^{n + 1}}} \right)} {\nabla \cdot \mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,{\rm{n}} + {\alpha ^f}}q{\rm{d}}\mathit{\Omega }} + \sum\limits_{{\rm{e}} = 1}^{{n_{{\rm{el}}}}} {\int_{{\mathit{\Omega }^e}} {\tau _m^f\left( {{\rho ^f}\left( {\mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,n + {\alpha ^f}} - w_{\rm{h}}^{f,{\rm{n}} + {\alpha ^f}}} \right) \cdot } \right.} } }\\ {\left. {\nabla {\phi ^f} + \nabla q} \right) \cdot \left( {{\rho ^f}{\partial _t}\mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,n + \alpha _{\rm{m}}^f} + {\rho ^f}\left( {\mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,n + {\alpha ^f}} - w_{\rm{h}}^{f,n + {\alpha ^f}}} \right) \cdot } \right.}\\ {\left. {\nabla \mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,n + {\alpha ^f}} - \nabla \cdot \bar \sigma _{\rm{h}}^{f,n + {\alpha ^f}} - \nabla \cdot \sigma _{\rm{h}}^{{\rm{sgs}},n + {\alpha ^f}} - {b^f}\left( {{t^{n + {\alpha ^f}}}} \right)} \right){\rm{d}}{\mathit{\Omega }^e} + }\\ {\sum\limits_{{\rm{e}} = 1}^{{n_{{\rm{el}}}}} {\int_{{\mathit{\Omega }_{\rm{e}}}} {\nabla \cdot {\phi ^f}\tau _c^f{\rho ^f}\nabla \cdot \mathit{\boldsymbol{\bar u}}_{\rm{h}}^{f,n + {\alpha ^f}}{\rm{d}}{\mathit{\Omega }^e}} } = }\\ {\int_{\mathit{\Omega }_{\rm{h}}^f\left( {{t^{n + 1}}} \right)} {{b^f}\left( {{t^{n + {\alpha ^f}}}} \right)} \cdot {\phi ^f}{\rm{d}}\mathit{\Omega } + \int_{\mathit{\Gamma }_{\rm{h}}^f\left( {{t^{n + 1}}} \right)} {{h^f} \cdot {\phi ^f}{\rm{d}}\mathit{\Gamma }} } \end{array} $ (11)

式中:第3~5行表示局部稳定项,τmfτcf为稳定参数,τmf定义为:

$ \tau _m^f = {\left[ {{{\left( {\frac{{2{\rho ^f}}}{{\Delta t}}} \right)}^2} + \left( {{}^\rho f} \right)2\left( {\bar u_{\rm{h}}^f - w_{\rm{h}}^f} \right) \cdot G\left( {\bar u_{\rm{h}}^f - w_{\rm{h}}^f} \right) + {C_I}\left( {{}^\mu {\rm{f}} + {\mu ^t}} \right)2G:G} \right]^{ - 1/2}} $ (12)

式中:CI 是常数,来自于元素的逆向估计;G是元素逆变度量张量的大小;μt是湍流粘性,其他各参数定义为:

$ \mathit{\boldsymbol{G}} = \frac{{\partial {\mathit{\boldsymbol{\xi }}^{\rm{T}}}}}{{\partial {\mathit{\boldsymbol{x}}^f}}}\frac{{\partial \mathit{\boldsymbol{\xi }}}}{{\partial {\mathit{\boldsymbol{x}}^f}}}\;\;\;\mathit{\boldsymbol{\tau }}_c^f = \frac{1}{{\left( {{\rm{tr}}\mathit{\boldsymbol{G}}} \right)\mathit{\boldsymbol{\tau }}_m^f}} $ (13)

式中xξ分别是空间坐标及其对应参数。稳定项有助于对流占优区域的速度场的稳定。

基于式(6),这里给出弹性体大幅变形的Galerkin有限元法离散的控制方程。这里使用超弹性属性的材料来描述大幅变形[11-12]。对于具有n个节点的单元,使用xia表示每个节点的坐标,其上标a表示节点范围,下标i表示空间范围(1~3),$\mathit{\boldsymbol{d}}_i^{{\rm{s}}, a}$表示每个节点的位移矢量,位移场和虚速度场分别表示为:

$ \mathit{\boldsymbol{d}}_i^s\left( {{\mathit{\boldsymbol{x}}^s}} \right) = \sum\limits_{a = 1}^n {{N^a}\left( {{\mathit{\boldsymbol{x}}^s}} \right)\mathit{\boldsymbol{d}}_i^{s,a}} ,\delta v_i^s\left( {{\mathit{\boldsymbol{x}}^s}} \right) = \sum\limits_{a = 1}^n {{N^a}\left( {{\mathit{\boldsymbol{x}}^s}} \right)\delta \mathit{\boldsymbol{v}}_i^{s,a}} $ (14)

式中:xs表示任意网格节点的坐标;Na表示形函数。将其代入式(6),可以得到非线性弹性体的控制方程:

$ \begin{array}{*{20}{c}} {\int_{\Omega _0^s} {{\rho ^s}{N^b}{N^a}\frac{{{\partial ^2}d_i^{s,b}}}{{\partial {t^2}}}{\rm{d}}\Omega } + \int_{\Omega _0^s} {\tau _{ij}^s\left[ {F_{pq}^s\left( {d_k^{s,b}} \right)} \right]\frac{{\partial {N^a}}}{{\partial {x_m}}}F_{mj}^{s, - 1}{\rm{d}}\Omega } - }\\ {\int_{\Omega _0^s} {{\rho ^s}N_i^s{N^a}{\rm{d}}\Omega } - \int_{{\mathit{\Gamma }^{fs}}} {t_i^s{N^a}{\eta ^s}{\rm{d}}\mathit{\Gamma }} = 0} \end{array} $ (15)

虚功方程(15)给出了n个非线性方程,包含了n个未知量。$\tau _{ij}^s\left[ {F_{pq}^s\left( {d_i^{s, a}} \right)} \right]$为Kirchhoff应力与变性梯度之间的函数关系。这组非线性方程在每个时间步内可以用Newton-Raphson迭代求解,时间推进使用标准Newmark方法处理,同样具有时间方向的二阶精度。需要注意的是方程组中的弹性矩阵需要在每个时间步内重构。

3 基于三维非匹配网格的Common-Refinement方法

流体和结构方程由速度和牵引力沿流体-结构耦合界面的连续性来进行数据交换,互相满足各自区域的边界条件。耦合界面的插值多项式的次数和流体与结构保持一致。

为了满足力的平衡条件,从流体的动量通量必须通过表面牵引力的形式转移到结构表面。假设NifNjs表示标准的有限元形函数,其中下标i表示在流体域中而下标j表示在固体域中,同时$\mathit{\boldsymbol{\tilde t}}_i^f \in {L^2}\left( {{\mathit{\Omega }^f}} \right)$$\mathit{\boldsymbol{\tilde t}}_i^s \in {L^2}\left( {{\mathit{\Omega }^s}} \right)$表示各自子区域的牵引力矢量,牵引力场tfts在边界ΓfΓs中可以表示为:

$ {t^f}\left( {{x^f}} \right) \approx \sum\limits_{i = 1}^{{m_f}} {N_i^ft_i^f,{t^s}\left( {{x^s}} \right)} \approx \sum\limits_{j = 1}^{{m_s}} {N_j^st_j^s} $ (16)

式中mfms分别是耦合界面上流体域和固体域的节点数目。一旦得到了tfNifNjs,进而便能够得到$\tilde t_j^s$。两侧同乘权函数积分后得:

$ \int_{{{\rm{\Gamma }}^{fs}}} {N_i^s{t^s}{\rm{d}}\mathit{\Gamma }} = \int_{{{\rm{\Gamma }}^{fs}}} {N_i^s{t^f}{\rm{d}}\mathit{\Gamma }} $ (17)

使用有限元近似后的牵引力平衡条件为:

$ \int_{{{\rm{\Gamma }}^{fs}}} {N_i^sN_j^s\tilde t_j^s{\rm{d}}\mathit{\Gamma }} = \int_{{{\rm{\Gamma }}^{fs}}} {N_i^sN_j^f\tilde t_j^f{\rm{d}}\mathit{\Gamma }} $ (18)

结构一侧的牵引力$\tilde t_j^s$为:

$ \tilde t_j^s = \mathit{\boldsymbol{M}}_{ij}^s{}^{ - 1}\mathit{\boldsymbol{f}}_i^s $ (19)

式中Ms表示耦合界面结构一侧的质量矩阵,表示为:

$ \mathit{\boldsymbol{M}}_{ij}^s = \int_{{{\rm{\Gamma }}^{fs}}} {N_i^sN_j^s{\rm{d}}\mathit{\Gamma }} $ (20)

节点力的矢量fjs为:

$ \mathit{\boldsymbol{f}}_i^s = \sum\limits_{j = 1}^{{m_f}} {\mathit{\boldsymbol{\tilde t}}_j^f\int_{{{\rm{\Gamma }}^{fs}}} {N_j^fN_i^s{\rm{d}}\mathit{\Gamma }} } $ (21)

当质量矩阵的构造只需要结构一侧的形函数时,荷载矢量积分由流体和固体的形函数共同组成。对于匹配的网格,这不会造成任何问题。然而对于非匹配网格,形状函数的不一致将导致两侧的不连续性。Common-Refinement方法是一类重要且特殊的数据结构,两套网格间Common-Refinement曲面是由同时对流体和固体子域的输入边界网格进行细分的多边形来构成。在有限元形式中,流体和固体界面网格的空间坐标可以写成:

$ {x^f} \approx \sum\limits_{i = 1}^{{m_f}} {N_i^f\left( x \right)x_i^f\;{\rm{on}}\;\mathit{\Gamma }_{\rm{h}}^f} $ (22)
$ {x^s} \approx \sum\limits_{j = 1}^{{m_s}} {N_j^s\left( x \right)x_f^s\;{\rm{on}}\;\mathit{\Gamma }_{\rm{h}}^s} $ (23)

Common-Refinement方法生成的网格拓扑由两侧输入的界面网格的交线来定义,或者被称作子单元。

在流固耦合模拟中,载荷矢量fjs在Common-Refinement网格中的子单元上计算方法为:

$ f_j^s = \sum\limits_{i = 1}^{{e_c}} {\int_{\sigma _l^c} {N_j^s{{\tilde t}^f}{\rm{d}}\mathit{\Gamma }} } $ (24)

式中:ec表示Common-Refinement网格上的子单元的总量;σic表示i个单元。利用高斯积分确定积分点位置及其权函数。Common-Refinement方法的实施方式为:首先根据输入的两侧的边界网格生成Common-Refinement网格的子单元,然后在各个子单元中遍历每个积分点,计算每个子单元的面积,然后将积分点与邻近的结构单元联系起来,对牵引力矢量按照式(24)进行积分。

在数值实现过程中,本文使用OpenMP结合MPI的方式,使用C语言编写了完整的数值并行流固耦合求解代码。整个流固耦合数值求解流程如图 1所示。流体域的载荷通过Common-Refinement网格传输给结构域,作为超弹性体的边界条件,超弹性体在当前载荷下产生的位移同样通过Common-Refinement网格传输给流体域,作为ALE网格求解和NS方程求解的边界条件,期间使用非线性力的迭代方式对当前流动载荷逐次修正。

Download:
图 1 流固耦合计算流程 Fig. 1 Flow chart of FSI simulation
4 非匹配网格误差分析

在典型的流固耦合模拟中,流体和固体的表面总是互相完整重合的。流体的载荷作用在到结构体表面上,而结构的位移作用在流体表面,进而在每一个时间步重复这样的数据传输。通常在数据传输中的误差来自于非匹配网格将数据从一个表面转移到另一个表面的误差。为了量化这个误差,建立了2个不同网格尺寸的相接触的完整曲面。把一个表面称为流体,而另一个表面为结构体。两者的直径D=1,高度H=100D,由图 3所示。

Download:
图 3 圆柱体-弹性板问题非匹配网格示意 Fig. 3 Non-matching mesh of cylinder-foil system

为了不失一般性,2个曲面在Z方向,也就是曲率无变化的方向划分wz均为36,而沿圆柱体的周向划分分别不同。这里定义单元数量为:

$ {A_s} = \frac{{{\rm{ \mathsf{ π} }}DH}}{{{w_z}{N_s}}},{A_f} = \frac{{{\rm{ \mathsf{ π} }}DH}}{{{w_z}{N_f}}} $ (25)

因此网格数量的非匹配的程度可以由As/Af来决定。根据这个标准,生成了如下参数范围的网格Ns, Nf∈{36, 54, 108}。这意味着网格非匹配的范围为As/Af∈[1/3, 3]。对于每一组网格,在位置上的流体表面网格的不同位置xf上施加指定的载荷,进而通过Common-Refinement方法传输到结构体表面。指定载荷为:

$ \begin{array}{*{20}{c}} {{t^s} = {t^s}\left( {\theta ,z} \right) = }\\ { - \left( {\frac{1}{2}{\rho ^f}U_\infty ^2\left( {1 - 4{{\sin }^2}\theta } \right) + {\rho ^f}gz} \right)\left[ {\begin{array}{*{20}{c}} {0.5\cos \theta }\\ {0.5\sin \theta }\\ 0 \end{array}} \right]} \end{array} $ (26)

式中:(θ, z)∈Γfs为当前圆柱面上的柱坐标位置矢量;U=1。柱坐标的示意图如图 2所示。指定载荷是通过沿圆柱体周围的流体按照静压力给出。假设(θjs, zjs)为节点j在圆柱表面的位置矢量,Tjs为传输后载荷,则传输相对误差ε1可以定义为:

Download:
图 2 误差分析圆柱体算例示意 Fig. 2 Presribed load transferred from fluid to solid across a cylinder surface
$ {\varepsilon _1} = \frac{{\sum\limits_j {T_j^s - {t^s}{{\left( {\theta _j^s,z_j^s} \right)}_2}} }}{{\sum\limits_j {{t^s}{{\left( {\theta _j^s,z_j^s} \right)}_2}} }} $ (27)

表 1给出了相对误差ε1根据不同的网格匹配比率(0.333 3~3.0)的计算结果。误差分析的结果表明Common-Refinement方法在任何网格比率的情况下都表现良好。值得注意的是误差在hs/hf>1和hs/hf < 1的情况下都是连续且良好的、保持对称的。

表 1 非匹配网格相对误差计算结果 Table 1 Relative error of non-matching meshes
5 圆柱体-弹性板问题

为了进一步评估基于非匹配网格的流固耦合求解器的可靠性,接下来将给出一个标准的基准问题:具有非常低的密度比的圆柱体-弹性板问题。首先根据文献中的可用数据验证的分区域流固耦合求解器的准确性,采用了文献[13]中提出的标准算法FSI-Ⅲ,一个标准的具有结构-流体低密度比的圆柱体-弹性板问题。圆柱体-弹性板问题为一个有限厚度的薄弹性板插入在一个固定圆柱体后,置于一个方形的流体域中所组成,根据圆柱体的直径,当前雷诺数Re=200。问题的相关参数见表 2

表 2 圆柱体-弹性板问题相关参数 Table 2 Parameters of cylinder-foil system

问题的边界条件设置和文献[13]保持一致。圆柱体表面,弹性板表面和流域顶部和底部为无滑移条件。入口处Γin指定抛物形速度剖面,具体为:

$ {u^f}\left( {0,y} \right) = 1.5\bar U\frac{{y\left( {H - y} \right)}}{{{{\left( {H/2} \right)}^2}}} = 1.5\bar U\frac{{4.0}}{{0.1681}}y\left( {0.41 - y} \right) $ (28)

式中:uf是速度矢量uf=[uf vf wf]的X方向分量, U为平均流速,H=4.1D为计算域高度,Γtop为顶部Γbottom为底部。为了描述流固耦合问题的物理特性,3个无量纲参数雷诺数,抗弯刚度和密度比分别定义如下:

$ \mathit{Re} = \frac{{{\rho ^f}\bar UD}}{{{\mu ^f}}},{K_B} = \frac{{E{w^3}}}{{12\left( {1 - {{\left( {{\nu ^s}} \right)}^2}} \right){\rho ^f}{{\bar U}^2}{L^3}}},{\rho ^r} = \frac{{{\rho ^s}}}{{{\rho ^f}}} $ (29)

式中:w为弹性板厚度;L为弹性板长度;E为杨氏模量;νs为泊松比。

使用有限元网格对流体和固体区域进行分解, 有一层边界层网格划分在结构体四周。3套不同密度的匹配网格用来进行网格收敛性的研究,其参数和计算结果如表 3所示。从表中可以看出M3网格的收敛性达到最好,其前端的最大位移和响应频率和标准算例的结果吻合良好。表明本文所采用的数值计算手段具有足够的精度和准确性。

表 3 网格收敛性计算结果 Table 3 Mesh convergence results

接下来研究不同网格比率的非匹配网格在流固耦合模拟中的准确度。典型的计算用非匹配网格见图 3。此时采用网格收敛性分析时的M3网格设置作为衡量基准。

选取网格比率在2.0,意味着选取从稠密到稀疏以及从稀疏到稠密一系列的网格来进行数值试验。不同网格比率hs/hf的计算结果如表 4所示。

表 4 不同网格比率下的流固耦合问题计算结果 Table 4 Numerical results of different mesh ratio

从计算结果能够明显的看出,不同网格比率算例hs/hf之间特征响应的结果非常接近。

6 结论

1) 当前的实现过程可以灵活的运用在不同领域的流固耦合问题中,同时因为可以任意划分和使用非匹配网格,因此具有足够的灵活性和稳健性。

2) 提出的三维数据传输方法基于L2模最小化了流体和结构网格之间的数据传输的误差,提供了一个准确的耦合方式。

3) 通过系统的误差分析,能够发现Common-Refinement方法对三维数据传输的一致性,不同网格比率的结算结果能够证明Common-Refinement方法具有的空间二阶精度。

4) 通过三维低密度比的基于非匹配网格的流固耦合问题的一系列数值实验,以及和标准算例的对比和验证,证明了这种方法的准确性、可靠性。

参考文献
[1]
HRON J, TUREK S. A monolithic FEM/Multigrid solver for an ALE formulation of fluid-structure interaction with applications in Biomechanics[M]//BUNGARTZ H J, SCHÄFER M. Fluid-Structure Interaction. Berlin, Heidelberg: Springer, 2006: 146-170. (0)
[2]
FELIPPA C, PARK K. Staggered transient analysis procedures for coupled mechanical systems:formulation[J]. Computer methods in applied mechanics and engineering, 1980, 24(1): 61-111. DOI:10.1016/0045-7825(80)90040-7 (0)
[3]
GURUGUBELLI P, JAIMAN R K. Self-induced flapping dynamics of a flexible inverted foil in a uniform flow[J]. Journal of fluid mechanics, 2015, 781: 657-694. DOI:10.1017/jfm.2015.515 (0)
[4]
JAIMAN R K, SEN S, GURUGUBELLI P S. A fully implicit combined field scheme for freely vibrating square cylinders with sharp and rounded corners[J]. Computers & fluids, 2015, 112: 1-18. (0)
[5]
JAIMAN R, GEUBELLE P, LOTH E, et al. Combined interface boundary condition method for unsteady fluid-structure interaction[J]. Computer methods in applied mechanics and engineering, 2011, 200(1/2/3/4): 27-39. (0)
[6]
De BOER A, VAN ZUIJLEN A H, BIJL H. Review of coupling methods for non-matching meshes[J]. Computer methods in applied mechanics and engineering, 2007, 196(8): 1515-1525. DOI:10.1016/j.cma.2006.03.017 (0)
[7]
DE BOER A, VAN ZUIJLEN A H, BIJL H. Comparison of conservative and consistent approaches for the coupling of non-matching meshes[J]. Computer methods in applied mechanics and engineering, 2008, 197(49/50): 4284-4297. (0)
[8]
JAIMAN R K, JIAO X, GEUBELLE P H, et al. Conservative load transfer along curved fluid-solid interface with non-matching meshes[J]. Journal of computational physics, 2006, 218(1): 372-397. DOI:10.1016/j.jcp.2006.02.016 (0)
[9]
JAIMAN R, JIAO X, GEUBELLE P H, et al. Assessment of conservative load transfer for fluid-solid interface with non-matching meshes[J]. International journal for numerical methods in engineering, 2005, 65(15): 2014-2038. (0)
[10]
JANSEN K E, WHITING C, HULBERT G. A generalized-α method for integrating the filtered Navier-Stokes equations with a stabilized finite element method[J]. Computer methods in applied mechanics and engineering, 2000, 190(3/4): 305-319. (0)
[11]
JANDRON M A, HURD R, BELDEN J, et al. Modeling of hyperelastic water-skipping spheres using abaqus/explicit[C]//Proceedings of the 2014 SIMULIA Community Conference. 2014. (0)
[12]
ZEMERLI C, LATZ A, ANDRÄ H. Constitutive models for static granular systems and focus to the Jiang-Liu hyperelastic law[R]. Fraunhofer (ITWM): Fraunhofer-Institut für Techno-und Wirtschaftsmathematik, 2012. (0)
[13]
TUREK S, HRON J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[M]//BUNGARTZ H J, SCHÄFER M. Fluid-Structure Interaction. Berlin, Heidelberg: Springer, 2006: 371-385. (0)