CN108920768A - 一种针对弹性薄壁结构的流固耦合算法 - Google Patents

Abstract

本发明公开了一种针对弹性薄壁结构的流固耦合算法：流体运动由光滑粒子流体动力学方法计算，并采用动力学边界粒子处理固壁边界，固体弹性响应由有限元方法计算；建立了有限元网格和动力学边界粒子间的空间位置对应关系，并提出了流固两项间流体力、速度和位移的信息传递方式；通过引入适当的结构刚度阻尼来增强流固耦合算的数值稳定性。本发明在流体力求解方面，省去了确定固壁边界法线方向的步骤，简化了对具有复杂物面形状的流固耦合问题的求解；在流固交界面处理方面，避免了固体域最外层网格因受力集中而发生网格畸变；在流固耦合算法的稳定性方面，通过设置结构刚度阻尼来稳定计算，避免了长时间计算所导致的物理量失真问题。

Description

一种针对弹性薄壁结构的流固耦合算法

技术领域

本发明涉及流体仿真及流固耦合的技术领域，特别涉及一种基于光滑粒子流体动力学方法和有限元方法的流固耦合计算方法。

背景技术

在土木、水利、海洋等工程领域常应用厚度相对长度能忽略不计的薄壁结构。这类结构在风、水流、波浪等荷载作用下产生大幅度弹性变形，结构变形又反过来影响上述环境因素，因此属于流固耦合问题。工程结构流固耦合问题的研究方法主要有三类：原型观测、模型试验和数值模拟。其中，原型观测成本高昂、观测耗时长，模型试验受比尺效应影响，难以同时满足所需的重力相似和弹性相似准则。近年来，随着理论研究的深入和计算机存储、计算能力的提升，数值模拟方法受到人们更多的关注。

常用的流体运动模拟方法有有限差分法、有限体积法、边界元法等。它们均是基于网格技术对流体域内部或边界进行离散，因此当流体域变形较大或出现流体分离及融合现象时易发生网格畸变，计算结果受到极大影响甚至导致计算失败。作为一种具有拉格朗日属性的无网格数值方法，光滑粒子流体动力学方法在处理上述问题时优势显著。

固体弹性变形的数值计算主要有模态叠加法和有限元法两类。前者计算简单且快速，但只能求解几何形状规则的结构的线弹性变形。后者的数学原理和求解过程虽相对复杂，但应用范围更为广泛。由于建筑工程结构往往受到强非线性流体荷载，其动力响应也就具有非线性，加之工程结构形式的多样性，因此适宜使用有限元方法求解。

但是，使用有限元方法同时求解流体运动和固体弹性变形时会遭遇网格畸变问题，使用光滑粒子流体动力学方法同时求解流体运动和固体弹性变形时会出现张力不稳定现象。

发明内容

本发明的目的是克服现有技术中的不足，提供一种针对弹性薄壁结构的光滑粒子流体动力学-有限元流固耦合算法，本发明算法首先采用光滑粒子流体动力学方法求解流体运动，并利用动力学边界粒子技术将流体作用力传递至固体域，之后采用有限元方法求解固体动力响应，并再次利用动力学边界粒子将有限元网格的位移和速度信息返还至流体域。依照以上步骤，实现对复杂流固耦合现象的数值模拟。

本发明所采用的技术方案是：一种针对弹性薄壁结构的流固耦合算法，其中，所述的弹性薄壁结构至少包括在土木、水利、海洋等工程领域应用的一切厚度相对较长度能忽略不计的结构构件，包括以下步骤：

步骤(1)、流体运动由流体连续性方程和纳维-斯托克斯方程控制，并利用流体状态方程闭合；采用光滑粒子流体动力学方法对由所述流体连续性方程、所述纳维-斯托克斯方程和所述流体状态方程组成的流体控制方程进行离散，得到常微分方程组；利用预报矫正数值积分方法求解所述常微分方程组，得到流体粒子的位移、速度和密度；利用动力学边界粒子方法施加流体固壁边界条件，基于流体粒子和动力学边界粒子间的动量守恒原理求得流体对固体的作用力；

