CN114757082B - 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法 - Google Patents

一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法 Download PDF

Info

Publication number
CN114757082B
CN114757082B CN202210215859.9A CN202210215859A CN114757082B CN 114757082 B CN114757082 B CN 114757082B CN 202210215859 A CN202210215859 A CN 202210215859A CN 114757082 B CN114757082 B CN 114757082B
Authority
CN
China
Prior art keywords
fluid
euler
lagrangian
solid
particles
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210215859.9A
Other languages
English (en)
Other versions
CN114757082A (zh
Inventor
王莉华
钱志浩
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tongji University
Original Assignee
Tongji University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tongji University filed Critical Tongji University
Priority to CN202210215859.9A priority Critical patent/CN114757082B/zh
Publication of CN114757082A publication Critical patent/CN114757082A/zh
Application granted granted Critical
Publication of CN114757082B publication Critical patent/CN114757082B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种基于拉格朗日‑欧拉稳定配点的流固耦合计算方法，包括：初始化处理；将流体信息从拉格朗日粒子映射到欧拉节点，并确定出自由表面与流固界面；将自由表面和流固界面附近的欧拉节点的流体速度外推；通过稳定配点法在欧拉节点上求解流固耦合问题的控制方程，并计算结构的速度；再次将自由表面和流固界面附近的欧拉节点的流体速度外推；将求解结果从欧拉节点映射至拉格朗日粒子，对粒子速度进行修正；更新拉格朗日粒子、重新分布流体粒子；判断当前时间步总时长是否大于设定的模拟时长，若判断为是，则进行后处理，否则更新当前时间步、返回迭代求解。与现有技术相比，本发明能够高效、准确、稳定地对流固耦合问题进行求解计算。

Description

一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法
技术领域
本发明涉及流固耦合计算技术领域，尤其是涉及一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法。
背景技术
涉及自由表面的流固耦合问题广泛存在于工程和自然界中，传统对该问题的研究主要包括理论分析和实验研究，而随着计算机技术的发展，数值模拟已成为当前研究该问题经济而有效的方法。在涉及流体自由表面的流固耦合问题中，一方面流体压力推动结构运动和变形，同时运动的结构又限制着流体流动；另一方面，流体域随着自由表面的演化和结构的运动而变化，同时自由表面又是由流体流动而决定；此外，流体在剧烈运动中还可能出现碎波现象。这些特点导致涉及自由表面的流固耦合问题模拟起来极具挑战性。
模拟流固耦合问题的数值方法主要可分为两类：整体求解方法和分区求解方法。其中，整体求解方法将流体和固体结合在同一求解域中，得到同时包括流体和固体的单一控制方程用于描述整个问题，并隐式地施加界面条件；此类方法使用统一的数值算法来求解流体和结构，使得该方法具有良好的精度和稳定性，然而流体和结构整体求解会增加计算规模，从而降低模拟效率。分区求解方法则将流体和结构作为两个独立的求解域，通过不同的数值方法进行求解，并通过迭代使其满足界面条件；该方法由于将流体域和结构域分开求解，减少了问题的计算规模，显著提高了模拟效率，但该类方法的缺点是精度和稳定性会受到迭代的影响，同时流体和结构的网格难以在界面上匹配。目前，分区求解方法比整体求解方法运用更为广泛。
根据问题域离散方式的不同，分区求解方法可进一步分为四大类：相容性网格方法、非相容性网格方法、粒子网格方法和粒子方法。在相容性网格方法中，流体网格和结构网格都沿界面构建并分别扩展到内部，界面处流体域和结构域共用节点。典型相容性网格方法是空间-时间方法和任意拉格朗日-欧拉方法，此类方法采用拉格朗日描述可以准确地捕捉流固界面，而在大变形或大范围运动问题中需要不断地进行网格重划分，计算成本高，网格质量难以保证，导致求解精度常常无法达到需求。相比之下，非相容性网格方法可以追溯到Peksin提出的浸没边界方法(IBM)。受IBM启发，众多学者开发了一系列非相容性网格方法，如改进的浸没边界法、浸没有限元法、虚拟域法、浸没界面法和浸没边界格子-玻尔兹曼方法。在这些方法中，流体域采用固定在空间上的欧拉网格进行离散化，而结构域采用固定在结构上的拉格朗日网格进行离散，并随着流体流动而运动。流固耦合的界面条件是通过向纳维尔-斯托克斯方程中添加等效力，同时流体对结构施加压力实现的。由于具有界面处容易施加无滑移条件，网格构造简单等优点，此类方法在求解流固耦合问题时简单有效。然而，此类方法需要非常小的切割单元来识别界面，导致时间步长过小而降低了计算效率，且界面附近需要非常细化的网格以提高界面附近的精度。综上可知，上述这些方法难以平衡精度和效率。
现有技术中，基于拉格朗日描述的粒子方法常用于极端变形和大范围运动问题的模拟，是模拟涉及自由表面的流体问题的一个很好的选择。因此，粒子-网格耦合方法在流固耦合问题中得到了发展，此类方法中流体采用拉格朗日点离散，固体采用拉格朗日网格离散，耦合界面信息交换通过插值实现。该种方法可以准确地捕捉自由表面，并成功地模拟了一些具有挑战性的流固耦合问题，包括结构入水和出水问题、海洋上浮体运动问题以及流体与柔性结构耦合问题等。然而，此类方法容易在界面附近发生粒子的非物理穿透问题。如果结构发生大变形时，此类方法仍然需要重新划分网格。另外，由于在两个区域中使用了不同的求解方法，此类方法中存在的时滞问题也限制着求解精度和效率。全粒子方法则能够彻底避免大变形下结构的网格重划分问题，该方法中，流体域和结构域都被离散为拉格朗日粒子，通过计算两种粒子之间的作用力，可以方便地考虑耦合效应[29]。此外，在不同域中使用相同的粒子方法可以很好地提高稳定性和改善时滞问题。全粒子方法中具有代表性的方法是光滑粒子水动力学法(SPH)。该方法在水下爆炸、接触爆炸、气泡坍塌和高速冲击等流固耦合问题中取得了众多成果。但该方法在数值模拟中存在拉伸不稳定性和精度低等问题。另一种典型方法是物质点法(MPM)，该方法也已应用于许多领域，例如高速冲击、不可压缩流动、流体-膜耦合问题和三维流固耦合问题。但MPM中的粒子在穿越背景网格时会出现不稳定性。总的来说，全粒子类方法在流固耦合领域取得了一定的成就，但仍然需要进一步改进以克服他们的缺点。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，以能够高效、准确、稳定地对流固耦合问题进行求解计算。
本发明的目的可以通过以下技术方案来实现：一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，包括以下步骤：
S1、初始化处理，包括初始化拉格朗日粒子和欧拉背景节点，其中，拉格朗日粒子用于描述流体运动；
S2、将流体信息从拉格朗日粒子映射到欧拉背景节点，并确定出自由表面与流固界面；
S3、将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S4、通过稳定配点法在欧拉背景节点上求解流固耦合问题的控制方程，并计算结构的速度；
S5、根据步骤S4的求解和计算结果，再次将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S6、将步骤S4的求解结果从欧拉背景节点映射至拉格朗日粒子，并对拉格朗日粒子速度进行修正；
S7、更新拉格朗日粒子、重新分布拉格朗日粒子；
S8、判断当前时间步总时长是否大于设定的模拟时长，若判断为是，则执行步骤S9，否则更新当前时间步，返回执行步骤S2；
S9、进行后处理，完成流固耦合计算过程。
进一步地，所述拉格朗日粒子包括流体粒子和结构粒子，其中，流体粒子用于附加流体信息，结构粒子用于附加刚体信息。
进一步地，所述步骤S2具体是将流体粒子映射到欧拉背景节点，使得欧拉背景节点附着在流体粒子上，并在一个时间步内随流体粒子一起发生运动。
进一步地，所述步骤S2中自由表面与流固界面具体为：
u＝uΓ(x,t) onΓ(x,t)
p＝0 onΠ(x,t)
ω＝ωk
其中，u＝[u1,u2]为二维问题中的流体速度矢量，p为流体压力，t为时间，v＝[v1,v2,ω]为固体结构旋转中心的速度矢量，v1、v2、ω分别为x方向速度、y方向速度和角速度，k为沿z方向的单位向量，ω为界面上的角速度矢量，Γ(x,t)、Π(x,t)和Λ(x,t)分别为固体边界、自由表面边界和流固界面，为界面上的速度矢量，Xr＝[Xr Yr]T为结构的旋转中心。
进一步地，所述步骤S4中流固耦合问题的控制方程由拉格朗日描述的N-S方程和牛顿定律公式组合而成，所述流固耦合问题的控制方程具体为：
Ωf∪Ωs＝Ω
ps＝[F,M]
其中，流固耦合问题域为Ω，流体域为Ωf，结构域为Ωs，u＝[u1,u2]为二维问题中的流体速度矢量，p为流体压力，t为时间，ρ为不可压缩流体的密度常数，ν为流体的运动粘度，f＝[0,-g]T为流体的体力矢量，g为重力加速度，v＝[v1,v2,ω]为固体结构旋转中心的速度矢量，v1、v2、ω分别为x方向速度、y方向速度和角速度，表示结构的质量矩阵，/>和/>分别为结构的质量和转动惯量，g＝[0,-g,0]T为结构的重力矢量，ps＝[F,M]表示施加在结构上的流体压力，F和M分别为作用在结构旋转中心上的等效轴向力和力矩，Xc＝[Xc Yc]T为结构质量中心，ns为结构的单位外法线向量，k为沿z方向的单位向量。
进一步地，所述步骤S4的具体过程为：
S41、将拉格朗日描述的N-S方程中的动量方程在与欧拉节点xl相关的规则子域中进行积分，可得积分方程：
其中，Ωl∈Ωf或
S42、引入向前差分对步骤S41得到的积分方程进行时域离散，并采用压力投影技术引入中间速度，可得：
根据压力投影技术，对上述方程进行分解，以分别求解中间速度求解流体压力p、求解最终速度un+1；
S43、根据流体压力p和最终速度un+1，计算得到流体对固体施加的压力再结合牛顿定律和界面上速度矢量的计算公式，进一步计算得到结构粒子的速度。
进一步地，所述步骤S42中求解过程具体包括以下步骤：
采用两组欧拉背景节点，其中一组用于流体速度计算，另一组用于压力计算；
通过基于相应源点的RK形函数，得到节点速度、中间速度和节点压力的近似值，将配点分配在与压力源点相同的位置，并定义和/>分别为流体域Ωf、固体边界Γ和流固界面Λ上的配点，NP，NQ和NR分别为对应的配点个数；
使用RK近似对流体速度和压力进行逼近，将结构受到的流体压力代入中间速度/>流体压力p、最终速度un+1对应的半离散方程及其相应的边界条件，并强制方程在与配点相关的子域中满足，以得到完全离散方程；
对完全离散方程求解得到中间速度流体压力p、最终速度un+1。
进一步地，所述求解中间速度对应的半离散方程及相应的边界条件为：
进一步地，所述求解流体压力p对应的半离散方程及相应的边界条件为：
其中，和/>均为子域Ωl的法向量，n＝[n1,n2]T为/>的外法线向量，θ为浸入流体中的子域一侧边界的长度分数，/>分别表示左边、右边、顶边、底边中心位置的值(·)。
进一步地，所述求解最终速度un+1对应的半离散方程及相应的边界条件为：
与现有技术相比，本发明针对具有自由表面的流固耦合问题，提出一种拉格朗日-欧拉稳定配点法进行求解计算，该方法中流体被离散为携带信息的拉格朗日粒子、覆盖整个运动空间的问题域被离散为均匀分布的欧拉背景节点；流体、固体和界面的耦合控制方程通过在欧拉节点上采用RK近似的无网格稳定配点方法求解；该方法中欧拉背景节点在每个时间步开始时都被重置为初始状态，不需要重构形函数；在流固耦合计算中，引入单元切割算法(Cell-cut method)，即在求解计算过程中使用不同算法，以同时求解流体压力和流固耦合力，从而避免传统耦合算法中繁琐的迭代求解过程，实现流体与结构的耦合求解。
本发明采用固定于空间的背景网格求解控制方程，避免了形函数的反复重构计算，效率高；而同时求解流体和固体控制方程，则避免了分区迭代求解，实现了强耦合计算，提高了精度和效率；此外，背景积分域内容易实现高阶积分，求解精度高，背景网格上应用稳定配点法获得的流场数值解满足局域守恒；本发明还能够实现自由表面和流固界面准确捕捉。由此保证了本发明方法能够高效、准确、稳定地对流固耦合问题进行求解计算。
附图说明
图1为本发明的方法流程示意图；
图2为本发明方法的具体应用过程示意图；
图3为拉格朗日-欧拉稳定配点法的离散方案示意图；
图4为流固耦合问题的几何模型示意图；
图5为稳定配点法中欧拉域离散方案示意图；
图6为实施例中二维圆柱入水问题的几何模型示意图；
图7a为t＝0.1825s时圆柱入水问题的流场分布示意图；
图7b为t＝0.2625s时圆柱入水问题的流场分布示意图；
图7c为t＝0.3325s时圆柱入水问题的流场分布示意图；
图7d为t＝0.4075s时圆柱入水问题的流场分布示意图；
图7e为t＝0.4675s时圆柱入水问题的流场分布示意图；
图7f为t＝0.6s时圆柱入水问题的流场分布示意图；
图8为自由表面演化的实验结果、LESCM结果和BEM结果的对比示意图；
图9为圆柱位移的实验结果、SPH结果和LESCM结果的比较示意图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
如图1所示，一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，包括以下步骤：
S1、初始化处理，包括初始化拉格朗日粒子和欧拉背景节点，其中，拉格朗日粒子用于描述流体运动，拉格朗日粒子包括流体粒子和结构粒子，流体粒子用于附加流体信息，结构粒子用于附加刚体信息；
S2、将流体信息从拉格朗日粒子映射到欧拉背景节点(具体是将流体粒子映射到欧拉背景节点，使得欧拉背景节点附着在流体粒子上，并在一个时间步内随流体粒子一起发生运动)，并确定出自由表面与流固界面；
S3、将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S4、通过稳定配点法在欧拉背景节点上求解流固耦合问题的控制方程，并计算结构的速度；
S5、根据步骤S4的求解和计算结果，再次将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S6、将步骤S4的求解结果从欧拉背景节点映射至拉格朗日流体粒子，并对拉格朗日流体粒子速度进行修正；
S7、更新拉格朗日粒子、重新分布拉格朗日流体粒子；
S8、判断当前时间步总时长是否大于设定的模拟时长，若判断为是，则执行步骤S9，否则更新当前时间步，返回执行步骤S2；
S9、进行后处理，完成流固耦合计算过程。
将上述方案应用于实际，其过程如图2所示。本技术方案提出了一种拉格朗日-欧拉稳定配点法(LESCM)，以对流固耦合问题进行模拟和求解计算。在LESCM中，问题域由一组拉格朗日粒子和一组欧拉背景网格节点进行离散，拉格朗日粒子用于描述流体运动，欧拉背景节点用于控制方程的求解，流固耦合问题也沿用该离散方式。如图3所示，流体由一组表示为的粒子进行离散，固体由一组表示为/>的粒子进行离散(该粒子仅位于固体表面)，其中Na和Nb分别表示离散流体和固体的拉格朗日粒子数。定义求解流体速度的欧拉背景节点为/>求解流体压力的背景节点为/>其中Nd和Nq是与流体速度和压力相关的节点总数。流体的所有信息(例如质量、速度、压力、位置)都附加到拉格朗日流体粒子上，而欧拉节点不携带长期信息。由于该问题中固体被视为刚体，因此这些刚体的信息(例如质量、转动惯量)被附加到每个固体结构上，且固体的拉格朗日粒子除了它们的坐标外不携带任何信息。在每个时间步上，流体的信息首先从粒子映射到背景网格节点。网格节点附着在粒子上，并在一个时间步内与它们一起运动。然后在网格节点上求解流体、结构和界面的耦合控制方程。在每个时间步结束时，背景网格节点被重置为初始状态。流体粒子通过从欧拉节点映射回来的数值解进行物理量的更新，结构粒子则通过求解出来的结构速度和加速度来更新。
如图4所示，定义流固耦合问题域为Ω，流体域为Ωf，结构域为Ωs，且满足Ωf∪Ωs＝Ω。拉格朗日描述的不可压缩N-S方程和牛顿定律表示的结构控制方程如下：
其中，u＝[u1,u2]是二维问题中的流体速度矢量，p是流体压力，t是时间，ρ是不可压缩流体的密度常数，ν是流体的运动粘度，f＝[0,-g]T是流体的体力矢量，g是重力加速度，v＝[v1,v2,ω]是固体结构旋转中心的速度矢量，其中v1，v2，ω分别为x方向速度、y方向速度和角速度，表示结构的质量矩阵，其中/>和/>分别是结构的质量和转动惯量。g＝[0,-g,0]T是结构的重力矢量，ps＝[F,M]表示施加在结构上的流体压力，F和M是作用在结构旋转中心上的等效轴向力和力矩，可以表示为：
其中Xc＝[Xc Yc]T是结构质量中心，ns是结构的单位外法线向量，k是沿z方向的单位向量。流动问题的固体边界条件和自由表面条件分别为：
u＝uΓ(x,t) onΓ(x,t) (5)
p＝0 onΠ(x,t) (6)
流体和结构之间的界面条件为
其中Γ(x,t)，Π(x,t)和Λ(x,t)分别代表固体边界、自由表面边界和流固界面，是界面上的速度矢量，表示为：
这里ω＝ωk，Xr＝[Xr Yr]T是结构的旋转中心。
在进行子域积分和时域离散时，由于每个时间步中，欧拉网格节点固结于拉格朗日粒子上，并在这个时间步中随拉格朗日粒子一起移动。这意味着拉格朗日描述中的N-S方程可以在欧拉节点上求解。如图5所示，SCM的求解过程将在欧拉节点上完成，将N-S方程中的动量方程在与欧拉节点xl相关的规则子域中进行积分，可得：
其中Ωl∈Ωf或引入向前差分对上诉积分方程进行时域离散，并采用压力投影技术引入中间速度，可得：
根据压力投影技术，可将上述方程分解为如下三步进行求解：
第一步：求解中间速度
第二步：求解流体压力p
pn＝0 on Π (15)
其中，
这里和/>是子域Ωl的法向量，n＝[n1,n2]T是/>的外法线向量，θ是浸入流体中的子域一侧边界的长度分数。/>分别表示左边、右边、顶边、底边中心位置的值(·)。
第三步：求解最终速度un+1
为了提高精度和稳定性，本技术方案在求解时采用两组背景节点，一组用于流体速度计算，另一组用于压力计算。节点速度、中间速度和节点压力的近似值可以通过基于相应源点的RK形函数近似。将配点分配在与压力源点相同的位置，并定义和/>分别为流体域Ωf、固体边界Γ和流固界面Λ上的配点。这里NP、NQ和NR是对应的配点个数，施加在结构上的流体压力可以改写为：
其中，矩阵J的计算基于公式(4)，使用RK近似对流体速度和压力进行逼近，并将结构受到流体压力代入上述压力投影技术中第一、二、三步的半离散方程及其相应的边界条件，并强制方程在与配点相关的子域中满足，可以得到以下完全离散方程：
其中：
α＝x或y且β＝Ⅰ,Ⅱ,Ⅲ或Ⅳ，/>和/>分别代表子域Ωl的左、右、上、下边的中心点。需要注意的是，压力的自由表面边界条件pn＝0并未直接在上述求解压力的步骤中直接施加，其具体施加方式可使用自由表面压力插值算法或基于RK近似的改进方法。在半离散方程(23)-(25)的积分中引入高斯数值积分，表示为：
其中，
分别是与配点Pl、Ql和Rl相关的积分点，ξL(L＝1,…Ni)是与积分点/>或/>或/>对应的权重，/>是界面片段/>上的第L个积分点。方程(27)中，/>可以从拉格朗日流体粒子到欧拉节点速度映射中获得。因此可以将方程(27)-(29)总结为：
其中，K1,K2,K3；F1，F2,F3的具体形式可通过方程(27)-(29)获得。通过求解方程(31)可获得流体的压力和速度对应的系数向量，进而可求得流体压力p和速度un+1，进而通过方程(22)可计算得流体对固体施加的压力，然后使用牛顿定律(3)和公式(8)可获得拉格朗日固体粒子的速度。
实施例
本实施例针对二维圆柱入水问题进行求解计算，如图6所示，一个半径为r0＝0.055m的二维半浮力圆柱体从Hc＝0.5m的高度落入最初平静的水中。“半浮力”是指圆柱体的重量等于该圆柱体完全浸没水中时浮力的一半。圆柱体和水的密度分别为ρc＝0.5×103kg/m3和ρw＝1.0×103kg/m3，水的运动粘度为υ＝1.01×10-6m2/s。由于圆柱体在接触水面之前自由下落，圆柱体在接触水面的瞬时速度可以由牛顿定律计算获得，即v0＝-2.955m/s。整个水池作为欧拉域被离散为302×364个欧拉背景节点，水池中的水使用235332个拉格朗日进行离散。根据稳定性条件，求解的时间步长被设置为Δt＝5.0×10-4s。
流场的粒子分布和对应的压力分布如7a～7f所示，该结果表明了LESCM得到的压力场光滑、且无压力不稳定问题。另外，LESCM可以很好地模拟水溅和碎波现象。图8展示了自由表面位置的实验结果、边界元方法(BEM)获得的数值结果和本发明LESCM获得的数值结果对比，该结果表明数值解与实验结果吻合良好。如图9所示，将圆柱位移的LESCM数值结果与实验结果、SPH数值结果进行了比较，表明本发明方法得到的数值解和实验数据能够良好吻合。
综上可知，本技术方案具有高精度、高效率和良好的稳定性，能够广泛地应用于流固耦合的工程应用中。

Claims (10)

1.一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，包括以下步骤：
S1、初始化处理，包括初始化拉格朗日粒子和欧拉背景节点，其中，拉格朗日粒子用于描述流体运动；
S2、将流体信息从拉格朗日粒子映射到欧拉背景节点，并确定出自由表面与流固界面；
S3、将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S4、通过稳定配点法在欧拉背景节点上求解流固耦合问题的控制方程，并计算结构的速度；
S5、根据步骤S4的求解和计算结果，再次将自由表面和流固界面附近的欧拉背景节点的流体速度外推；
S6、将步骤S4的求解结果从欧拉背景节点映射至拉格朗日粒子，并对拉格朗日粒子速度进行修正；
S7、更新拉格朗日粒子、重新分布拉格朗日粒子；
S8、判断当前时间步总时长是否大于设定的模拟时长，若判断为是，则执行步骤S9，否则更新当前时间步，返回执行步骤S2；
S9、进行后处理，完成流固耦合计算过程。
2.根据权利要求1所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述拉格朗日粒子包括流体粒子和结构粒子，其中，流体粒子用于附加流体信息，结构粒子用于附加刚体信息。
3.根据权利要求2所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述步骤S2具体是将流体粒子映射到欧拉背景节点，使得欧拉背景节点附着在流体粒子上，并在一个时间步内随流体粒子一起发生运动。
4.根据权利要求1所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述步骤S2中自由表面与流固界面具体为：
其中，为二维问题中的流体速度矢量，p为流体压力，t为时间，/>为固体结构旋转中心的速度矢量，v1、v2、ω分别为x方向速度、y方向速度和角速度，/>为沿z方向的单位向量，/>为界面上的角速度矢量，/> 和/>分别为固体边界、自由表面边界和流固界面，/>为界面上的速度矢量，/>为结构的旋转中心。
5.根据权利要求1所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述步骤S4中流固耦合问题的控制方程由拉格朗日描述的N-S方程和牛顿定律公式组合而成，所述流固耦合问题的控制方程具体为：
其中，流固耦合问题域为流体域为/>结构域为/> 为二维问题中的流体速度矢量，p为流体压力，t为时间，ρ为不可压缩流体的密度常数，ν为流体的运动粘度，为流体的体力矢量，g为重力加速度，/>为固体结构旋转中心的速度矢量，v1、v2、ω分别为x方向速度、y方向速度和角速度，/>表示结构的质量矩阵，/>和/>分别为结构的质量和转动惯量，/>为结构的重力矢量，/>表示施加在结构上的流体压力，/>和/>分别为作用在结构旋转中心上的等效轴向力和力矩，/>为结构质量中心，/>为结构的单位外法线向量，/>为沿z方向的单位向量。
6.根据权利要求5所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述步骤S4的具体过程为：
S41、将拉格朗日描述的N-S方程中的动量方程在与欧拉节点相关的规则子域中进行积分，可得积分方程：
其中，或/>
S42、引入向前差分对步骤S41得到的积分方程进行时域离散，并采用压力投影技术引入中间速度，可得：
根据压力投影技术，对上述方程进行分解，以分别求解中间速度求解流体压力p、求解最终速度/>
S43、根据流体压力p和最终速度计算得到流体对固体施加的压力/>再结合牛顿定律和界面上速度矢量的计算公式，进一步计算得到结构粒子的速度。
7.根据权利要求6所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述步骤S42中求解过程具体包括以下步骤：
采用两组欧拉背景节点，其中一组用于流体速度计算，另一组用于压力计算；
通过基于相应源点的RK形函数，得到节点速度、中间速度和节点压力的近似值，将配点分配在与压力源点相同的位置，并定义和/>分别为流体域/>固体边界Γ和流固界面Λ上的配点，NP，NQ和NR分别为对应的配点个数；
使用RK近似对流体速度和压力进行逼近，将结构受到的流体压力代入中间速度/>流体压力p、最终速度/>对应的半离散方程及其相应的边界条件，并强制方程在与配点相关的子域中满足，以得到完全离散方程；
对完全离散方程求解得到中间速度流体压力p、最终速度/>
8.根据权利要求7所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述求解中间速度对应的半离散方程及相应的边界条件为：
9.根据权利要求7所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述求解流体压力p对应的半离散方程及相应的边界条件为：
其中，和/>均为子域/>的法向量，/>为/>的外法线向量，θ为浸入流体中的子域一侧边界的长度分数，/>中i＝I、II、III、IV分别表示左边、右边、顶边、底边中心位置的值。
10.根据权利要求7所述的一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法，其特征在于，所述求解最终速度对应的半离散方程及相应的边界条件为：
CN202210215859.9A 2022-03-07 2022-03-07 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法 Active CN114757082B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210215859.9A CN114757082B (zh) 2022-03-07 2022-03-07 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210215859.9A CN114757082B (zh) 2022-03-07 2022-03-07 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法

Publications (2)

Publication Number Publication Date
CN114757082A CN114757082A (zh) 2022-07-15
CN114757082B true CN114757082B (zh) 2024-04-12

Family

ID=82325203

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210215859.9A Active CN114757082B (zh) 2022-03-07 2022-03-07 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法

Country Status (1)

Country Link
CN (1) CN114757082B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104268943A (zh) * 2014-09-28 2015-01-07 北京航空航天大学 一种基于欧拉-拉格朗日耦合方法的流体仿真方法
CN104850689A (zh) * 2015-04-30 2015-08-19 昆明理工大学 一种基于固定网格技术的流固耦合计算方法
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法
CN111339658A (zh) * 2020-02-25 2020-06-26 河海大学 基于拉格朗日无网格粒子法的水力瞬变模拟方法及设备

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10114911B2 (en) * 2010-05-24 2018-10-30 Fujitsu Limited Fluid structure interaction simulation method and apparatus, and computer-readable storage medium
US20160004802A1 (en) * 2014-07-03 2016-01-07 Arizona Board Of Regents On Behalf Of Arizona State University Multiscale Modelling of Growth and Deposition Processes in Fluid Flow

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104268943A (zh) * 2014-09-28 2015-01-07 北京航空航天大学 一种基于欧拉-拉格朗日耦合方法的流体仿真方法
CN104850689A (zh) * 2015-04-30 2015-08-19 昆明理工大学 一种基于固定网格技术的流固耦合计算方法
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法
CN111339658A (zh) * 2020-02-25 2020-06-26 河海大学 基于拉格朗日无网格粒子法的水力瞬变模拟方法及设备

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
A hybrid variational-collocation immersed method for fluid-structure interaction using unstructured T-splines;Hugo Casquero;《NUMERICAL METHODS IN ENGINEERING》;全文 *
A Lagrangian meshless finite element method applied to fluid–structure interaction problems;S.R. Idelsohn;《Computers and Structures 》;全文 *
A meshfree stabilized collocation method (SCM) based on reproducing kernel approximation;Lihua Wang;《ScienceDirect》;全文 *
中性悬浮大颗粒对湍槽流影响的数值研究;余钊圣;王宇;邵雪明;吴腾虎;;浙江大学学报(工学版)(第01期);全文 *
无网格稳定配点法及其在弹性力学中的应用;王莉华;《计算机力学学报》;全文 *
极端变形问题的物质点法研究进展;张雄;刘岩;张帆;陈镇鹏;;计算力学学报(第01期);全文 *
流固耦合理论与算法评述;陈锋;王春江;周岱;;空间结构(第04期);全文 *

Also Published As

Publication number Publication date
CN114757082A (zh) 2022-07-15

Similar Documents

Publication Publication Date Title
Bathe et al. Finite element developments for general fluid flows with structural interactions
CN112364571B (zh) 大型复杂耦合航天器动力学模型建模方法
CN106950853B (zh) 一种球形贮箱微重力环境下液体晃动的建模方法
CN111737897B (zh) 一种深海网箱高密度养殖鱼群的数值模拟方法
US20070043544A1 (en) Method for simulating stable but non-dissipative water
Roccia et al. Computational dynamics of flapping wings in hover flight: a co-simulation strategy
CN109446581B (zh) 一种波浪作用下浮体的水动力响应的测量方法及系统
CN110262526B (zh) 一种水下机器人空间6自由度自航操纵预报的类物理数值方法
He et al. A full-Eulerian solid level set method for simulation of fluid–structure interactions
CN108920768A (zh) 一种针对弹性薄壁结构的流固耦合算法
Qian et al. A highly efficient and accurate Lagrangian–Eulerian stabilized collocation method (LESCM) for the fluid–rigid body interaction problems with free surface flow
CN112507600B (zh) 一种移动粒子半隐式法的对称边界条件的构建方法
CN106896817A (zh) 一种基于粘滞阻尼振荡模型的多auv编队控制方法
Wang et al. Steady motion of underwater gliders and stability analysis
CN115310339A (zh) 基于物质点法的具有表面张力效应的固液耦合模拟方法
CN114757082B (zh) 一种基于拉格朗日-欧拉稳定配点的流固耦合计算方法
Wadi et al. Modeling and system identification of an autonomous underwater vehicle
CN116011350B (zh) 面向模块船运动特性分析的深度强化学习黑箱辨识方法
Liu et al. Hydrophobic and Hydrophilic Solid-Fluid Interaction
Li et al. A general computing platform for offshore renewable energy systems (OREGEN)
Deng et al. A level set based boundary reconstruction method for 3-D bio-inspired flow simulations with sharp-interface immersed boundary method
CN114692446A (zh) 一种矢量转动球坐标参数化的非线性壳有限元建模方法
CN109376461B (zh) 一种预测空泡迁移方向的方法
CN111931438A (zh) 基于化学势的晶格Boltzmann模型模拟液滴润湿现象的方法
CN104656445A (zh) 混合悬浮环境中的阻力效应补偿方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant