气泡包络仿真方法和装置及有仿真程序的可读存储介质 本发明涉及一种用于通过用气泡包络船舶来减少表面摩擦的数字仿真的方法和装置，以及具有记录的仿真程序的计算机可读的存储介质，尤其涉及一种通过在边界层中射入微气泡，以便包络船的浸在水中的表面来减少摩擦的技术。
关于减少表面摩擦的技术在日本专利申请首次公开S50-83992,S53-136289,S60-139586,S61-71290，以及实用新型申请S61-39691,S61-128185中披露了。在这些公开文本中所述的减少摩擦的船基于从航行的船向水中射入气体例如空气，以便在船和水之间的边界层中产生微气泡(下面称为气泡)，利用这些包络浸在水中的表面地气泡的作用，减少表面摩擦。
一种已知的仿真技术是计算流体动力学(CFD)方法。该方法根据解三维Navier-Stokes方程，该方程通过数值分析直接地或近似地描述各种粘滞流场。通过把包括紊流边界层的流场划分为许多正交的三维格子(方格)，使得在每个方格中可以解三维Navier-Stokes方程，从而计算每个方格中的紊流能量和能量耗散系数(由于摩擦的热而损失的能量的部分)，对在船表面上的紊流进行数值仿真(紊流仿真)。这种船表面附近的紊流场的数值仿真研究使得能够利用最少量的计算设计减少摩擦的船。
然而，现有的仿真方法被设计用于提供船表面附近的定性的流场图形，因而，不能够提供关于存在于船和水之间的边界层中的气泡分布的定量的数据，即不能提供说明气泡在实际上是如何包络船的数据。因此，在理论上，难于评价关于表面摩擦效应的各种技术，这样，使得研究者除去实验方法之外别无选择。这种减少摩擦效应的实验研究在花费的时间和劳动方面是极其昂贵的。尤其是，迄今为止不能进行由于气泡的存在而产生的减少摩擦效应的理论计算，不能提供为驱动减少摩擦的船所需的功率的理论估算。
本发明提供一种用于解决许多上述问题的方法，其几个目的如下：
a提供一种方法和装置，用于通过按照一种应用程序计算由从航行中的船射入气体而产生的气泡的分布数据，仿真包络船的气泡，并提供一种包括仿真程序的计算机可读的存储介质。
b提供一种用于仿真包络船的气泡的方法和装置，使得能够按照一种应用程序进行减少表面摩擦技术的理论评价，并提供一种包括仿真程序的计算机可读的存储介质。
c提供一种用于仿真包络船的气泡的方法和装置，使得能够以低的费用按照一种应用程序进行减少表面摩擦技术的理论评价，并提供一种包括仿真程序的计算机可读的存储介质。
d提供一种用于按照一种应用程序仿真包络任何形状的船的气泡的方法和装置，并提供一种包括仿真程序的计算机可读的存储介质。
因此，一种仿真方法包括以下步骤：
a把包括在船的水下表面周围形成的紊流边界层的流场划分成三维正交的网格而形成方格，并计算每个方格中的紊流能量、能量耗散系数和平均流体速度；
b在给定的时刻从给定的射流喷口以给定的初始气泡速度向紊流层中射入气泡；
c按照紊流能量、能量耗散系数、平均流体速度和随机数计算每个方格中的定向气泡速度；
d按照定向气泡速度和平均流体速度通过解运动方程计算在相继的单位时间的气泡的位置；
e按照在相继的单位时间的气泡位置计算每个方格中的空隙比；
f对于给定的迭代次数重复步骤c到步骤e，直到通过使时间参数增加一个单位时间而达到一个整数时间；以及
g计算所有方格内的空隙比的分布，以便在所述整数时间表示所有的气泡分布，并且当确定气泡分布收敛时结束计算。
一种实施所述仿真方法的装置，其通过执行以下的分析步骤构成：
a把包括在船的水下表面周围形成的紊流边界层的流场划分成三维正交的网格而形成方格，并计算每个方格中的紊流能量、能量耗散系数和平均流体速度；
b在给定的时刻从给定的射流喷口以给定的初始气泡速度向紊流层中射入气泡；
c按照紊流能量、能量耗散系数、平均流体速度和随机数计算每个方格中的定向气泡速度；
d按照定向气泡速度和平均流体速度通过解运动方程计算在相继的单位时间的气泡的位置；
e按照在相继的单位时间的气泡位置计算每个方格中的空隙比；
f对于给定的迭代次数重复步骤c到步骤e，直到通过使时间参数增加一个单位时间而达到一个整数时间；以及
g计算所有方格内的空隙比的分布，以便在所述整数时间表示所有的气泡分布，并且当确定气泡分布收敛时结束计算。
一种存储介质，其包括执行以下步骤的应用程序：
a把包括在船的水下表面周围形成的紊流边界层的流场划分成三维正交的网格而形成方格，并计算每个方格中的紊流能量、能量耗散系数和平均流体速度；
b在给定的时刻从给定的射流喷口以给定的初始气泡速度向紊流层中射入气泡；
c按照紊流能量、能量耗散系数、平均流体速度和随机数计算每个方格中的定向气泡速度；
d按照定向气泡速度和平均流体速度通过解运动方程计算在相继的单位时间的气泡的位置；
e按照在相继的单位时间的气泡位置计算每个方格中的空隙比；
f对于给定的迭代次数重复步骤c到步骤e，直到通过使时间参数增加一个单位时间而达到一个整数时间；以及
g计算所有方格内的空隙比的分布，以便在所述整数时间表示所有的气泡分布，并且当确定气泡分布收敛时结束计算。
上面概括的方法能够利用定量的数据仿真包络船的气泡，这用现有的方法是不可能的。因而，现在能够在理论上评价各种减少表面摩擦的技术，从而以高的效率和低的费用实现减少摩擦的船。采用κ-ε模型(κ：紊流能量，ε：能量耗散系数)，能够使用CFD数据仿真水下的船体表面上的流场。
图1表示在基本理论中使用的坐标系统；
图2表示在计算时使用的气泡分布；
图3是方法的流程图；
图4是表示在低的迭代次数之后空隙比分布的收敛状态的曲线；
图5是表示收敛的空隙比分布的曲线；
图6是对于气泡直径为500μm的第一分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图7是对于气泡直径为500μm的第二分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图8是对于气泡直径为500μm的第三分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图9是对于气泡直径为500μm的第四分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图10是表示用气泡包络的船的计算结果的曲线；
图11是对于气泡直径为1000μm的第一分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图12是对于气泡直径为1000μm的第二分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图13是对于气泡直径为1000μm的第三分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
图14是对于气泡直径为1000μm的第四分布曲线，表示在船的围长方向在各个纵向位置空隙比分布的计算结果；
首先说明在本发明的仿真方法中使用的基本理论。
本方法的根据是基于在船的船体表面形成的给定的方格内的流场数据的κ-ε模型的CFD分析，其中κ是紊流能量，ε是能量耗散系数。κ和ε的导出的值被用于下面所示的动力学方程中，以便计算存在于围绕水下表面的各个方格中的气泡的数量，借以得到船体表面上的全部气泡分布。
1．证明应用Monte Carlo方法的合理性
在本发明中，采用根据使用随机数的直接仿真方法的MonteCarlo方法的形式作为在包括紊流边界层的流场中的动力学气泡的分散方程的解。首先说明使用Monte Carlo方法解分散方程的合理性。
假定在单位时间内气泡从一个方格随机地运动到一个方格。在给定时刻在点(i,j,k)的气泡在经过一个单位时间之后将运动到以下6个点之一：
点①(i-1,j,k)
点②(i+1,j,k)
点③(i,j-1,k)
点④(i,j+1,k)
点⑤(i,j,k-1)
点⑥(i,j,k+1)
如果假定在单位时间内气泡将运动到另一个方格而与过去的路径历史无关，则从原点(i=0,j=0,k=0)开始的气泡在n个单位时间之后到达点(i,j,k)的几率可以认为是在该给定时刻n的浓度值。对于在n+1个单位时间的时刻到达点(i,j,k)的气泡，该浓度被表示为Pni,j,k，在n个单位时间的时刻气泡必然位于上述6个点中的一个上。
气泡从6个点中的一个到达点(i,j,k)的几率是1/6，因此，可以建立以下的递归公式：Pn+1i,j,k=(Pni,j,k+1+Pni,j,k-1+Pni,j+1,k+Pni,j-1,k+Pni+1,j,k+Pni-1j,k)/6 (1)
式(1)和在较早的一个单位时间的浓度Pni,j,k之间的差，即式(2)被得到：pn+1i,j,k-Pni,j,k={(Pni,j,k+1-Pni,j,k)+(pni,j,k-1-Pni,j,k)}/6+{(Pni,j+1,k-Pni,j,k)+(Pni,j-1,k-Pni,j,k)}/6+{(pni+1,j,k-Pni,j,k)+(Pni-1,j,k-Pni,j,k)}/6 (2)
如果气泡速度为0，则对于可变的气泡浓度P和分散系数D，在各向同性的紊流中分散方程可以认为如式(3)：∂P∂t=D(∂2P∂x2+∂2P∂y2+∂2P∂z2)---(3)]]>
当式(3)被重写为在时刻n的浓度Pni,j,k和在时刻n+1的浓度Pn+ 1i,j,k之间的差分方程时，得到式(4)：Pn+1i,j,k-Pni,j,kΔt]]>=DΔx{(Pni+1,j,k-Pni,j,kΔx)-(Pni,j,k-Pni-1,j,kΔx)}]]>+DΔy{(Pni,j+1,k-Pni,j,kΔy)-(Pni,j,k-Pni,j-1,kΔy)}]]>+DΔz{(Pni,j,k+1-Pni,j,kΔz)-(Pni,j,k-Pni,j,k-1Δz)}---(4)]]>其中，如果紊流是各向同性的，由分散系数的定义，建立方程(5)D≡<u·1> =Δx·Δx/6Δt =Δy·Δy/6Δt =Δz·Δz/6Δt (5)
其中括号<>表示在分散过程中涉及的所有气泡的集体平均(集平均)，“u”是与气泡的平均气泡速度的偏差，“1”是位移偏差。
把式(5)代入式(4)得到式(6)。式(6)和式(2)匹配描述气泡的随机运动。因此，可以使用采用随机数直接仿真气泡运动的Monte Carlo方法解各向同性紊流中的分散方程。Pn+1i,j,k-Pni,j,k=16{(Pni+1,j,k-Pni,j,k)-(Pni,j,k-Pni-1,j,k)}]]>+16{(Pni,j+1,k-Pni,j,k)-(Pni,j,k-Pni,j-1,k)}]]>+16{(Pni,j,k+1-Pni,j,k)-(Pni,j,k-Pni,j,k-1)}----(6)]]>
2．在方格内气泡的运动和空隙比
在本部分将探讨通过利用κ-ε表示紊流而使用CFD数据计算流场。假定在方格内的紊流数据是均匀的和各向同性的，并且在船周围的边界层内的整个流场可以用在各个方格内的流场的集描述。在本研究中，使用在方格内紊流是均匀的和各向同性的条件用以下的动力学方程表示气泡的运动。
图1用于描述气泡运动方程的基本坐标系统。这XYZ-0系统是一个连续的坐标空间，Z轴是垂直方向(plumb verticaldirection)。速度参数ux,uy,uz分别是沿X-,Y-,Z方向的平均流体速度(在本例中为围绕船的水的速度)。下面使用这一坐标系统说明：(1)当气泡远离船的表面时的运动方程；(2)当气泡黏附在船的表面时的运动方程；(3)紊流分散的影响；以及(4)空隙比的定义。
(1)当气泡远离船时的运动方程
在一个方格内，假定气泡流数据是均匀的和各向同性的，但是在邻近的方格内，速度梯度和压力梯度对气泡运动的影响将被考虑。考虑这种影响的一种方法由Meng J.C.S和Uhlman J.S.Jr.在“Micro-bubble formulation and splitting in a turbulentboundary layer for Turbulence Reduction.A Symposium in Honorof Maurice Holt on His 70th Birthday,Willamsburg,Va.pp,168-217”,1989中描述了。
按照这种方法，速度梯度的影响由作用在实心球上的升力L(Saffman升力)近似，在本研究中实心球用气泡代替。式(7)根据液体密度p(kg/m3)，流体速度U(m/s)，气泡半径R(m)，运动黏度系数ν(m2/s)表示升力L。关于Saffman升力的细节，请参阅Saffman的“The lift on a small sphere in a slow shearflow,JFM,vol.22”。L=6.46ρ·U·R2(ν∂U∂n)1/2----(7)]]>
压力梯度的影响由式(8)近似，其中Fi是沿i方向由压力梯度施加的力；V是一个气泡的容积(m3)；Pr是时间平均压力；xi是沿i方向的位移。Fi=-V∂Pr∂xi----(8)]]>
因而，在本坐标系统中沿X-Y-Z方向的气泡的运动方程由式(9)-(11)表示，其中ux,uy,uz是沿X-Y-Z方向的水的平均速度(m/s),ux’,uy’,uz’是定向气泡流的平均速度，m是气泡的质量(kg)，mA是附加质量(kg),μ是静态黏度系数(m2/s),g是重力加速度(m/s2)。(m+mA)X··=6πμ(ux+ux′-X·)R]]>+6.46ρR2ν1/2{(uy+uy′-Y·)(∂uy∂x)1/2+(uz+uz′-Z·)(∂u2∂x1/2}-V∂Pr∂x----(9)]]>(m+mA)Y··=6πμ(uy+uy′-Y·)R]]>+6.46ρR2ν1/2{(uz+uz′-Z·)(∂uz∂y)1/2+(ux+ux′-X·)(∂ux∂y)1/2}-V∂Pr∂y---(10)]]>(m+mA)Z··=6πμ(uz+uz′-Z·)R-ρgV+mg]]>+6.46ρR2ν1/2{(ux+ux′-X·)(∂ux∂z)1/2+(uy+uy′-Y·)(∂uy∂z)1/2}-V∂Pr∂z----(11)]]>
换句话说，在三个坐标方向中的气泡的运动方程通过把定向紊流速度表示为平均流体速度和定向气泡流速度的和的函数被建立了。在这些方程(9)-(11)中，方程右边的第一项通过使用用于实心球的Stokes近似进行选择。
(2)当气泡黏附在船的表面时的运动方程；
当气泡黏附在船上时，在气泡上的反作用力是有效的垂直力和摩擦力。因此，可从方程(9)-(11)得到运动方程，如式(12)-(14)所示。
在方程(12)-(14)中，NX,NY,NZ是有效垂直力(单位牛顿)，FX,FY,FZ是摩擦力(单位牛顿)。(m+mA)X··=6πμ(ux+ux′-X·)R]]>+6.46ρR2ν1/2{(uy+uy′-Y·)(∂uy∂x)1/2+(uz+uz′-Z·)(∂uz∂x)1/2}-V∂Pr∂x-Nx+Fx----(12)]]>(m+mA)Y··=6πμ(uy+uy′-Y·)R]]>+6.46ρR2ν1/2{(uz+uz′-Z·)(∂uz∂y)1/2+(ux+ux′-X·)(∂ux∂y)1/2}-V∂Pr∂y+Ny+Fy----(13)]]>(m+mA)Z··=6πμ(uz+uz′-Z·)R-ρgV+mg]]>+6.46ρR2ν1/2{(ux+ux′-X·)(∂ux∂z)1/2+(uy+uy′-Y·)(∂uy∂y)1/2}-V∂Pr∂z+Nz+Fz----(14)]]>
在这一研究中，为了简化，在每个轴向的摩擦力FX,FY,FZ取为0。用于获得有效垂直力NX,NY,NZ的约束条件假定通过把船表面坐标XS,YS,ZS表示为下式(15)和(16)而被给出。dYsdXs=Y·X·---(15)]]>dZsdYs=Z·Y·---(16)]]>
(3)紊流分散的影响
在各向同性流下紊流分散对气泡运动的影响通过使用随机数随机地改变定向气泡速度ux’,uy’,uz’对气泡轨迹施加扰动而被包括在上面概述的Monte Carlo方法中。例如，定向气泡速度ux’,uy’,uz’可以用扰动的幅值u’(绝对值)和运动角φ,θ表示，如下式(17)-(19)所示。u’x= u’sinθ (17)u’y= u’cosθcosφ (18)u’z= u’cosθsinφ (19)
其中，假定扰动的幅值u’跟随一阶马尔可夫过程，并且由[-180°和180°]之间的均匀的随机数给出在-180°和180°之间的运动角φ，由[-90°和90°]之间的均匀的随机数给出在-90°和90°之间的运动角θ。马尔可夫过程是这样一种过程，其中在某一时刻一个事件的几率受前一个单位时间的过去几率的影响，一阶马尔可夫过程通过在前一个单位时间的随机现象的一阶方程描述在某一时刻的随机现象。
那就是说，在(n+1)单位时间的扰动u’n+1的幅值由在n单位时间的扰动幅值的一阶方程给出，如下式(20)所示。
u′n+1=χ·u′n+γn (20)
其中，χ是0＜χ＜1中的一个常数，γn是0平均值和标准偏差σ0的正态随机数。由式(20)获得的扰动的幅值u’n的自相关函数R(m,Δt)如式(21)所示。R(m,Δt)=un′·un+m′‾(u′n)2]]>=χm(un′)2‾+χm+1(un′·γn)...‾.+χ(un′·γn+m-1)‾(un′)2----(21)]]>=χm
其中，让我们考虑包括积分时标T*，其定义为衰减时间常数或在紊流边界层中的扰动的寿命。积分时标T*的定义如式(22)所示，其中Δt是时间增量，τ是微分时间。T*=∫0∞R(τ)dτ=∫0∞χτΔtdτ=-Δtlogχ---(22)]]>因而，常数χ可以用时间积分T*表示如下。χ=exp(-Δt/T*) (23)
如果m=1,Δt=T(连续的间隔)，则常数χ成为连续系统中的自相关函数，并且等效于泰勒分散理论中使用的拉格朗日自相关函数。众所周知，拉格朗日自相关函数和实验结果吻合很好。这些考虑使得我们相信，用马尔可夫过程表示在气泡流场中的扰动是正确的。
从式(20)获得的幅值u’n的扰动分散在式(24)中给出。σu2=χ2·σu2+σ02 (24)
其中σu2是定向气泡速度的方均值，其可以表示为紊流能量κ，如式(25)所示。σu2=2κ (25)由式(24)和式(25)，给出标准偏差如下。σ0=(1-χ2)1/2·σu (26)
继续研究积分时标T*的表达式，试图根据能量耗散系数ε(m2/s3)和紊流能量κ表示T*。为此，参考Tomomasa Tatsumi编的“Scienceof Turbul ent Flow”，东京大学出版社，第372页，其中表明能量耗散系数ε可以根据气泡流速u的流向x和扰动分量u’(m/s)表示，如式(27)所示。ϵ=15ν(∂u′∂x)2‾----(27)]]>
泰勒标度λ(m)被定义为流速u的偏差分量(deviationcomponent)的曲率半径，由式(28)给出。λ2=u′2‾(∂u′∂x)2‾----(28)]]>紊流能量κ由式(29)给出。κ≡12u′2‾----(29)]]>由式(28)和(29)，建立以下方程。λ2=30ν·κ/ε (30)
在另一方面，按照上述的参考文献“Science of TurbulentFlow”，泰勒标度λ和积分时标可通过下式(31)相关。
T*=λ2/10ν (31)因此，由式(30)和(31)，积分时标T*由下式给出。
T*=3κ/ε (32)
(4)空隙比
下面说明空隙比的定义。在按照上述方程计算气泡运动方程之后，通过用在给定时刻在方格内存在的气泡的容积vi除以方格的容积Δv，可以得到空隙比。例如，在图2所示的情况下，在一个方格内有10个气泡，因而由式(33)得到空隙比
(5)仿真方法的说明
下面参照图3所示的流程图说明仿真方法的一个例子。所述流程图表明被记录在计算机的仿真程序中的各个步骤。
仿真根据使用已知的κ-ε紊流模型利用CFD数据仿真船体表面上的流场数据。具体地说，根据船的尺寸和形状以及航行速度等，使用模型预先计算基本流场参数。基本流场参数包括：包括紊流边界层的流场被划分成的合适的方格数；定向流体流速ux,uy,uz的平均值；时间平均压力Pr，紊流能量κ和能量耗散系数ε。在下面的仿真方法的例子中，船的尺寸是L×B×d=7×0.6×0.1m，航行速度为3m/s。气泡射流出口是矩形的，宽度为0.3m，从船头沿着围长方向以间隔0.12m设置一排25个射流出口。
通过把由仿真模型获得的能量耗散系数ε和紊流能量κ代入式(32)，计算得到积分时标T*为0.1-2.9s。单位时间(计算增量Δt)被这样选择，使其比最小的积分时标(0.1s)足够小。在这种情况下，取1/10，从而得到0.01s作为单位时间。
计算由初始化气泡的位置和初始气泡速度开始(步S1)。例如，在船的围长方向多个射流出口(本例中为25个)由方格结构的交点表示。选择这些交点中的一个(例如端点)作为气泡的初始注射。对于气泡位置参数“Pr”，时间参数“n”两者都被设为〖0〗，但初始速度可以设为0或某个其它的值。
接着，假定直径为500μm的第一气泡被射入水中(边界层)，并按照产生的一个随机数计算对于第一气泡的扰动(步S2)。具体地说，由每个随机数发生器产生在〖-180°,180°〗之间的均匀的随机数作为运动角φ，在〖-90°,90°〗之间产生均匀的随机数作为运动角θ，并把所得的值代入式(17)-(19)中的φ,Q。
另一个随机数发生器产生具有0平均和标准偏差σ0的正态随机数，并被输入式(20)的γn中。积分时标T*和根据积分时标T*计算的单位时间Δt被输入式(23)中，以便计算常数χ。所算得常数χ和从式(29)获得的过去时间(一个单位时间之前)的扰动u’n的幅值一道被输入式(20)中，从而产生在当前时间u’n+1的扰动的幅值。这样按照由各个算法产生的随机数便得到定向气泡速度u’x,u’y,u’z。
继续进行该过程，检查在步S1射入水中的气泡是否已经离开船的表面(步S3)。定向气泡速度u’x,u’y,u’z的计算结果表明现在气泡的位置，即现在气泡占据哪个方格。
在步S3，如果确定气泡已经离开船的表面，则使用运动方程(9)-(11)计算下一个时刻(相继的单位时间)气泡的位置(方格位置)。
在另一方面，在步S3，如果确定气泡没有离开船的表面，即仍然黏附于表面上，则进行在坐标方向NX,NY,NZ上有效(effective)垂直力的初始化(步S5)。然后，在下一时刻的气泡位置(方格位置)按照运动方程(12)-(14)进行计算(步S6)，然后检查由运动方程(12)-(14)计算的结果是否满足约束条件方程(15)、(16)(步S7)。
如果判定约束条件不满足，则在步S5中的有效的垂直力NX,NY,NZ的初始值被校正(S8)，继续进行步S6、S7的迭代处理，直到计算结果满足约束条件方程(15)、(16)。
当通过一系列的迭代步骤S6-S8使得满足约束条件方程(15)、(16)时，步S7的决定成为〖yes〗，过程前进到检查有效的垂直力NX,NY,NZ是否小于0(步S9)。如果在步S9判断结果为〖yes〗，则该数据流被送到步S4，如果判断结果为〖no〗，则进行步S10的处理。
来自步S9中的〖no〗的数据流和步S4中的计算结果被送到步S10，检查气泡是否在静止水面以下，即气泡是否在水中。如果在步S10的确定结果是〖yes〗，则按照式(33)计算每个方格中的空隙比(步S11)。如果在步S10的确定结果是〖no〗，则不计算空隙比。
在计算射入的气泡的运动到一个整倍数n=1之后(使得n成为1的单位时间Δt的倍数)，则准备使时标增加另一个整数时间。整数时间n被从“0”增加到“1”(步S12)，使得对第二个气泡重复步S2。
在步S2，一个新的气泡(第二个气泡)被射入水中，现在共有两个气泡(第一和第二气泡)位于紊流边界层内。在后续过程中，使用随机数发生器进行定向气泡速度的扰动计算。对于第一个气泡，由上述的步骤确定的旧的位置在相继单位时间Δt被改变为新的位置。对于第二个气泡，根据上述的第一气泡的初始位置获得其在单位时间Δt的位置。
通过对预先确定的整数时间数重复这些步骤，获得在每个方格内的气泡的总数，便使得能够确定在某一时刻的空隙比(步S11)。图4和图5表示计算的在沿围长方向的各个方格内空隙比的分布。射入紊流边界层中的气泡沿着和船的运动方向相反的方向运动，沿着底面形成气泡的纵向流。因此，图4和图5表示在船的不同的纵向位置上在x轴上计算的空隙比数据和在y轴上沿围长方向(横向)的不同方格的位置。图4表示在1-5次迭代步骤之后的结果，而图5表示在96-100次迭代之后的结果。
如图4所示，对于1-5次迭代，计算结果十分分散，但是，如图5所示，在足够的迭代次数之后，计算结果则开始收敛。由图5可见，在96-100迭代次数之后，在每个方格内的空隙比似乎收敛于一个常数。
这样，通过进行多次迭代(例如100次)使得足够收敛而确定气泡位置(步S3或S6)，并计算相应的空隙比(步S11)。其中，初始气泡位置(射流出口)被改变(步S13)。这意味着，初始气泡射流出口选择了沿围长方向排列的一个射流出口(方格交点)，但是这被改变为和旧的交点相邻的一个交点。对于第二交点重复步S2-S11。通过对100个气泡对全部相邻的25个射流出口重复这一过程，计算所有的横向射流出口的相应空隙比。
当对所有射流出口进行这种计算时，在步S14计算船的全部水下表面中的空隙比，并检查结果是否收敛。从船射入的气泡沿船的纵向和船相反地运动，因而需要在步S14检查在不同的纵向位置上空隙比的计算结果是否收敛。如果确定结果是〖yes〗，则认为整个仿真过程被完成，如果确定结果是〖no〗，则继续进行仿真过程，直到在全部纵向方格中空隙比收敛。
图6-9表示对于显示S.S值的不同的仿真部分空隙比计算的结果。在图6-9中，S.S值越大，空隙比的分布越接近射流出口。换句话说，对于小的S.S值，空隙比的分布向船的后部移动。在轮廓曲线中时间增量的最小值是1×10-5，而较高的值处于0.002-0.02(单位时间为0.002)的范围内。在任何S.S值下，观察到空隙比的峰值出现在底面上靠近船的后部，还观察到，气泡趋于向着横向的船体表面分散，并且峰值较高。这被认为是沿垂直于船表面方向隔开的方格和流体流动速度的影响。
图10是由计算过程获得的气泡轨迹的底视图。结果表明，大部分气泡在和底面接触的状态下流动。这被认为是沿着底面流动的水的低的相对速度的结果，使得Saffman升力几乎不能影响气泡的行为，并且在紊流边界层薄的底部区域，水的静压梯度低，因而对气泡分散的影响很小，使得升力占主要地位。结果没有表明引起大部分底部气泡从船的中部向横向船体表面移动的这种大的分散。
以下的表1表示在不同的S.S值下在船的后脊骨区域附近没有气泡的流场的计算结果。在表1中，示出的数据用于比较在垂直于船脊骨(水中的)的方向一个方格间距的平均方格间距和同一个方格间距的下游气泡速度。在图10中，可以看到本地的空隙比在沿着底面的一些区域是高的，因为气泡被强迫压向底面。在计算的结果中，空隙比按方格容积平均，从而在垂直方向上方格间距大的底部的中心区域空隙比峰值较小。为了改善这一结果，认为需要匹配方格间距(在水中的)和气泡直径。
表1S.S 垂直于水中的脊骨 方向的一个方格间 距的平均方格间距 垂直于水中脊骨方 向的一个方格间距 的下游流速5 0.0058m 2.258m/s1 0.0042m 2.114m/s1/2 0.0032m 2.057m/s1/4 0.0023m 1.979m/s
图11-14表示半径为1000μm的空隙比分布的纵截面图，用于和表示半径为500μm的气泡的图6-9比较。在船的底面，沿围长方向半径为1000μm的空隙比分布比半径为500μm的气泡呈现较小的分散性，1000μm气泡的峰值比500μm的气泡的峰值较高。这被认为是射流出口数量的影响。气泡黏附于底面的趋势对两种直径似乎没有明显的不同。