一种水下探测装置无动力运动预报自动计算方法与流程

本发明属于水下运动预报领域，更具体地，涉及一种水下探测装置无动力运动预报自动计算方法。

背景技术：

潜艇的下潜和上浮利用了物理学中的浮力原理，进行无动力下潜上浮。在潜艇的内外壳体之间，通常设有10来个水柜，用它来控制潜艇的下潜和上浮。下潜时，可以往水柜中灌水，使潜艇沉人水中；上浮时，用高压气体将水柜中的水压出水柜，潜艇便可浮出水面。为了调整潜艇平衡，还设有若干个小水柜，通过调整水量的大小来控制潜艇的稳定和平衡。为了保持和改变航向，潜艇通常在潜艇的尾部装有方向舵(垂直舵)；为了保持和改变深度，通常分别在潜艇的前部和尾部还装有升降舵(水平舵)。

基于气球原理设计的第一代和第二代载人潜水器均是采用无动力下潜和上浮的方式，这样做最大好处就是可以大大地节省潜水器上的能源，但需要解决无动力潜浮过程中水动力学问题。

目前在现有技术中，已经有技术在对相关的水动力学问题进行研究，例如在蛟龙号设计阶段刘正元等就开始研究载人潜水器的无动力下潜上浮运动给出了一套简化分析方法并应用于蛟龙号载人潜水器的无动力潜浮设计，但是上述技术存在如下问题：利用牛顿力学定理和微分方程进行数学分析，并未考虑装置的动力学模型在水下所受水动力的影响。事实上，流体带来的水动力作用对水下物体的运动有着非常重要的影响，复杂的水动力学系数是分析流场中物体运动规律的重要参数。随着Fluent软件功能的日渐完善，张军等人采用非结构动网格方法已完成对含多介质流场的数值模拟，李迎华、吴宝等人,应用动网格技术对一典型水下低速航行体缩比模型主体及带附体的水动力进行了计算,结果与风洞模型试验结果吻合较好。

但是，在上述现有技术中，存在如下有待解决的问题：(1)动网格生成时边界网格点容易发生交错和网格点丢失，而且根据新的物面生成网格计算量大，需要耗费大量的时间；(2)装置在水下的连续运动，流体阻力与装置运动速度相关，采用离散化的方法单步计算每一分段的阻力和速度，十分耗费时间。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种水下探测装置无动力运动预报自动计算方法，该方法兼顾装置的动力学模型及海流水动力影响，使用简单、可靠、自适应性好的非结构化网格生成技术，能够准确、高效地实现无动力深海装置的水下运动预报。

为实现上述目的，按照本发明提供了一种水下探测装置无动力运动预报自动计算方法，其特征在于，所述自动计算方法包括如下步骤：

(1)自定义流场速度传递给Fluent

通过用户自定义函数完成速度自动输入；

(2)Fluent阻力自动传出

设置所述水下探测装置模型和网格模型，利用等时长叠加算法通过自

动计算方法实现Fluent软件连续计算每一分段的阻力、速度、位移由

此完成所述水下探测装置无动力运动预报。

进一步地，所述等时长叠加算法为：

取所述水下探测装置无动力下潜上浮垂直方向上的等时间步长为Δt，当Δt取得极其小，变加速运动的各个分段可视为匀加速运动，运用改进的Euler公式可计算出第K个分段的初末速度V_k-1和V_k、加速度a_k、位移S_k；

第K个分段的初速度V_k-1已知，由于时间步长极小，则第K分段的预报加速度a_k⁰可表示为f(V_k-1)，运用显示欧拉法，可得第K分段结束时的预报速度V_k⁰＝V_k-2+2Δtf(V_k-1)，再运用欧拉矫正公式可得：

第K分段末速度

$V_k=V_{k-1}+\frac{\mathrm{\Δ}t}{2}\[f\left(V_{k-1}\right)+f\left({V_k}^0\right)\]$

第K分段加速度

a_k＝f(V_k)

第K分段位移

$S_k=\frac{V_k+V_{k-1}}{2}\mathrm{\Δ}t$

每段依次叠加计算，可得出每一分段的速度、加速度、位移，从而得到装置上浮下潜垂直方向上的的运动模拟过程。

进一步地，所述自动计算方法按照如下方法实现：

通过在Scheme文件中写入文件循环语句，使得Journal文件不断被调用，所述Journal文件中包含所述等时长叠加算法的操作命令。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：充分考虑水动力系数和装置动力学模型的影响进行仿真，仿真结果更有可靠性；采用非结构体网格划分技术，自适应水下探测装置复杂外轮廓的性能好；采用等时长叠加算法，能够将复杂连续的水下运动离散化，简化了计算；采用自动计算的方法，不仅能够大大缩减工作时间，还能很大程度上地减少由于人为导致的错误。

