CN106547944A - 一种抑制可压缩两相流晃荡的隔板优化方法 - Google Patents

Abstract

本发明提供一种抑制可压缩两相流晃荡的隔板优化方法，包括：基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力；基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制，所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值；根据所述最优渗水率确定所述阻晃隔板形状。本发明实现了考虑带自由液面的可压缩两相流剧烈运动情况下，液舱内抑制可压缩流体运动的隔板形状的优化设计。

Description

一种抑制可压缩两相流晃荡的隔板优化方法

技术领域

本发明实施例涉及可压缩流体阻晃技术领域，尤其涉及一种抑制可压缩两相流晃荡的隔板优化方法。

背景技术

在海洋和航空航天领域，液态货物及液态燃料储液舱都是运载装备的重要组成构件之一，如海上浮式液态天然气开采终端(FLNG)的液化天然气储液舱(LNG)、运载火箭的燃料剂液态氢储液箱等。其共同特点都是船舶或航空器在运动过程中姿态在不停地改变，且储液舱内流体载液率处于不饱和状态，因此储液舱内具有自由液面的可压缩流体在外激励作用下会发生剧烈的运动，使得可压缩流体对液舱结构产生剧烈的冲击。同时，在砰击过程中，自由液面破碎时会包裹气体进入流场，气液混合易导致局部砰击荷载增大。大量事实说明：可压缩流体晃荡产生的局部砰击荷载对储液舱的结构安全和性能、对船舶及航空器的整体运行姿态和稳定性都会产生严重的影响。因此，为了保证海洋工程特种装备和航天运载器具有良好的结构安全性和工作稳定性，就必须对其内部流体的晃荡加以抑制，尽量避免液舱内可压缩流体发生剧烈的晃荡，避免流体晃荡与整体结构的固有振动发生耦合，使之不发生共振现象，减少储液舱结构的安全风险，保证船舶及航空器的正常工作。

目前阻晃隔板是用来研究晃荡抑制的较为有限的方法之一，但是针对阻晃隔板形状设计的研究较为匮乏，且板上挖空区域形状均为规则形状。为了更准确反映真实运载装备储液舱中可压缩流体的运动，亟待将具有自由液面的可压缩流体运动与阻晃隔板形状的优化设计系统地结合，准确且高效分析液舱内具有自由液面的可压缩流体剧烈运动下阻晃隔板的最优形状。

发明内容

本发明实施例提供一种抑制可压缩两相流晃荡的隔板优化方法，以克服上述技术问题。

本发明一种抑制可压缩两相流晃荡的隔板优化方法，包括：

基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力；

基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制，所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值；

根据所述最优渗水率确定所述阻晃隔板形状。

进一步地，所述基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力，包括：

建立添加阻晃隔板的液舱计算模型并采用六面体结构化网格划分计算域，隔板上每个网格均具有所述渗水率；

设置两种等温不混溶可压缩流体计算域、自由液面高度、气体和流体的密度、两相流的运动粘度系数物理参量；

根据所述流体计算域、液面高度和液舱主尺度计算自由液面晃荡的固有频率，流经隔板的流体相对速度始终为0，给出边界条件和初始条件；

根据所述自由液面晃荡的固有频率设置液舱晃动激励频率，使得激励频率在逼近共振频率，采用动网格技术模拟液舱受迫往复振荡的工况，阻晃隔板与液舱具有相同的运动；

建立所述可压缩流体运动的纳维斯托克斯方程组，选择可压缩流体体积法追踪和捕捉可压缩两相流界面，使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解；

对所述液舱侧壁处压力场积分，进而给出所述液舱侧壁所感受到的合力。

进一步地，所述可压缩流体运动的可压缩纳维斯托克斯方程组，包括：

(1)连续性方程

其中，ρ为密度[kgm^-3]，U为流速[ms^-1]，t为时间[s]；

针对可压缩流体，不同液相或气相的连续性方程：

其中，下标i为相，α为气体或流体在计算单元内的体积分数；

(2)可压缩纳维斯托克斯方程(动量方程)：

其中，P为压力[Pa]，μ为动力粘性系数[kgm^-1s^-1]，g为重力加速度[ms^-2]，Γ为两相流交界面，δ(x-x_s)为三维狄拉克δ函数，x是坐标向量，x_s是两相流交界面处坐标向量，表面张力项采用连续表面力(continuum surface force)计算，表面张力项的积分求解按下式计算：