步骤(2)、固体弹性变形由结构动力学方程控制；采用有限元方法离散所述结构动力学方程，得到常微分方程；利用纽马克数值积分方法求解离散后的所述常微分方程，得到固体上各有限元单元的位移、速度和加速度；

步骤(3)、建立起有限元网格和动力学边界粒子间的空间位置对应关系；

步骤(4)、建立起有限元网格和动力学边界粒子间的信息传递关系。

进一步的，步骤(1)所述的光滑粒子流体动力学方法是基于弱可压缩流体计算模式，所述的弱可压缩流体计算模式具体表现为引入流体状态方程来补充流体密度与压强间的关系；所述的光滑粒子流体动力学方法中还采用了Shepard密度过滤法来减弱弱可压缩流场中的高频压力振荡。

进一步的，步骤(2)所述的有限元方法是利用整齐排列的四节点四边形网格来离散固体域；所述的有限元方法的具体求解方法为：采用纽马克数值积分方法求解离散后的结构动力学方程，通过调整积分参数，使数值解法无条件稳定且具有二阶精度；时间步长的选取应保证一个结构自振周期内至少包含8个时间步，其中所述的结构自振周期指二阶或二阶以上；所述的有限元方法中设置结构刚度阻尼来吸收于弱流固耦合系统交界面处产生的多余能量。

进一步的，步骤(3)所述的有限元网格和动力学边界粒子间的空间位置对应关系为在每一个有限元网格的几何形心位置均设置一个动力学边界粒子，动力学边界粒子的位置随有限元网格位置的改变而改变。

进一步的，步骤(4)所述的建立起有限元网格和动力学边界粒子间的信息传递关系为利用动力学边界粒子将流体作用力传递至固体域，同时将固体域的变形和运动速度传递至流体域，具体包括以下步骤：

步骤(4-1)、在流体域内应用光滑粒子流体动力学方法求解得到流体粒子对动力学边界粒子的作用力；

步骤(4-2)、将步骤(4-1)中所求得的作用力散布在动力学边界粒子所在有限元网格的各个节点上，接着应用有限元方法求解固体网格节点的位置和速度；

步骤(4-3)、将步骤(4-2)中所求得的网格节点位置和速度集中在动力学边界粒子上，并传递至流体域；

至此完成了一个求解步骤，重复步骤(4-1)至步骤(4-3)，直至到达所设定的模拟时间上限。

本发明的有益效果是：

1、在流体力求解方面，相比传统的压力积分方法，本发明提出的流体力新算法省去了确定固壁边界法线方向的步骤，极大简化了对具有复杂物面形状的流固耦合问题的求解。

2、在流固交界面处理方面，相比已被提出的主从接触式耦合算法，本发明提出的光滑粒子流体动力学-有限元流固耦合算法能够将流体作用力更加深入、牢固地传递至固体域，避免固体域最外层网格因受力集中而发生网格畸变。

3、在光滑粒子流体动力学-有限元流固耦合算法的稳定性方面，相比已被提出的通过引入数值耗散来增强计算稳定性的方法，本发明提出的通过设置结构刚度阻尼来稳定计算的方法更加合理和贴近真实情况，同时有效避免了长时间计算所可能导致的物理量失真问题。

附图说明

图1：光滑粒子流体动力学固壁边界处理示意图；

图2：光滑粒子流体动力学-有限元流固耦合算法的流固交界面处理示意图；

图3：光滑粒子流体动力学-有限元流固耦合算法流程图。

图4：光滑粒子流体动力学-有限元流固耦合算法在淹没水平板防波堤工程中的一个具体应用示例。

具体实施方式

为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下：

一种针对弹性薄壁结构的流固耦合算法，其中，所述的弹性薄壁结构至少包括在土木、水利、海洋工程领域等应用的一切厚度相对长度能忽略不计的结构构件，包括以下步骤：