附图说明

图1为按照本发明实现的基于等时长叠加算法的水下运动模拟图；

图2为按照本发明实现的流体一体化模型下潜运动计算过程图；

图3为按照本发明实现的流体一体化模型自动计算流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明提供了一种装置的动力学模型仿真的方法，在于模型的绘制、模型的网格划分、流场设置、模型阻力求解；所述模型的绘制，即流场和装置水下一体化模型，包括流场和装置的绘制、流场装置一体化对称模型绘制；所述网格划分，包括模型导入、网格类型选择、参数设置；所述流场设置，包括网格导入、流体类型设置、流场设置；所述模型阻力求解，包括求解方程类型设置、初始参数输入、迭代求解次数设置、阻力值输出。

网格划分，采用简单、可靠、自适应性能好的非结构网格生成技术，非结构网格对模型的自适应性好，工程人员工作量小，但是计算量大，对计算机要求高，适用于结构复杂的模型，深海探测装置外轮廓较为复杂，故选择非结构网格生成技术来划分网格，其次需要确定划分网格的类型，分析装置在流场中的运动，需对装置外轮廓及周围流体子域划分网格，生成网格必须是体网格。

按照本发明的另一方面，提供了一种自动计算的方法，其特征在于等时长叠加算法、自动计算方法；所述等时长叠加算法，是基于预报-矫正系统改进的Euler公式和梯形公式所建立的一种累加算法，其通过自动计算的方式实现；所述自动计算方法，包括自定义流场速度传递、阻力自动传出、计算过程自动进行三个部分构成。

优选地，自定义流场速度传递、阻力自动传出、计算过程自动进行三个过程，采用用户自定义函数UDFs和Journal文件批处理实现。

本发明提出一种基于等时长叠加算法的装置水下运动预报的自动计算方法，目的在于模拟预报各种水下装置的无动力下潜上浮运动。

首先，按照本发明实现的等时长叠加算法如下：如图1所示为水下装置无动力下潜上浮垂直方向上的等时间步长运动模拟图，时间步长取为Δt，当Δt取得极其小，变加速运动的各个分段可视为匀加速运动，运用改进的Euler公式可计算出个分段的初末速度(V_k-1和V_k)、加速度(a_k)、位移(S_k)。

以图中分段的第K分段为例来进行说明每个分段的初末速度、加速度和位移值计算如下：

初速度V_k-1(也是前一分段的末速度)已知，由于时间步长极小，则第K分段的预报加速度a_k⁰可表示为f(V_k-1)，运用显示欧拉法，可得第K分段结束时的预报速度V_k⁰＝V_k-2+2Δtf(V_k-1)(公式中的Δt显示不出，请修正)，再运用欧拉矫正公式可得：

第K分段末速度

$V_k=V_{k-1}+\frac{\mathrm{\Δ}t}{2}\[f\left(V_{k-1}\right)+f\left({V_k}^0\right)\]---\left(1\right)$

加速度

a_k＝f(V_k) (2)

位移

$S_k=\frac{V_k+V_{k-1}}{2}\mathrm{\Δ}t---\left(3\right)$

依次叠加计算，可得出每一分段的速度、加速度、位移，从而得到装置上浮下潜垂直方向上的的运动模拟过程。这就是等时长叠加算法。同理可以根据这个方法，模拟预报装置在水平方向上在水流的作用下的运动过程。f(v_k-1)表示装置在自身重力浮力和流体阻力的作用下的加速度，装置的重力和浮力已知，因为流体阻力和装置的速度有关，所以我们利用FLUENT软件，将装置网格模型输入到FLUENT软件中，设置相关参数设置流体子域并初始化流场，输入速度值，进行迭代计算，得出装置的相应的流体阻力。

所以为了得到装置的流体阻力，我们需将装置与流体域的一体化模型用Solidworks软件画出来，导入ANSYS中进行网格划分，划分完毕之后将网格导入Fluent软件之中计算。得出分段的流体阻力之后，便可计算出其加速度，利用本发明提及的改进的Euler法，可获得下一分段的速度，再利用Fluent软件，可计算出这一分段的流体阻力，依次类推，得出所有分段的初末速度、加速度、位移，包括速度达到稳定时候的值，也就模拟出了装置在水下的运动。