其中，Ω_i为第i个计算单元，Γ∩Ω_i为Ω_i中的两相流交界面部分，σ为两相流交界面张力，κ为两相流交界面的曲率系数，n为两相流交界面单位法向矢量，α₁为相1的体积分数，两相流交界面的曲率根据体积分数α₁计算：

(3)状态方程：

ρ_i＝ρ_0i+ψ_iP(5)

其中，ρ₀为标准密度[kgm^-3]，ψ为压缩性[s²m^-2]，ψ＝1/c²，c为在流体中声音传播的速度[ms^-1]，根据不同相在空间的体积分数，有：

其中，α₁+α₂＝1，α₂为相2的体积分数；

(4)相连续方程：

其中，U_r为相对速度，U_r＝U₁-U₂，dgdt为压缩项，

进一步地，所述使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解，包括：

—在每一时刻的计算中，将上一时刻计算得到的压力场代入可压缩纳维斯托克斯方程，求出临时速度场；

—将动量方程代入连续性方程得到压力方程，并将临时速度场代入压力方程，求解出所述临时速度对应的压力场；

—求解每个单元的体积分数、密度和运动粘度系数；

—利用所述临时速度对应的压力场代入动量方程更新速度场；

—判断当前时刻，所述临时速度对应的压力场和更新后的速度场是否收敛；若不收敛，将更新后的速度场作为临时速度，返回至所述的“将动量方程带入……求解出所述临时速度对应的压力场”步骤，迭代计算；若收敛，结束当前时刻计算，进入下一时刻，重复上述计算步骤。

进一步地，所述基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制，所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值，包括：

三维液舱阻晃隔板形状优化的目标函数是两平行于隔板的垂直侧壁所感受到的合力的泛函，优化的最终目标是使两垂直侧壁所受压力合力最小，隔板上每个网格上的渗水率作为优化参数，优化约束条件是渗水率必须在0.0与1.0之间；

基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，形成隔板拓扑优化模块，对每一时刻，隔板渗水率进行优化，进而实现对所述阻晃隔板的渗水率的最优控制。

进一步地，所述基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，形成隔板拓扑优化模块，对每一时刻，隔板渗水率进行优化，进而实现对所述阻晃隔板的渗水率的最优控制，包括：

—设定初始渗水率η和优化开始时刻T_start；

—当激励时间t≤T_start(t＝t+Δt，Δt是时间步长，t为激励时长)，根据物理参数、初始条件和边界条件，调用带自由液面的可压缩计算流体力学模块计算流场参数，记录t＝T_start时，流场参数和初始目标函数J₀，J_min＝J₀；

—当激励时间t>T_start，在每一时刻上的计算中，首先调入优化算法，优化渗水率，用优化后的渗水率更新渗水率，再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*；

—判断当前时刻，渗水率是否最优，若J^*≥J_min，返回到所述的“当激励时间t>T_start……再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*”步骤，迭代计算，若J^*<J_min，J_min＝J^*，结束当前时刻计算，进入下一时刻(t＝t+Δt)；

—对每一时刻的渗水率进行优化，用上一时刻优化得到的渗水率替换原有的渗水率，再重复到计算当前时刻目标函数J^*步骤，直到结束。

本发明实施例，采用有限体积法离散可压缩流体纳维斯托克斯方程组，并结合可压缩流体体积法追踪和捕捉可压缩两相流界面，实现了带自由液面的可压缩两相流运动的数值算法；实现外激激励频率逼近自由液面固有频率，使舱内流体晃荡处于共振状态，模拟出可压缩裹气流体的砰击状态。基于带自由液面的可压缩两相流的砰击，采用拓扑优化和最优控制方法，实现一步一优化，进而对逼近共振工况下能有效抑制可压缩裹气流体晃荡隔板形状的优化设计。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作以简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明抑制可压缩两相流晃荡的隔板优化方法流程图；

图2为本发明具有阻晃隔板的三维液舱的示意图；

图3为本发明带自由液面的可压缩流体算法与隔板拓扑优化方法的计算流程图；

图4为本发明经过几个周期拓扑优化后液舱阻晃隔板形状正视图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明抑制可压缩两相流晃荡的隔板优化方法流程图，如图1所示，本实施例方法，包括：

步骤101、基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力；

进一步地，基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力，包括：

建立添加阻晃隔板的液舱计算模型并采用六面体结构化网格划分计算域，隔板上每个网格均具有所述渗水率；

设置两种等温不混溶可压缩流体，即气体和流体，气体和流体的计算域、液面高度、气体和流体的密度、两相流的运动粘度系数物理参量；

根据所述流体计算域、液面高度和液舱主尺度计算自由液面晃荡的固有频率，流经隔板的流体相对速度始终为0，给出边界条件和初始条件；

根据所述自由液面晃荡的固有频率设置液舱晃动激励频率，使得激励频率在共振频率附近。采用动网格技术模拟液舱受迫往复振荡的工况，阻晃隔板与液舱具有相同的运动；

建立所述可压缩流体运动的纳维斯托克斯方程组，选择可压缩流体体积法追踪和捕捉可压缩两相流界面，使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解；

对所述液舱侧壁处压力场积分，进而给出所述液舱侧壁所感受到的合力。

具体来说，本实施例中建立添加阻晃隔板的液舱计算模型(所有计算单位均为无量纲)，液舱长(x轴)×宽(y轴)×高(z轴)分别为2×2×2，坐标原点设在液舱底板中心处，并采用六面体结构化网格划分计算域网格，通过网格无关性检验确定壁面设置100×100×100个网格，隔板设置100×100个网格。六个壁面和隔板的边界条件均为固壁无滑移，类型均为动边界。具有阻晃隔板的三维液舱的示意图如图2所示。

设置流体属性为水(α₁＝1)，流体密度为1，运动粘度系数为10^-6，设置气体属性为空气(α₂＝0)，气体密度为0，运动粘度系数为1.589×10^-6。流体自由液面高度为0.7。液舱内自由液面的最低阶共振频率根据如下公式计算：

其中，i为阶数，σ_i为第i阶自由液面固有圆频率，T_i是第i阶自由液面固有周期，f_i是第i阶自由液面固有频率。流体晃荡的能量主要集中于自由液面最低阶固有圆频率附近，因此人们通常关注σ₁(即f₁或T₁)。a为液舱长度，h为液面高度。本实施例中，液舱晃动的速度激励为y＝0.1cos(2πft)＝0.1cos(3.76t)，t为激励时间，激励时长为N×T₁(N为周期数)，本例N＝10，即t＝16.7s，时间步长Δt＝0.05s。

隔板与液舱具有相同的运动属性，且流经隔板的流体相对速度始终为0，其他参量不作要求。设置所述阻晃隔板的渗水率η(指板上挖空面积与板面积的比值)的作用，渗水率设置在0.0至1.0区间内，即当η＝0.0时，表示板上完全不能透过水流，即一块完整的隔板将原始液舱分成容积和形状完全相同的两个独立舱室，η＝1.0时，表示完全透水，即原始矩形液舱内部光滑且无任何结构。

建立所述可压缩流体运动的纳维斯托克斯方程组，选择可压缩流体体积法追踪和捕捉可压缩两相流界面，使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解；

对所述液舱侧壁处压力场积分，进而给出所述液舱侧壁所感受到的合力。

进一步地，所述可压缩流体运动的可压缩纳维斯托克斯方程组，包括：

(1)连续性方程

其中，ρ为密度[kgm^-3]，U为流速[ms^-1]，t为时间[s]；

针对可压缩流体，不同液相或气相的连续性方程：

其中，下标i为相，α为气体或流体在计算单元内的体积分数；

(2)动量方程，也即可压缩纳维斯托克斯方程：

其中，P为压力[Pa]，μ为动力粘性系数[kgm^-1s^-1]，g为重力加速度[ms^-2]，Γ为两相流交界面，δ(x-x_s)为三维狄拉克δ函数，x是坐标向量，x_s是两相流交界面处坐标向量，表面张力项采用连续表面力(continuum surface force)计算，表面张力项的积分求解按下式计算：

其中，Ω_i为第i个计算单元，Γ∩Ω_i为Ω_i中的两相流交界面部分，σ为两相流交界面张力，κ为两相流交界面的曲率系数，n为两相流交界面单位法向矢量，α₁为相1的体积分数，两相流交界面的曲率根据体积分数α₁计算：

(3)状态方程：

ρ_i＝ρ_0i+ψ_iP(5)

其中，ρ₀为标准密度[kgm^-3]，ψ为压缩性[s²m^-2]，ψ＝1/c²，c为在流体中声音传播的速度[ms^-1]，根据不同相在空间的体积分数，有：

其中，α₁+α₂＝1，α₂为相2的体积分数；

(4)相连续方程：

其中，U_r为相对速度，U_r＝U₁-U₂，dgdt为压缩项，