步骤(1)、流体运动由流体连续性方程和纳维-斯托克斯(Navier-Stokes，N-S)方程控制，并利用流体状态方程闭合；采用光滑粒子流体动力学(Smoothed ParticleHydrodynamics，SPH)方法对由所述流体连续性方程、所述纳维-斯托克斯(Navier-Stokes，N-S)方程和所述流体状态方程组成的流体控制方程进行离散，得到常微分方程组；利用预报矫正(Predictor-corrector)数值积分方法求解所述常微分方程组，得到流体粒子的位移、速度和密度；利用动力学边界粒子(DBPs)方法施加流体域固壁边界条件，基于流体粒子和动力学边界粒子间的动量方程求得流体对固体的作用力。此为模拟流体运动的方法。

其中，所述的光滑粒子流体动力学方法是利用一组携带独立物理量且可任意移动的粒子来离散流体域。本发明中，采用的光滑粒子流体动力学方法是基于弱可压缩流体计算模式。所述的弱可压缩流体计算模式具体表现为引入流体状态方程来补充流体密度与压强间的关系，继而回避了求解复杂的压力泊松(Poisson)方程。所述的流体动力学方法中还采用了Shepard密度过滤法来减弱弱可压缩流场中的高频压力振荡。

步骤(2)、结构弹性变形由结构动力学方程控制；采用有限元(Finite ElementMethod，FEM)方法离散所述结构动力学方程，得到常微分方程；利用纽马克(Newmark)数值积分方法求解离散后的所述常微分方程，得到固体网格单元的位移、速度和加速度。此为模拟结构弹性变形的方法。

其中，所述的有限元方法是利用整齐排列的四节点四边形网格来离散固体域。所述的有限元方法的具体求解方法为：采用纽马克法求解结构动力方程，通过调整积分参数，使数值解法无条件稳定且具有二阶精度；时间步长的选取应保证一个结构自振周期内至少包含8个时间步，其中所述的结构自振周期指二阶或二阶以上。所述的有限元方法中适当设置了适量结构刚度阻尼来吸收于弱流固耦合系统交界面处产生的多余能量。

步骤(3)、建立起有限元网格和动力学边界粒子间的空间位置对应关系。其中，所述的有限元网格和动力学边界粒子间的空间位置对应关系为在每一个有限元网格的几何形心位置均设置一个动力学边界粒子，动力学边界粒子的位置随有限元网格位置的改变而改变。

步骤(4)、建立起有限元网格和动力学边界粒子间的信息传递关系。所述的有限元网格和动力学边界粒子间的信息传递关系的具体方法为：利用动力学边界粒子将流体作用力传递至固体域，同时将固体域的变形和运动速度传递至流体域，具体包括以下步骤：

步骤(4-1)、在流体域内应用光滑粒子流体动力学方法求解得到流体粒子对动力学边界粒子的作用力；

步骤(4-2)、将步骤(4-1)中所求得的作用力散布在动力学边界粒子所在有限元网格的各个节点上，接着应用有限元方法求解固体网格节点的位置和速度；

步骤(4-3)、将步骤(4-2)中所求得的网格节点位置和速度集中在动力学边界粒子上，并传递至流体域；

至此完成了一个求解步骤，重复步骤(4-1)至步骤(4-3)，直至到达所设定的模拟时间上限。

本发明的原理在于：

(1)为准确求解流体对固体的作用力，本发明提出了基于动力学边界粒子方法和流体粒子与动力学边界粒子间动量守恒原理的流体力新算法。

(2)为实现光滑粒子流体动力学方法和有限元方法的耦合计算，本发明建立了动力学边界粒子与有限元网格的空间位置对应关系，动力学边界粒子所受流体力被平均分布在其邻近的有限元网格节点上，同时动力学边界粒子的位移与速度信息通过对邻近有限元网格节点的位移和速度信息的平均得到。

(3)为保证光滑粒子流体动力学-有限元耦合方法的稳定性，本发明提出了通过设置结构刚度阻尼来消除在流固交界面上因交错推进算法所产生的多余能量的方法。

下面结合具体公式对本发明做进一步解释：

图3给出了基于光滑粒子流体动力学方法和有限元方法的弹性薄壁结构流固耦合算法流程图，如图3所示，本发明提供一种基于光滑粒子流体动力学和有限元耦合方法的弹性薄壁结构动力响应求解技术，主要包括以下步骤：

1.流体运动的求解

本发明采用光滑粒子流体动力学方法计算流体运动。流体运动由连续性方程和纳维-斯托克斯控制，为回避求解压力泊松方程，假定流体微弱可压，继而引入流体状态方程。经光滑粒子流体动力学方法离散后的流体控制方程表示如下：

其中，下标i与j表示相互作用的粒子，ρ,p,m,c,u,r,g分别为流体粒子的密度、压强、质量、数值声速、速度、位移和重力加速度，ρ₀为流体参照密度，c₀为流体数值声速，W_ij为光滑函数，本发明中选用Wendland五次核函数形式，Π_ij为人工粘性，γ为常数，本例中取为7，t为时间。

本发明采用具有二阶精度的预报校正数值积分方法求解离散后的流体控制方程。为保证数值积分的稳定性，采用了可变时间步长计算模式，光滑粒子流体动力学方法的时间步长Δt_SPH由CFL条件、流体扩散速度和粒子加速度联合决定：

其中，c_i是流体粒子的数值声速，h是光滑长度，v是流体的运动黏度。

本发明应用Shepard密度过滤法来减弱弱可压缩流场内的非物理性高频压力振荡，即计算过程中每30个时间步利用下式重新指定流体粒子密度：

其中，为修正的流体粒子密度。

2.固体弹性响应的求解

本发明采用有限元方法计算固体弹性响应，网格选用四节点四边形单元。经有限元方法离散后的结构动力学方程表示如下：

[M]{α}+[C]{v}+[K]{δ}＝{F}

其中，[M]、[C]和[K]分别表示结构的质量、阻尼和刚度矩阵，[F]表示结构外力荷载矩阵，{a},{v},{δ}分别表示节点的加速度、速度和位移矢量。

本发明采用纽马克数值积分方法求解结构动力学方程：

其中，β_N和γ_N分别取0.25和0.5，在此条件下纽马克方法具有二阶精度、且无条件稳定。有限元方法的时间步长Δt_FEM按下式选取，以便至少能够捕捉结构的二阶振动形态：

其中，T⁽²⁾＝2π/ω⁽²⁾，ω⁽²⁾是结构的二阶固有频率。

本发明引入结构刚度阻尼来吸收弱于流固耦合系统交界面处产生的多余能量，刚度阻尼的计算式为：

[C]＝c_k[K]

其中，系数c_k按下式计算：

其中，ζ^(r)和ω^(r)分别是结构r阶振动模态下的材料阻尼比和固有频率。

3.光滑粒子流体动力学和有限元的耦合计算

在每一个有限元网格的几何形心位置均设置一个动力学边界粒子，即建立有限元网格和动力学边界粒子间的空间位置对应关系。计算过程中，首先利用光滑粒子流体动力学方法求解流体粒子对动力学边界粒子的作用力：

其中，f_i为动力学边界粒子受到的流体力。

接着，将动力学边界粒子受到的流体力向其所在的有限元网格的四个节点传递：

其中，f_e为网格节点受到的作用力。

利用有限元方法求解固体网格节点的位移和速度，并按照以下规则集成至动力学边界粒子：

其中，δ_e和v_e分别为网格节点的位移和速度，r_i和u_i分别为动力学边界粒子的位移和速度。

至此完成了一个求解步骤，重复上述步骤，直至到达所设定的模拟时间上限。