但是我们可以看到，在上述的叠加算法中，为了得到最后稳定阻力值，需多次手动输入速度，初始化流场、迭代计算，并且每一步都要手动输出阻力到表中且结算下一分段初速度，非常消耗时间，并且这段时间工作人员不能离开。

在本发明中，更进一步地，为针对上述算法的执行过程，对fluent进行二次开发，使其能够自动计算。

在软件的执行过程中，自动计算过程包括主要包括如下的三个主要的步骤:

自定义流场速度传递给Fluent；

FLUENT阻力自动传出；

自动执行仿真过程。

(1)自定义流畅速度传递给Fluent

自定义流场速度自动传入FLUENT的关键是用户自定义函数(User－Defined Functions，即UDFs)，UDFs可以提高FLUENT程序的标准计算功能。它可以自动读取文本中速度值，并输入到FLUENT速度端口。

它是用C语言书写的，有两种执行方式：interpreted型和compiled型。Interpreted型比较容易使用，但是可使用代码(C语言的函数等)和运行速度有限制。Compiled型运行速度快，而且也没有代码使用范围的限制，但使用略为繁琐。将C语言程序在记事本中写好之后，保存为.C文件，在FLUENT软件UDF功能中进行编译、下载。在确定编译没有错误之后，就可以自定义流场速度的自动输入了。UDF函数加载完毕之后，点击边界条件，对流场入口“in”设置速度，选择“udf inlet_x_velocity：：labudf”。

(2)FLUENT阻力自动传出

速度自动输入之后，就需要实现仿真过程的自动执行，以及当前分段的阻力自动传出。

自动计算传出参数为装置外轮廓上的流体阻力，该流体阻力包括压差阻力和粘性阻力，要想将该参数传递到FLUENT外部做进一步处理，只需使用FLUENT软件中的Reports功能下的write功能即可。write功能可以将仿真结果也就是阻力值传出到外部文本中。由于每一步仿真结果都要执行write，所以将其写到仿真过程自动执行Journal文件当中，实现每一分段的阻力值自动传出。

而Journal文件的作用就像windows下的批处理文件一样，文件中包含了一连串的操作命令。Fluent读入journal文件后，它就会按照顺序执行一系列的命令。因为Fluent的journal相对比较难编写，所以一个简单的方法就是，设置的过程中把设置过程保存成journal，相当于先在线录制设置操作，下次直接调用而不用自行编写journal文件。具体方法：

(1)FLUENT—File—write—start journal—输入文件名，

(2)读入case—设置相关参数—stop journal。这样journal就生成了，下次使用，直接读入就行了，Fluent—File—read—journal。

本实施例中录制的journal文件包含以下内容：流场入口速度初始化、仿真步数设置、Profile Update Interval设置、输出阻力值。

Journal文件与scheme文件的配合使用促成了仿真的自动化，二者区别在于scheme文件只能通过编写scheme程序得到，而journal文件不仅可以通过编程得到，还可以直接通过录制操作得到，更加方便。引入scheme文件的原因在于journal文件自身不能循环也不能够嵌套调用，所以只能通过在scheme文件中写入使得journal文件循环执行的语句来实现不间断的计算。自动计算过程之所以能够实现自动计算，主要是由于引进了scheme文件，其功能是实现journal文件的循环，以致于整个仿真过程能够连续执行下去，无需人来干预。

如图1所示为基于等时长叠加算法的装置水下运动模拟图，具体操作如下，即如上所述的当速度自动输入之后，仿真过程的自动执行按照下述的步骤开始执行：

(1)已知装置在初始零时刻的运动速度V₀＝0，可以计算装置此时刻加速度a₀＝f(V₀)；

(2)用欧拉法计算装置第一分段结束时装置速度V₁＝V₀+f(V₀)Δt，此刻加速度a₁＝f(V₁)，第一分段位移

(3)用欧拉公式计算第二分段结束时预报速度V₂⁰＝V₀+2Δtf(V₁)，获取第二分段的预报加速度用梯形公式校正该速度得到同样可求得此刻加速度a₂＝f(V₂)，第二分段位移为

(4)按照第(3)步方法可求出第3分段结束时预报速度V₃⁰＝V₁+2Δtf(V₂)，校正后速度加速度a₃＝f(V₃)，该分段位移

(5)以此类推，可求得第k段结束时预报速度V_k⁰＝V_k-2+2Δtf(V_k-1)，校正速度、加速度、第k分段位移分别为

(式中为第k段结束由欧拉法得到预报速度，f(V_k⁰)为预报速度对应加速度)；