进一步地，使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解，包括：

—在每一时刻的计算中，将上一步计算得到的压力场代入可压缩纳维斯托克斯方程，求出临时速度场；

—将动量方程代入连续性方程得到压力方程，并将临时速度场代入压力方程，求解出所述临时速度对应的压力场；

—求解每个单元的体积分数和其他相关物理量；

—利用所述临时速度对应的压力场代入动量方程更新速度场；

—判断当前时刻，所述临时速度对应的压力场和更新后的速度场是否收敛；若不收敛，将更新后的速度场作为临时速度，返回至所述的“将动量方程带入……求解出所述临时速度对应的压力场”步骤，迭代计算；若收敛，结束当前时刻计算，进入下一时刻，重复上述计算步骤。

步骤102、基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制，所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值；

三维液舱阻晃隔板形状优化的目标函数是两平行于隔板的垂直侧壁所受到的压力合力P的泛函，优化的最终目标是使两垂直侧壁所受压力合力最小。优化目标函数：

其中，J为目标函数，P_l0为优化前左侧壁感受到的压力合力，P_l为优化后左侧壁感受到的压力合力，P_r0为优化前右侧壁感受到的压力合力，P_r为优化后右侧壁感受到的压力合力，a是液舱长度，c和d分别为所述左侧壁面网格单元数，j和k分别为所述右侧壁面网格单元数，c和j分别为y轴方向，d和k分别为z轴方向。隔板的渗水率η为优化参数，优化约束条件是：

其中，m和n所述隔板网格划分单元数目，m为y轴方向，n为z轴方向。

本实施例中阻晃隔板形状的优化流程是：

—设定初始渗水率η和优化开始时刻T_start；

—当激励时间t≤T_start(t＝t+Δt，Δt是时间步长，t为激励时长)，根据物理参数、初始条件和边界条件，调用带自由液面的可压缩计算流体力学模块计算流场参数，记录t＝T_start时，流场参数和初始目标函数J₀，J_min＝J₀；

—当激励时间t>T_start，在每一时刻上的计算中，首先调入优化算法，优化渗水率，用优化后的渗水率更新渗水率，再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*；

—判断当前时刻，渗水率是否最优，若J^*≥J_min，返回到所述的“当激励时间t>T_start……再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*”步骤，迭代计算，若J^*<J_min，J_min＝J^*，结束当前时刻计算，进入下一时刻(t＝t+Δt)；

—对每一时刻的渗水率进行优化，用上一时刻优化得到的渗水率替换原有的渗水率，再重复到调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*步骤，直到结束。

根据所述最优渗水率确定所述阻晃隔板形状。整体计算流程如图3所示。最优控制问题是随时间而变的动态最优化问题。本实施例中针对抑制可压缩两相流运动的阻晃隔板形状的最优控制的具体求解过程如下：

已知受控系统的状态方程：

其中，x＝(x₁,x₂,···,x_n)^T为n维状态向量，u＝(u₁,u₂,···,u_r)^T为r维控制向量。t为自变量，f＝(f₁,f₂,···,f_n)^T为x、u和t的n维函数向量。

系统的起始状态(初态)如下，通常是已知的：

x(t₀)＝x₀(11)

其中，x(t₀)是初始时刻的值，t₀是起始时刻。

对系统所要达到的状态(末态)的要求可用如下末态约束条件表示：

其中，x(t_f)是结束时刻的值，t_f是结束时刻。

满足上述末态约束的状态集合的目标集，记为M，表述如下：

M＝{x(t_f):x(t_f)∈Rⁿ,g₁(x(t_f),t_f)＝0,g₂(x(t_f),t_f}≤0}(13)

从给定初态x(t₀)到目标集M的转移可通过不同的控制律u(t)来实现。为了在各种可行的控制律中找到一种效果最好的控制，采用性能指标函数J评价控制效果好坏，表述如下：

其中，S(x(t_f),t_f)为末态型性能指标。为积分型性能指标。

求一个容许控制u(t)∈U，t∈[t₀,t_f]，使系统(10)从给定初态(11)出发，在t_f>t₀时刻转移到目标集(13)的转移，使得性能指标(14)为最小。这就是最优控制问题的思想和程序编写过程，其u^*(t)为最优控制，其性能指标J^*＝J[u^*(·)]则被称为最优性能指标。

步骤103、根据所述最优渗水率确定所述阻晃隔板形状。

经过本发明通过几个周期拓扑优化后的液舱阻晃隔板形状示意图，如图4所示。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (6)