在海岸及海洋工程领域，淹没水平板被认为是一种环境友好的新型防波堤结构，并通常采用桩柱锚泊的方式。图4为波浪与淹没水平板相互作用的一个示例，此例中水平板由固定在其两端的桩柱锚泊。由于板长远大于板厚度，因此结构具有弹性特征，在波浪荷载作用下水平板中部发生挠曲变形。

尽管上面结合附图对本发明的优选实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，并不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护范围之内。

Claims (5)

1.一种针对弹性薄壁结构的流固耦合算法，其中，所述的弹性薄壁结构至少包括在土木、水利、海洋工程领域应用的一切厚度相对较长度能忽略不计的结构构件，其特征在于，包括以下步骤：

步骤(1)、流体运动由流体连续性方程和纳维-斯托克斯方程控制，并利用流体状态方程闭合；采用光滑粒子流体动力学方法对由所述流体连续性方程、所述纳维-斯托克斯方程和所述流体状态方程组成的流体控制方程进行离散，得到常微分方程组；利用预报矫正数值积分方法求解所述常微分方程组，得到流体粒子的位移、速度和密度；利用动力学边界粒子方法施加流体固壁边界条件，基于流体粒子和动力学边界粒子间的动量守恒原理求得流体对固体的作用力；

步骤(2)、固体弹性变形由结构动力学方程控制；采用有限元方法离散所述结构动力学方程，得到常微分方程；利用纽马克数值积分方法求解离散后的所述常微分方程，得到固体上各有限元单元的位移、速度和加速度；

步骤(3)、建立起有限元网格和动力学边界粒子间的空间位置对应关系；

步骤(4)、建立起有限元网格和动力学边界粒子间的信息传递关系。

2.根据权利要求1所述的一种针对弹性薄壁结构的流固耦合算法，其特征在于，步骤(1)所述的光滑粒子流体动力学方法是基于弱可压缩流体计算模式，所述的弱可压缩流体计算模式具体表现为引入流体状态方程来补充流体密度与压强间的关系；所述的光滑粒子流体动力学方法中还采用了Shepard密度过滤法来减弱弱可压缩流场中的高频压力振荡。

3.根据权利要求1所述的一种针对弹性薄壁结构的流固耦合算法，其特征在于，步骤(2)所述的有限元方法是利用整齐排列的四节点四边形网格来离散固体域；所述的有限元方法的具体求解方法为：采用纽马克数值积分方法求解离散后的结构动力学方程，通过调整积分参数，使数值解法无条件稳定且具有二阶精度；时间步长的选取应保证一个结构自振周期内至少包含8个时间步，其中所述的结构自振周期指二阶或二阶以上；所述的有限元方法中设置了结构刚度阻尼来吸收于弱流固耦合系统交界面处产生的多余能量。

4.根据权利要求1所述的一种针对弹性薄壁结构的流固耦合算法，其特征在于，步骤(3)所述的有限元网格和动力学边界粒子间的空间位置对应关系为在每一个有限元网格的几何形心位置均设置一个动力学边界粒子，动力学边界粒子的位置随有限元网格位置的改变而改变。

5.根据权利要求1所述的一种针对弹性薄壁结构的流固耦合算法，其特征在于，步骤(4)所述的建立起有限元网格和动力学边界粒子间的信息传递关系为利用动力学边界粒子将流体作用力传递至固体域，同时将固体域的变形和运动速度传递至流体域，具体包括以下步骤：

步骤(4-1)、在流体域内应用光滑粒子流体动力学方法求解得到流体粒子对动力学边界粒子的作用力；

步骤(4-2)、将步骤(4-1)中所求得的作用力散布在动力学边界粒子所在有限元网格的各个节点上，接着应用有限元方法求解固体网格节点的位置和速度；

步骤(4-3)、将步骤(4-2)中所求得的网格节点位置和速度集中在动力学边界粒子上，并传递至流体域；

至此完成了一个求解步骤，重复步骤(4-1)至步骤(4-3)，直至到达所设定的模拟时间上限。