(6)如此计算下去，设第n分段为最后一段变加速运动，第n分段以后装置达到动态平衡做匀速运动，以速度V_n完成剩下m-n段分运动，则装置运动总位移可表示如下：

根据以上公式，可求出装置下潜时间、上浮时间、回收预计时间以及在海平面装置回收范围。

整个上述的稳态叠加算法的在程序上大的实现主要包括Fluent单步计算和数据处理两部分内容。

Fluent单步计算：Fluent单步计算由前处理、求解和后处理三个阶段组成：前处理包括模型建立、网格划分、边界条件设立等步骤；求解基本思路是将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积，再对待解微分方程的每一个控制体积积分以完成对方程的离散化进而求解；后处理能够对仿真结果进行收敛判断，还可以输出流体阻力。

数据处理:仅通过Fluent软件不能模拟装置的上浮过程，故处理Fluent单步计算结果不可或缺.数据处理的目的是将前一分段的流体阻力通过欧拉预报梯-形校正数值积分法得出后一分段的流体速度，然后以该速度为下一分段初速度求得下一分段的流体阻力，使整个运动过程连续。

在一个具体的实施例中，如图2、3所示，按照本发明实现的水下探测装置无动力运动预报及其自动计算方法主要执行如下步骤来进行相应的计算：

(a)建立模型与流体一体化模型

本发明主要利用Solidworks建立模型(其他绘图软件均可)，在模型建立过程中，我们将装置模型在流场中的运动状态与流场画成一体化模型，取其对称面的一侧，这样有利于在ANSYS软件中进行流场出入口确定。流场区域要足够大，以防流场边缘湍流影响装置阻力的计算。

(b)将模型导入ANSYS中进行网格划分

将几何模型文件另存为后缀名为IGS的文件，导入ANSYS ICEM软件中划分网格。使用ICEM划分网格前先对几何模型进行检查和修正，并将模型划分为INLET(流体入口)、OUTLET(流体出口)、WALL(流体边界)、SYMMETRY(流体子域对称面)、BODY(装置外轮廓)五个组成部分。

(c)网格划分完之后，将网格文件导入到fluent中，设置相关参数

设置流体类型和流体域，接着选择流体子域，双击“water”，选择water-liquid；边界条件设置，边界in选为速度入口velocity-inlet,并点击“Edit”后，在“Velocity Magnitude”栏输入速度大小，边界out选为outflow，其他边界不改变；solution initialization，computer from“in”,再点击“Initialize”；然后，点击“Run Calculation”，设置迭代次数，点击计算。在流体域、流体子域、边界条件设置完成以后，将文件保存为CAS类型，下次打开时可直接输入速度进行求解流体阻力。

(d)手动输入速度，初始化流场，计算阻力

将速度值输入到velocity-inlet中，初始化流场入口，设置迭代参数，计算阻力值，将的到的阻力值导入表格。

(e)判定阻力是否达到稳定，若稳定，则结束计算；若不稳定，则跳回

步骤(d)继续计算。

由图3所示，自动计算的主要步骤如下：

(a)编写自定义函数(User－Defined Functions，即UDFs)，其功能是实现从文档中读取速度自动输入至Fluent中；

(b)编写Journal文件可scheme文件。自动计算过程之所以能够实现自动计算，主要是由于引进了scheme文件，其功能是实现journal文件的循环，以致于整个仿真过程能够连续执行下去，无需人来干预。而Journal程序的功能则是实现流场入口速度初始化、仿真步数设置、Profile Update Interval设置、输出阻力值。Journal文件与scheme文件的配合使用促成了仿真的自动化，二者区别在于scheme文件只能通过编写scheme程序得到，而journal文件不仅可以通过编程得到，还可以直接通过录制操作得到，更加方便。引入scheme文件的原因在于journal文件自身不能循环也不能够嵌套调用，所以只能通过在scheme文件中写入使得journal文件循环执行的语句来实现不间断的仿真。

(c)在程序中我们需要使用四个文档，用于存放变量，分别如下：

用于存放原始阻力数据的文件：a.txt

用于存放当前速度的文件：temp.txt

用于存放循环次数m的文件：b.txt

用于存放过程中产生的所有阻力和速度的记录文件：c.txt

所以，我们在仿真进行前需要在指定路径建立这四个文档。

(d)打开Fluent软件，打开手动计算保存的.CAS文件，导入编译UDFs函数，第一步手动执行仿真计算，将初始速度手动输入到Fluent输入端口，手动启动仿真过程，将第一步的阻力手动输入到a.txt；调用scheme文件，执行自动动计算。

将得到的结果导入excel表格，生成折线图与手动计算结果对比。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。