1.一种抑制可压缩两相流晃荡的隔板优化方法，其特征在于包括：

基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到液舱内加入阻晃隔板后可压缩流体晃荡的流速和液舱侧壁感受到的合力；

基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制；所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值；

根据所述最优渗水率确定所述阻晃隔板形状。

2.根据权利要求1所述的方法，其特征在于所述基于有限体积法并结合流体体积法，求解可压缩纳维斯托克斯方程组，得到所述液舱内加入阻晃隔板后可压缩流体晃荡的流速和所述液舱侧壁感受到的合力，包括：

建立添加阻晃隔板的液舱计算模型并采用六面体结构化网格划分计算域，隔板上每个网格均具有所述渗水率；

设置两种等温不混溶可压缩流体计算域、自由液面高度、气体和流体的密度、两相流的运动粘度系数物理参量；

根据所述流体计算域、液面高度和液舱主尺度计算自由液面晃荡的固有频率，流经隔板的流体相对速度始终为0，给出边界条件和初始条件；

根据所述自由液面晃荡的固有频率设置液舱晃动激励频率，使得激励频率在逼近共振频率，采用动网格技术模拟液舱受迫往复振荡的工况，阻晃隔板与液舱具有相同的运动；

建立所述可压缩流体运动的纳维斯托克斯方程组，选择可压缩流体体积法追踪和捕捉可压缩两相流界面，使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解；

对所述液舱侧壁处压力场积分，进而给出所述液舱侧壁所感受到的合力。

3.根据权利要求2所述方法，其特征在于，所述可压缩流体运动的可压缩纳维斯托克斯方程组，包括：

(1)连续性方程

$\frac{\&part;\mathrm{\&rho;}}{\&part;t}+\&dtri;\&CenterDot;\left(\mathrm{\&rho;}U\right)=0---\left(1\right)$

其中，ρ为密度[kgm-³]，U为流速[ms-¹]，t为时间[s]；

针对可压缩流体，不同液相或气相的连续性方程：

$\frac{\&part;{\mathrm{\&rho;\&alpha;}}_i}{\&part;t}+\&dtri;\&CenterDot;\left({\mathrm{\&rho;\&alpha;}}_i{U}_i\right)=0---\left(2\right)$

其中，下标i为相，α为气体或流体在计算单元内的体积分数；

(2)可压缩纳维斯托克斯方程(动量方程)：

$\begin{array}{c}\frac{\&part;\mathrm{\&rho;}U}{\&part;t}+\&dtri;\&CenterDot;\left(\mathrm{\&rho;}UU\right)=-\&dtri;(P+\frac{2}{3}\mathrm{\&mu;}\&dtri;\&CenterDot;U)+\&lsqb;\&dtri;\&CenterDot;(\mathrm{\&mu;}\&dtri;U)+\&dtri;U\&CenterDot;\&dtri;\mathrm{\&mu;}\&rsqb;\\ +\mathrm{\&rho;}g+\underset{\mathrm{\&Gamma;}}{\&Integral;}\mathrm{\&sigma;}\mathrm{\&kappa;}\mathrm{\&delta;}(x-x_s)nd\mathrm{\&Gamma;}\left(x_s\right)\end{array}---\left(3\right)$

其中，P为压力[Pa]，μ为动力粘性系数[kgm-¹s-¹]，g为重力加速度[ms-²]，Γ为两相流交界面，δ(x-x_s)为三维狄拉克δ函数，x是坐标向量，x_s是两相流交界面处坐标向量，表面张力项采用连续表面力(continuum surface force)计算，表面张力项的积分求解按下式计算：

$\underset{\mathrm{\&Gamma;}\&cap;{\mathrm{\&Omega;}}_i}{\&Integral;}\mathrm{\&sigma;}\mathrm{\&kappa;}nd\mathrm{\&Gamma;}\left(x_s\right)=\underset{{\mathrm{\&Omega;}}_i}{\&Integral;}\mathrm{\&sigma;}\mathrm{\&kappa;}\&dtri;{\mathrm{\&alpha;}}_{1}dV---\left(4\right)$

其中，Ω_i为第i个计算单元，Γ∩Ω_i为Ω_i中的两相流交界面部分，σ为两相流交界面张力，κ为两相流交界面的曲率系数，n为两相流交界面单位法向矢量，α₁为相1的体积分数，两相流交界面的曲率根据体积分数α₁计算：

(3)状态方程：

