舰艇在遭受武器攻击或发生海损事故后，很可能因为破损进水造成不沉能力的严重丧失，导致舰艇发生倾覆甚至沉没．进水平衡系统是舰艇发生破损进水后，恢复舰艇平衡的一种重要手段[1-2]．通过在左右两舷的舱室之间设置连通管，当一侧舱室破损进水时，海水通过连通管导移至另一侧舱室，从而平衡舰体[3-4]．连通管、通气管、左右舷对称舱室共同构成了进水平衡系统．对于不同的破损进水模式，破口大小和位置、连通管直径和位置、通气管直径、舱室容积的大小及船体静水力性能，共同决定了进水的平衡时间及平衡过程中不沉性指标．舰艇进水平衡是“船体-进水-气流”多个变量之间相互作用、动态耦合的复杂过程．国际海事组织《国际海上人命安全公约(SOLAS)》会议、美国海军的稳性标准(DDS079-1)和美国海岸警卫队设计与构造标准(SWBS079)均对舰艇进水的平衡时间和最大横倾角提出了要求［1，5］．针对该要求，文献[6]对舰艇进水平衡过程中横倾角随时间减小的过程开展了仿真计算；文献[7]提出了流量系数和舰艇平衡时间的解析计算方法．但是文献[6-7]以破损舱立即灌满为前提，仅研究了舰艇横倾角从最大到平衡的变化过程，忽略了破损舱水位上升与平衡舱进水同步进行的进水过程，缺乏“舱内进水、空气压缩、船体平衡”三者耦合关系的分析和建模计算．目前的研究主要集中在单个效应的建模分析，缺少进水平衡过程中“船体-进水-气流”诸变量统一的、相互迭代计算的时域微分方程组，舰艇进水平衡过程的耦合研究尚未形成有机整体．考虑到目前理论方法的局限性，当进行连通管和通气管的设计时，无法综合分析各个设计参数对进水平衡过程的整体影响，因而现有成果无法预测船体开始破损后(而非破损舱灌满水后)的进水平衡全时域过程．综上所述，须要拓展现有舰艇进水平衡过程的研究范围，以快速有效地开展不同进水平衡系统设计方案平衡能力的对比分析．1 舰艇进水平衡过程计算模型1.1　舰艇进水平衡过程耦合关系分析舰艇进水平衡过程中的耦合关系如图1和图2所示．舰艇的初始水线为WL，随着海水通过舷侧破口向破损舱涌入，破损舱进水量vw,i增加[8]，破损舱水位Hi、连通管流量Qdct、平均吃水深度Tm和10.13245/j.hust.238445.F001图1舰艇进水平衡过程示意图10.13245/j.hust.238445.F002图2舰艇进水平衡过程耦合关系图横倾角θ均会增大，破口流量Qdmg会减小，新的水线为W1L1；平衡舱进水量vw,j增大后，平衡舱水位Hj、平衡舱气压pa,j、平均吃水深度Tm均会增大，而连通管流量Qdct、平衡舱的气体体积va,j和横倾角θ均会减小；平衡舱气压pa,j增大后，通气管气流质量流量m˙a,j会增大，连通管流量Qdct会减小；平均吃水深度Tm和横倾角θ增大，均会造成破口流量Qdmg增大．根据图2所示的耦合关系，主要须要建立水流流量计算模型、空气压缩效应计算模型和“船-水-气”动态耦合关系模型．1.2　水流流量计算模型a. 破口流量计算模型基于伯努利方程，可得舷侧破口的流量方程为[9]Qdmg=Cd,iAdmg2gΔH，(1)式中：Cd,i为破口流量系数(依据文献[9]设定为0.6）；Admg为破口横截面积；∆H为破口两侧液面高度差，与破损舱水位Hi、船舯吃水深度Tm及横倾角θ均相关；g为重力加速度．由于破口流量Qdmg会随着时间变化，因此可将因变量破口流量Qdmg与自变量破损舱水位Hi、船舯吃水深度Tm、横倾角θ、破口横截面积Admg、破口高度Hdmg和时间t之间的耦合关系记为Qdmg(t)=F1(Hi(t),Tm(t),θ(t),t,Admg,Hdmg), (2)式中：变量Hi，Tm，θ会随着时间t的变化而变化；变量Admg和Hdmg在进水仿真前设置完成．b. 连通管水流计算若破损舱内水流的速度可忽略不计，基于伯努利方程，则通过横截面积为Adct的连通管体积流量为[10]Qdct(t)=Cd,jAdct2gΔH-2pa,j/ρw，(3)式中：Cd，j为连通管流量系数；ΔH为连通管两侧的水位差；pa，j为平衡舱气压．由于连通管流量Qdct会随时间变化，因此可将因变量连通管流量Qdct与自变量破损舱水位Hi、平衡舱水位Hj、平衡舱气压pa，j、连通管流量系数Cd，j、连通管横截面积Adct、连通管中心点高度Hdct和时间t之间的耦合关系记为    Qdct(t)=F2(Hi(t),Hj(t),pa,j(t),Cd,j,Adct,Hdct,t), (4)式中：变量Hi，Hj，pa，j会随着时间t的变化而变化；变量Cd，j和Adct在进水仿真前设置完成．c. 舱室水位上升微分方程破损舱水位上升的时域微分方程为：dvw,i(t)/dt=Qdmg(t)-Qdct(t);Si(t)dHi(t)/dt=Qdmg(t)-Qdct(t), (5)式中vw,i(t)和Si(t)分别为破损舱t时刻的舱室进水量和自由液面面积．平衡舱水位上升的时域微分方程为：dvw,j(t)/dt=Qdct(t);Sj(t)dHj(t)/dt=Qdct(t), (6)式中vw,j(t)和Sj(t)分别为平衡舱t时刻的舱室进水量和自由液面面积．舱室自由液面面积与舱室进水量和舱室型线有关，因而将破损舱自由液面面积记为关于进水量的隐函数Si=F3(vw,i)，(7)将平衡舱自由液面面积记为关于进水量的隐函数Sj=F4(vw,j)．(8)1.3　空气压缩效应计算模型a. 平衡舱总体积守恒微分方程舱室总体积守恒方程将进水与空气占据的舱室体积之和不变作为耦合计算的基本条件，由于舱室减少的气体体积等于舱室流入的液体体积，对于平衡舱而言[11]，dva,j(t)/dt=-dvw,j(t)/dt, (9)式中va，j为平衡舱的气体体积．b. 舱室出气量平衡微分方程随着海水向平衡舱内漫延，舱内气体体积受到压缩，舱室气压升高．当平衡舱内气压高于大气压时，平衡舱内空气会溢出至舱外，最终平衡舱气压与大气压相等．根据波义尔定律，可得平衡舱的出气量平衡微分方程[12]va,jρatmpatmdpa,jdt+(pa,j+patm)ρatmpatmdva,jdt=-m˙a,j，(10)式中：m˙a,j为平衡舱的通气管质量流量；pa，j为平衡舱的气压；patm为大气压；ρatm为大气密度．c. 通气管气流计算模型由于忽略舱内气流的速度及气流势能的变化，根据伯努利方程，可得平衡舱通气管的流量[11]    m˙a,j=Cd,aAa{2ρatmpatm∙[ln(pa,j+patm)-ln patm]}1/2, (11)式中：Aa为平衡舱通气管的有效横截面积；Cd，a为平衡舱通气管流量系数．由于破损舱的气密性完全丧失，因此可假设其处于完全通风状态．1.4　“船-水-气”动态耦合关系模型对舰艇进水平衡过程中的微分方程组求解，联立式(2)和式(4)～(8)，可得在进水平衡过程中的进水时域微分方程组：Qdmg(t)=F1(Hi,Tm,θ,Admg,Hdmg,t);Qdct(t)=F2(Hi,Hj,pa,j,Cd,j,Adct,Hdct,t);dvw,i(t)/dt=Qdmg(t)-Qdct(t);Si(t)dHi(t)/dt=Qdmg(t)-Qdct(t);Si=F3(vw,i);dvw,j(t)/dt=Qdct(t);Sj(t)dHj(t)/dt=Qdct(t);Sj=F4(vw,i). (12)联立式(6)和式(9)～(11)，可得平衡舱在进水过程中气液流动的时域微分方程组：dva,j(t)/dt=-Qdct(t);va,jρatmpatmdpa,jdt+(pa,j+patm)ρatmpatmdva,jdt=-m˙a,j;m˙a,j=Cd,aAa{2ρatmpatm[ln(pa,j+patm)-lnpatm]}1/2. (13)由于舰艇的平均吃水深度Tm、横倾角θ和纵倾角ϕ确定，并且破损舱进水量vw，i和平衡舱进水量vw，j确定后，破损舰艇的水下容积V、浮心坐标(XC，YC，ZC)和重心坐标(XG，YG，ZG)均为定值，因此可得舰艇进水平衡过程中的平衡方程式[13]：V(Tm,θ,ϕ,vw,i,vw,j,t)-Vini=0;    Mzx(Tm,θ,ϕ,vw,i,vw,j,t)+Mxy(Tm,θ,ϕ,vw,i,vw,j,t)tan θ=0;    Myz(Tm,θ,ϕ,vw,i,vw,j,t)+Mxy(Tm,θ,ϕ,vw,i,vw,j,t)tan ϕ=0. (14)式(12)～(14)中共包括时间t和14个随时间t变化的变量(破口流量Qdmg、连通管流量Qdct、破损舱进水量vw，i、破损舱水位Hi、破损舱自由液面Si、平衡舱进水量vw，j、平衡舱水位Hj、平衡舱自由液面Sj、平衡舱气压pa，j、平衡舱空气体积va，j、平衡舱通气管质量流量m˙a,j、平均吃水深度Tm、横倾角θ、纵倾角ϕ)，然而式(12)和(13)中共有14个方程．除时间变量t以外，每个变量在耦合方程组中均既为自变量，又为因变量，构成了相互耦合的关系．因此，必须将时间t离散化，通过设定时间仿真步长Δt，计算特定时刻下各个变量的值．1.5　时域耦合方程组求解流程针对实际舰艇进水平衡过程较为缓慢的特点，基于准静态假设的舰艇进水平衡过程时域耦合方程组求解流程如图3所示．先确定进水过程初始条件和结束条件，然后求解破口流量与连通管流量，并求解平衡舱的通气管质量流量．通过基于四阶龙格-库塔法求解式(12)和(13)中的微分方程，可以得出任意时刻t时，破损舱和平衡舱的水位高和进水量，以及平衡舱的气压和气体体积，进而求解舰艇新的自由液面面积和舰艇浮态与稳性参数，并将此作为下一仿真时刻的准静态．如此循环计算，直至破口流量与连通管流量为零．10.13245/j.hust.238445.F003图3时域耦合方程组求解流程当舰艇破损进水达到平衡时，舰艇姿态不再变化，因此可以将吃水深度和横倾角是否发生变化作为仿真是否结束的判别条件，即Tmk+1-Tmk≤ε且θk+1-θk≤δ，(15)式中：Tmk和Tmk+1分别为第k个和第k+1个时刻的平均吃水深度；θk和θk+1分别为第k个和第k+1个时刻的横倾角；ε和δ均为足够小的数．假设将舰艇在任意时刻都看作是准静止状态，也就是说对于任意时刻t，Tm和θ可以认为是不变的．当时刻t固定时，Tm，θ和ϕ可以看作是定值．此时，微分方程组(12)和(13)就只包含了11个随时间t变化的变量并且由11个非线性微分方程构成，具备了求解的基本要求．2 求解算法的实验验证2.1　船模实验设计方案2.1.1　船模舱室划分与不沉性参数船模船体与舱室的形状尺寸和实物如图4所示．该船模由7个舱室组成，钢板的厚度均为2 mm，其中3号舱作为破损舱，4号舱作为平衡舱．4号舱顶部用透明可拆卸的有机玻璃密封好，并可安装通气管．船模典型设计参数(不沉性设计指标)如下：初始排水体积Vini=0.233 m3，水线面面积S=0.6 m2，船模长度LS=1.5 m，横倾角初值θini=0°，船模型深DS=0.8 m，船模型宽BS=0.4 m，吃水深度初值Tini=0.426 m，初稳度初值hini=0.023 m．10.13245/j.hust.238445.F004图4船模中横剖面视图和俯视图(m)2.1.2　通气管、连通管及其流量系数计算可拆卸式弯头形状的通气管外形如图5所示，3号舱与4号舱之间设有3个直线式连通管．3个通气管和3个连通管的尺寸规格与根据文献[14]计算得出的流量系数理论值计算结果见表1，表中：D为内径；tn为壁厚；L为长度；Cd,j为流量系数．10.13245/j.hust.238445.F005图5通气管(mm)10.13245/j.hust.238445.T001表1通气管与连通管流量系数理论值管D/mmtn/mmL/mmCd,j1号通气管72.04000.4852号通气管101.54000.4943号通气管121.54000.4991号连通管302.01000.6332号连通管602.01000.6253号连通管802.01000.6172.1.3　船模进水过程实验数据的采集船模艏、艉、左舷和右舷及舱室底部均安装有水位传感器，水位传感器量程为1 m，精度为5 mm．2.2　进水过程实验案例3号舱破口面积Admg=36 cm2，破口中心高度Hdmg=0.1 m．4号舱顶部安装有10 mm通气管，2号连通管(直径为60 mm，中心点高度Hdct=130 mm)将3号舱内的进水漫延至4号舱，1号连通管和3号连通管封闭．进水漫延过程中3号舱水位变化曲线如图6所示，4号舱水位变化曲线如图7所示，船模横倾角变化曲线如图8所示．由于纵倾角ϕ的变化可忽略不计，因此仅分析舰艇横倾的变化；同时，船模吃水深度由0.43 m缓慢上升至0.52 m，船模进水过程共有以下5个进水阶段．10.13245/j.hust.238445.F006图6破损舱(3号舱)水位变化曲线a. 当t＜1.0 s时，3号舱水位和船模横倾角迅速升高．由于该阶段3号舱水位未淹没连通管，连通管流量为0，因此4号舱进水量为0．10.13245/j.hust.238445.F007图7平衡舱(4号舱)水位变化曲线10.13245/j.hust.238445.F008图8船模横倾角变化曲线b. 当1.0 s≤t＜2.8 s时，3号舱水位已淹没连通管和舷侧破口．随着3号舱水位不断升高，破口流量不断减小，连通管流量不断增大，造成3号舱水位上升趋慢，4号舱水位上升趋快．c. 当t=2.8 s时，3号舱进水速率与4号舱相等，船模的横倾角达到最大值7.3°．当2.8 s＜t＜8.6 s时，船模开始逐步趋于平衡．随着4号舱水位上升和4号舱气压产生反压，连通管流量开始减小，4号舱水位上升趋慢．d. 当8.6 s≤ t ≤27.7 s时，4号舱水位上升速度进一步变慢，导致船模横倾平衡速率明显趋缓．e. 当t＞27.7 s时，3号舱水位、4号舱水位与平均吃水深度均为0.52 m，进水平衡过程结束，船模完全恢复平衡．根据图6~8的模拟值与实验值的对比结果可知：由于破口尺寸与连通管尺寸适中，因此通气管通过增大舱室气压能够起到限制水流的作用，水流较为平稳，符合仿真准静态的条件假设．仿真数据与实验数据的变化趋势完全相同，数值大小的符合度较高，能够验证水流流量计算模型、空气压缩效应计算模型、破损浮态与稳性计算模型之间耦合关系的正确性．当1.0 s≤t＜2.8 s时，由于4号舱开始迅速进水，使得4号舱内空气急剧压缩，出现4号舱实际进水量比理论计算量小的情况，因此存在对船模横倾角估计偏低的现象．3 案例对比分析设置3号舱破口面积Admg=10 cm2，进行多组案例对比计算，结果见表2~4，归纳得出以下规律．10.13245/j.hust.238445.T002表230 mm直径连通管仿真结果参数7 mm通气管10 mm通气管12 mm通气管完全通风(Aa/Adct)/%5.411.116.0—平衡时间/s60.244.342.140.5最大横倾角/(°)11.39.99.69.410.13245/j.hust.238445.T003表360 mm直径连通管仿真结果参数7 mm通气管10 mm通气管12 mm通气管完全通风(Aa/Adct)/%1.42.84.0—平衡时间/s47.825.522.119.2最大横倾角/(°)9.86.24.93.410.13245/j.hust.238445.T004表480 mm直径连通管仿真结果参数7 mm通气管10 mm通气管12 mm通气管完全通风(Aa/Adct)/%0.81.62.3—平衡时间/s46.323.720.117.2最大横倾角/(°)9.75.84.12.7a. 舱室通气管内径对平衡能力的影响随着通气管内径的增大，舰艇平衡时间缩短，最大横倾角变小，该规律表明增大通气管有利于增强舰艇的平衡能力．根据表2的计算结果，当通气管面积与连通管面积的比值超过10%时，空气压缩效应可忽略．b. 进水平衡系统设计方案优选舰艇设计人员须要在舰艇平衡能力与占用总体资源之间取最优解．在表2~4的仿真计算结果中，10 mm通气管的平衡能力显著优于7 mm通气管，12 mm通气管的平衡能力略强于10 mm通气管，因而在3种通气管中10 mm通气管为最优解．在表2～4的仿真计算结果中，60 mm连通管的平衡能力显著优于30 mm连通管，80 mm连通管的平衡能力略强于60 mm连通管，因而在3种连通管中60 mm连通管为最优解．4 结语为解决舰艇进水平衡过程的时域耦合建模问题，首先分析了船体舱室、进水和气流之间的动态耦合效应，然后构建了时域耦合微分方程组，最后通过船模实验验证了该时域耦合模型的准确性，并揭示了进水平衡过程经历的五个阶段．仿真研究发现：随着连通管与通气管横截面积增大，其平衡能力的边际效应更显著，即最大横倾角的减小和舰艇平衡时间的缩短不明显；因此，在舰艇设计过程中，应当权衡进水平衡系统占用空间资源与其具备的平衡能力，从而实现连通管和通气管的优化设计．与此同时，舰艇进水平衡过程时域耦合模型的计算结果也能够为开展舰艇平衡态势预判和舰艇平衡指挥决策提供快速精准的支撑和依据．