ρ_i＝ρ_0i+ψ_iP(5)

其中，ρ₀为标准密度[kgm^-3]，ψ为压缩性[s²m-²]，ψ＝1/c²，c为在流体中声音传播的速度[ms-¹]，根据不同相在空间的体积分数，有：

ρ＝α₁ρ₁+α₂ρ₂

μ＝α₁μ₁+α₂μ₂ (6)

U＝α₁U₁+α₂U₂

其中，α₁+α₂＝1，α₂为相2的体积分数；

(4)相连续方程：

$\frac{\&part;{\mathrm{\&alpha;}}_1}{\&part;t}+\&dtri;\&CenterDot;\left({\mathrm{\&alpha;}}_{1}U\right)+\&dtri;\&CenterDot;\left({\mathrm{\&alpha;}}_1\right(1-{\mathrm{\&alpha;}}_1\left)U_r\right)={\mathrm{\&alpha;}}_1(1-{\mathrm{\&alpha;}}_1)dgdt+{\mathrm{\&alpha;}}_1\&dtri;\&CenterDot;U---\left(7\right)$

其中，U_r为相对速度，U_r＝U₁-U₂，dgdt为压缩项：

$dgdt=({\mathrm{\&psi;}}_2/{\mathrm{\&rho;}}_2-{\mathrm{\&psi;}}_1/{\mathrm{\&rho;}}_1)\frac{DP}{Dt}.$

4.根据权利要求3所述的方法，其特征在于，所述使用有限体积法离散可压缩纳维斯托克斯方程组，通过压力修正法求解离散后的代数方程组，迭代得到液舱内可压缩流体晃荡的流速场和压力场的收敛解，具体包括如下步骤：

—在每一时刻的计算中，将上一时刻计算得到的压力场代入可压缩纳维斯托克斯方程，求出临时速度场；

—将动量方程代入连续性方程得到压力方程，并将临时速度场代入压力方程，求解出所述临时速度对应的压力场；

—求解每个单元的体积分数、密度和运动粘度系数；

—利用所述临时速度对应的压力场代入动量方程更新速度场；

—判断当前时刻，所述临时速度对应的压力场和更新后的速度场是否收敛；若不收敛，将更新后的速度场作为临时速度，返回至所述的“将动量方程带入……求解出所述临时速度对应的压力场”步骤，迭代计算；若收敛，结束当前时刻计算，进入下一时刻，重复上述计算步骤。

5.根据权利要求1所述的方法，其特征在于，所述基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，以所述合力为目标函数、所述隔板的渗水率为优化变量对所述阻晃隔板的渗水率进行最优控制，所述渗水率为所述阻晃隔板上挖空面积与阻晃隔板总面积的比值，包括：

三维液舱阻晃隔板形状优化的目标函数是两平行于隔板的垂直侧壁所感受到的合力的泛函，优化的最终目标是使两垂直侧壁所受压力合力最小，隔板上每个网格上的渗水率作为优化参数，优化约束条件是渗水率必须在0.0与1.0之间；

基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，形成隔板拓扑优化模块，对每一时刻上，隔板渗水率进行优化，进而实现对所述阻晃隔板的渗水率的最优控制。

6.根据权利要求5所述的方法，其特征在于，所述基于拓扑优化方法并结合带有自由液面的可压缩两相流运动算法，形成隔板拓扑优化模块，对每一时刻，隔板渗水率进行优化，进而实现对所述阻晃隔板的渗水率的最优控制，包括：

—设定初始渗水率η和优化开始时刻T_start；

—当激励时间t≤T_start(t＝t+Δt，Δt是时间步长，t为激励时长)，根据物理参数、初始条件和边界条件，调用带自由液面的可压缩计算流体力学模块计算流场参数，记录t＝T_start时，流场参数和初始目标函数J₀，J_min＝J₀；

—当激励时间t>T_start，在每一时刻上的计算中，首先调入优化算法，优化渗水率，用优化后的渗水率更新渗水率，再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*；

—判断当前时刻，渗水率是否最优，若J^*≥J_min，返回到所述的“当激励时间t>T_start……再调用带自由液面的可压缩计算流体力学模块计算当前时刻目标函数J^*”步骤，迭代计算，若J^*<J_min，J_min＝J^*，结束当前时刻计算，进入下一时刻(t＝t+Δt)；

—对每一时刻上的渗水率进行优化，用上一时刻优化得到的渗水率替换原有的渗水率，再重复所述的计算当前时刻目标函数步骤，直到结束。