鐘強(qiáng)，何勇，姚凱學(xué)，張慶銘，胡加德

(1.貴州大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，貴州貴陽 550025;2.貴州民族大學(xué)數(shù)據(jù)科學(xué)與信息工程學(xué)院，貴州貴陽550025；3.貴州科海新技術(shù)發(fā)展有限公司，貴州貴陽550000)

0 前言

對(duì)閥門的建模往往是建立在閥芯受力分析的基礎(chǔ)上，由于大部分的閥門都采用液壓傳動(dòng)技術(shù)，而活門建模中對(duì)壓力的處理是至關(guān)重要的。以往的研究中，主要采用流量連續(xù)性定理分析壓力[1-6]，將閥口的壓力近似為閥芯處所受的壓力[4]或者根據(jù)泵的工作狀態(tài)將壓力近似為某一函數(shù)[5]，亦或是將它作為未知參數(shù)，設(shè)計(jì)控制策略以對(duì)它進(jìn)行調(diào)節(jié)[6]。某型航空發(fā)動(dòng)機(jī)燃油計(jì)量系統(tǒng)的模型中存在多個(gè)未知參數(shù)，這些參數(shù)的調(diào)節(jié)過程復(fù)雜，并且系統(tǒng)中的流場是湍流場，系統(tǒng)內(nèi)的壓力變化基本沒有規(guī)律，因此將閥口或泵的壓力近似為閥芯處的壓力存在偏差。在實(shí)際中，為獲得閥芯處的具體數(shù)據(jù)，通常非常困難而且花費(fèi)巨大。近年來隨著計(jì)算機(jī)技術(shù)和CFD(Computational Fluid Dynamics)計(jì)算方法的迅速發(fā)展，CFD技術(shù)取得了很大的成就，這使得流場數(shù)值模擬變得高效且可靠[7]。眾多學(xué)者已經(jīng)使用CFD技術(shù)對(duì)多種工藝的各種特性進(jìn)行了研究[8-11]，并將計(jì)算值與實(shí)驗(yàn)所得的結(jié)果進(jìn)行比較，驗(yàn)證了CFD技術(shù)的可靠性。對(duì)于具有隨機(jī)性的時(shí)間序列數(shù)據(jù)，可以采用時(shí)間序列分析方法進(jìn)行分析。時(shí)間序列數(shù)據(jù)本質(zhì)上反映的是某個(gè)或者某些隨機(jī)變量隨時(shí)間不斷變化的趨勢(shì)，而時(shí)間序列預(yù)測(cè)問題的核心就是從數(shù)據(jù)中挖掘出這種規(guī)律，并利用這些規(guī)律作出估計(jì)[12]。近些年來,時(shí)間序列分析已廣泛應(yīng)用于心理學(xué)、地球科學(xué)、數(shù)據(jù)挖掘、數(shù)字化誤差等各方面，多位學(xué)者已經(jīng)對(duì)該方法的有效性進(jìn)行了驗(yàn)證[13]。因此，本文作者采用CFD方法對(duì)活門的流場進(jìn)行數(shù)值模擬，得到其閥芯處的壓力數(shù)據(jù)，結(jié)合時(shí)間序列分析方法對(duì)該數(shù)據(jù)進(jìn)行分析，最終建立計(jì)量活門閥芯的壓力時(shí)間序列預(yù)測(cè)模型。

1 系統(tǒng)工作原理與CFD數(shù)學(xué)模型

1.1 系統(tǒng)工作原理

燃油計(jì)量系統(tǒng)是發(fā)動(dòng)機(jī)的“心臟”，燃油計(jì)量系統(tǒng)的使命是通過控制燃油來控制發(fā)動(dòng)機(jī)。某型燃油計(jì)量系統(tǒng)的工作原理如圖1所示，電子控制器接收位移傳感器反饋信號(hào)，并與輸入信號(hào)作對(duì)比，經(jīng)控制器PI調(diào)節(jié)后輸出計(jì)量活門控制腔壓力，調(diào)節(jié)電液伺服閥的流量，進(jìn)而改變計(jì)量活門控制腔壓力，使計(jì)量活門移動(dòng)，最終改變計(jì)量活門的過流面積，實(shí)現(xiàn)燃油計(jì)量。電液伺服閥的流量同時(shí)也進(jìn)入了指令活門控制腔，可通過改變指令活門的過流面積實(shí)現(xiàn)指令油的釋放[14]。

圖1 燃油計(jì)量系統(tǒng)工作原理

1.2 基本控制方程

某型航空發(fā)動(dòng)機(jī)燃油計(jì)量系統(tǒng)的流場是一個(gè)復(fù)雜的黏性不可壓縮湍流流場，在定常條件下，可以得到流場的連續(xù)性方程：

(1)

1.3 湍流模型

該燃油計(jì)量系統(tǒng)的內(nèi)部存在分離流現(xiàn)象，為更加精確地描述流場，本文作者使用大渦模擬模型進(jìn)行瞬態(tài)模擬計(jì)算。大渦模擬模型對(duì)流場的初值要求較高，因此先采用κ-ε模型對(duì)流場進(jìn)行穩(wěn)態(tài)模擬計(jì)算，再將得到的結(jié)果作為初值進(jìn)行瞬態(tài)模擬計(jì)算[15]。大渦模擬的思路可以簡述為：對(duì)大尺度紊流運(yùn)動(dòng)直接求解N-S方程；利用次網(wǎng)格尺度模型模擬小尺度紊流運(yùn)動(dòng)對(duì)大尺度紊流運(yùn)動(dòng)的影響。

對(duì)于黏性流體，在直角坐標(biāo)系下時(shí)，其運(yùn)動(dòng)規(guī)律受N-S方程控制，簡化不可壓縮流體的N-S方程可以得到式(2):

(2)

從而得到：

(3)

式中：ui、uj分別為i、j方向的速度分量；p為流體壓力；xi、xj分別為i、j方向的位移；t為時(shí)間。

2 流場數(shù)值模擬計(jì)算

2.1 網(wǎng)格劃分及邊界條件設(shè)置

由于燃油計(jì)量系統(tǒng)結(jié)構(gòu)較為復(fù)雜，選擇在AutoCAD軟件中進(jìn)行系統(tǒng)三維建模，隨后將該三維模型導(dǎo)入CFD軟件中進(jìn)行網(wǎng)格劃分。圖2所示為燃油計(jì)量系統(tǒng)的三維模型。

圖2 燃油計(jì)量系統(tǒng)的三維模型 圖3 系統(tǒng)網(wǎng)格劃分結(jié)果

煤油從進(jìn)油口進(jìn)入系統(tǒng)，然后通過內(nèi)部進(jìn)油口進(jìn)入計(jì)量活門控制腔,作用到計(jì)量活門閥芯和指令活門閥芯上，通過出油口流出。因?yàn)樵撓到y(tǒng)結(jié)構(gòu)復(fù)雜，所以采用非結(jié)構(gòu)網(wǎng)格進(jìn)行劃分，共生成約156萬個(gè)網(wǎng)格。劃分好的網(wǎng)格模型如圖3所示。

網(wǎng)格劃分完畢后，需進(jìn)行邊界條件設(shè)置。該研究中的流體介質(zhì)為煤油，為不可壓縮流體，密度為780 kg/m3、動(dòng)力黏度為2.4×10-3Pa·s；進(jìn)口為壓力進(jìn)口，進(jìn)口壓力恒為7 MPa，進(jìn)口處的流速為10 m/s；出口為壓力出口，出口壓力恒為6 MPa；壁面均為無滑動(dòng)壁面。

2.2 仿真結(jié)果

將劃分好的網(wǎng)格模型導(dǎo)入Fluent軟件，設(shè)置好邊界條件，對(duì)流場進(jìn)行穩(wěn)態(tài)模擬計(jì)算。選用κ-ε湍流模型求出流場的初值，在該初值的基礎(chǔ)上選用LES湍流模型進(jìn)行瞬態(tài)求解，并在計(jì)量活門的閥芯處設(shè)置監(jiān)測(cè)器，每一個(gè)time step便計(jì)算一次閥芯處的壓力，共得到4 600個(gè)數(shù)據(jù)。最終得到的系統(tǒng)流場壓力分布云圖和計(jì)量活門閥芯表面壓力分布云圖分別如圖4(a)(b)所示。

圖4 仿真結(jié)果

從圖4(b)可知：在計(jì)量活門閥芯表面存在明顯的漩渦，并且監(jiān)測(cè)器得到的數(shù)據(jù)表明閥芯處的壓力在(6.7±0.2)MPa范圍內(nèi)隨機(jī)變動(dòng)。

3 時(shí)間序列原理與壓力模型的建立

3.1 時(shí)間序列建模原理

設(shè){xt,t∈T}是一個(gè)時(shí)間序列，稱滿足如式(4)所示的結(jié)構(gòu)模型為自回歸移動(dòng)平均模型，簡記為ARMA(p,q)，

(4)

式中:εt為均值為零的白噪聲序列；φ0為常數(shù)；p、q為ARMA模型的階數(shù)；φi為待估計(jì)的自回歸參數(shù)；θj為待估計(jì)的滑動(dòng)平均參數(shù)。

自回歸滑動(dòng)平均(Autoregressive Moving Average，ARMA)模型是研究時(shí)間序列的主要方法，由自回歸(Autoregressive，AR)模型和滑動(dòng)平均(Moving Average，MA)模型混合而成。對(duì)于具有時(shí)間趨勢(shì)和周期特征的數(shù)據(jù)可以進(jìn)行差分，使它轉(zhuǎn)化為平穩(wěn)序列。

在自回歸移動(dòng)平均模型中，假設(shè)殘差序列是滿足方差齊性條件的，即殘差的方差始終為一個(gè)常數(shù)。若某一模型的殘差服從均值為零，方差是一個(gè)隨時(shí)間變化的量且這個(gè)隨時(shí)間變化的方差是該時(shí)間點(diǎn)過去有限噪聲值平方的線性組合的正態(tài)分布時(shí)，該模型稱為自回歸條件異方差模型(ARCH)。ARCH(q)模型的結(jié)構(gòu)如下：

(5)

式中:系數(shù)β0>0,βk≥0,k=1,2,…,q;Xt為前定解釋變量，包括被解釋變量的滯后項(xiàng);η為回歸參數(shù);εt為獨(dú)立同分布，εt～N(0,1)；ut為殘差項(xiàng)?？梢姉l件方差ht隨{ut}過去值的變化而變化。

3.2 時(shí)間序列建模步驟

對(duì)時(shí)間序列建模一般遵循4個(gè)步驟：模型識(shí)別、模型估計(jì)、模型檢驗(yàn)和模型應(yīng)用。

(1)模型識(shí)別。根據(jù)數(shù)據(jù)的特征，判斷所研究的時(shí)間序列屬于哪一類，主要檢驗(yàn)?zāi)Ｐ偷钠椒€(wěn)性，只有當(dāng)序列通過平穩(wěn)性檢驗(yàn)時(shí)才能使用時(shí)間序列分析模型。檢驗(yàn)的方法主要有DF檢驗(yàn)、ADF檢驗(yàn)、PP檢驗(yàn)、協(xié)整檢驗(yàn)等[17]。

(2)模型估計(jì)。依照樣本信息進(jìn)行模型識(shí)別后，得到所分析的時(shí)間序列大概的模型類型和模型結(jié)構(gòu)，模型的最終形式還需要估計(jì)模型的參數(shù)后才能確定。常用的參數(shù)估計(jì)方法：矩估計(jì)、極大似然估計(jì)和最小二乘估計(jì)等[18]方法。

(3)模型檢驗(yàn)。在模型識(shí)別時(shí)，為簡化問題提出了一些假設(shè)，因此必須對(duì)模型進(jìn)行檢驗(yàn)。時(shí)間序列模型的檢驗(yàn)有兩類：一類是模型的顯著性檢驗(yàn)；另一類是模型參數(shù)的顯著性檢驗(yàn)，這兩類檢驗(yàn)統(tǒng)稱為模型的診斷性檢驗(yàn)。此外，對(duì)某些序列，還需要進(jìn)行ARCH效應(yīng)檢驗(yàn)。

(4)模型應(yīng)用。時(shí)間序列模型的應(yīng)用主要包括變量動(dòng)態(tài)結(jié)構(gòu)分析、預(yù)測(cè)和控制。其中預(yù)測(cè)是時(shí)間序列建模的最主要目的，是用已經(jīng)估計(jì)出參數(shù)的模型預(yù)報(bào)變量未來的變化。

3.3 壓力時(shí)間序列模型的建立

通過流場仿真得到了計(jì)量活門閥芯的壓力數(shù)據(jù)共4 600個(gè)，該數(shù)據(jù)是在相同時(shí)間間隔內(nèi)經(jīng)軟件計(jì)算得到的具有隨機(jī)性的一組數(shù)據(jù)，滿足時(shí)間序列的特征，可以用時(shí)間序列分析方法對(duì)該數(shù)據(jù)進(jìn)行分析，取前4 500個(gè)數(shù)據(jù)進(jìn)行時(shí)間序列的建模，剩下的部分?jǐn)?shù)據(jù)用于模型的預(yù)測(cè)對(duì)比。本文作者使用R軟件進(jìn)行時(shí)間序列的建模。

(1)對(duì)原始數(shù)據(jù)進(jìn)行平穩(wěn)性檢驗(yàn)，常用的方法為ADF(Augmented Dickey-Fuller)檢驗(yàn)。ADF檢驗(yàn)可以用于3種自回歸模型的單位根檢驗(yàn)，即不含漂移項(xiàng)和時(shí)間趨勢(shì)項(xiàng)的自回歸模型、只含漂移項(xiàng)的自回歸模型以及兩者都含的自回歸模型。表達(dá)式分別為

xt=φ1xt-1+…+φpxt-p+εt

(6)

xt=μ+φ1xt-1+…+φpxt-p+εt

(7)

xt=μ+βt+φ1xt-1+…+φpxt-p+εt

(8)

以式(6)為例，將式(6)變形為

xt-xt-1=φ1xt-1+…+φpxt-p-xt-1+εt=

(φ2+…+φp)xt-1+φ1xt-1-xt-1-(φ2+…+φp)xt-1+

φ2xt-2+(φ3+…+φp)xt-2-(φ3+…+φp)xt-2+

φ3xt-3+(φ4+…+φp)xt-3+…-φpxt-p+1+φpxt-p+εt整理得：

(9)

其中：ρ=φ1+φ2+…+φp-1，βj=-φj+1-φj+2-…-φp,j=1,2,…,p-1。

相似地，也可以得到式(7)(8)的單位根檢驗(yàn)式為

(10)

(11)

單位根檢驗(yàn)的假設(shè)條件:

原假設(shè)H0：ρ=0(序列非平穩(wěn))?H1:ρ<0(序列平穩(wěn))

由此可以構(gòu)造ADF統(tǒng)計(jì)量：

通過蒙特卡洛方法，可以得到τ檢驗(yàn)統(tǒng)計(jì)量，并進(jìn)行檢驗(yàn)。共有3個(gè)ADF檢驗(yàn)式，但是針對(duì)實(shí)際應(yīng)用中應(yīng)采用哪一個(gè)檢驗(yàn)式，DOLADO等建議從式(11)開始檢驗(yàn)，若得出的結(jié)論為接受單位根假設(shè)時(shí)，還需謹(jǐn)慎考慮所用檢驗(yàn)式是否包含多余的時(shí)間趨勢(shì)項(xiàng)，即也應(yīng)對(duì)φ3統(tǒng)計(jì)量進(jìn)行檢驗(yàn)[19]。若φ3統(tǒng)計(jì)量的檢驗(yàn)結(jié)果為接受原假設(shè)，則認(rèn)為用式(11)進(jìn)行單位根檢驗(yàn)不合適，改用式(10)重新進(jìn)行單位根檢驗(yàn)。同理，若式(10)檢驗(yàn)結(jié)果接受原假設(shè)，仍需對(duì)φ1進(jìn)行檢驗(yàn)，根據(jù)所得結(jié)果進(jìn)一步考慮是否使用式(9)重新進(jìn)行單位根檢驗(yàn)。當(dāng)序列通過平穩(wěn)性檢驗(yàn)后，便可對(duì)該序列進(jìn)行定階，以求得其自相關(guān)系數(shù)圖(ACF)和偏自相關(guān)系數(shù)圖(PACF)進(jìn)而推斷模型的階數(shù)，還可以借助AIC準(zhǔn)則和BIC準(zhǔn)則來進(jìn)行定階。用含漂移項(xiàng)和時(shí)間趨勢(shì)項(xiàng)的檢驗(yàn)式(11)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)，檢驗(yàn)結(jié)果如表1所示。

表1 用式(11)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)的結(jié)果

從表1可以看出:φ3的檢驗(yàn)值只比顯著性水平為10%的臨界值略高一點(diǎn)，這表明用含漂移項(xiàng)和時(shí)間趨勢(shì)項(xiàng)的檢驗(yàn)式(11)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)是不合適的。因此，進(jìn)一步選擇只含漂移項(xiàng)的檢驗(yàn)式(10)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)，檢驗(yàn)結(jié)果如表2所示。

表2 用式(10)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)的結(jié)果

從表2中可以看出:φ1的檢驗(yàn)值比顯著性水平為1%的臨界值低，這表明用含漂移項(xiàng)的檢驗(yàn)式(10)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)是不合適的。因此，進(jìn)一步選擇不含漂移項(xiàng)和趨勢(shì)項(xiàng)的式(9)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)，檢驗(yàn)結(jié)果如表3所示。

表3 用式(9)對(duì)數(shù)據(jù)進(jìn)行ADF單位根檢驗(yàn)的結(jié)果

從表3可以看出:τ1的檢驗(yàn)值遠(yuǎn)高于顯著性水平為1%的臨界值，則接受原假設(shè)，即認(rèn)為該序列是非平穩(wěn)序列。

此外，本文作者還對(duì)原始數(shù)據(jù)的純隨機(jī)性進(jìn)行了檢驗(yàn)，結(jié)果表明數(shù)據(jù)并非白噪聲序列，即原始序列具有密切的相關(guān)關(guān)系，需要對(duì)它進(jìn)一步分析。

(2)由于原始序列非平穩(wěn)，對(duì)序列進(jìn)行一階差分，對(duì)差分后的序列再進(jìn)行單位根檢驗(yàn)，檢驗(yàn)結(jié)果表明該序列為平穩(wěn)序列。求出該平穩(wěn)序列的自相關(guān)系數(shù)(ACF)和偏自相關(guān)系數(shù)(PACF)分別如圖5、圖6所示。

圖5 一階差分后序列的自相關(guān)系數(shù)(ACF) 圖6 一階差分后序列的偏自相關(guān)系數(shù)(PACF)

從圖5可以看出:序列的自相關(guān)系數(shù)在3階截尾。

從圖6可以看出:序列的偏自相關(guān)系數(shù)在2階截尾。因此，擬合的ARMA(p,q)模型中,參數(shù)p=2、q=3。使用極大似然估計(jì)方法對(duì)參數(shù)進(jìn)行估計(jì)，所得的自回歸移動(dòng)平均模型為

xt=0.109 4xt-1+0.732 3xt-2+ut-0.560 4ut-1-

0.817ut-2+0.525 1ut-3

(12)

式中：ut為殘差序列。

(3)對(duì)該模型的殘差進(jìn)行白噪聲檢驗(yàn)，結(jié)果表明該序列是白噪聲序列。對(duì)該序列進(jìn)行ARCH效應(yīng)檢驗(yàn)，常用的ARCH效應(yīng)檢驗(yàn)方法是LM檢驗(yàn)和Portmanteau Q檢驗(yàn)。Portmanteau Q檢驗(yàn)的思想：如果殘差序列方差非齊且具有群集效應(yīng)，那么殘差平方序列通常具有自相關(guān)性，故可將方差非齊次的檢驗(yàn)轉(zhuǎn)化為殘差平方序列的自相關(guān)檢驗(yàn)。該檢驗(yàn)的假設(shè)條件[20]為

H0:ρ1=ρ2=…=ρq=0?H1:ρ1,ρ2,…,ρq不全為0。

其中：ρk表示殘差平方延遲k階的自相關(guān)系數(shù)。

Portmanteau Q檢驗(yàn)統(tǒng)計(jì)量的構(gòu)造為

其中：n為序列長度;ρi為殘差序列延遲i階的自相關(guān)系數(shù)，表達(dá)式為

其中：

若原假設(shè)成立，則統(tǒng)計(jì)量近似服從χ2分布，其自由度為q-1。因此，當(dāng)檢驗(yàn)統(tǒng)計(jì)量Q(q)的P值小于顯著性水平時(shí)，拒絕原假設(shè)，殘差平方序列具有自相關(guān)關(guān)系。使用Portmanteau Q方法對(duì)該序列進(jìn)行ARCH效應(yīng)檢驗(yàn)，檢驗(yàn)結(jié)果如表4所示。

表4 ARMA(2，3)殘差平方的Q檢驗(yàn)結(jié)果

從表4可知:ARMA(2,3)的Q檢驗(yàn)P值均遠(yuǎn)小于顯著性水平0.05，因此拒絕原假設(shè)，殘差平方序列具有自相關(guān)關(guān)系，序列的ARCH效應(yīng)顯著。對(duì)差分后的平穩(wěn)序列擬合ARMA(2,3)-ARCH(1)模型，最終模型的各參數(shù)如表5所示。

由表5可知:所有參數(shù)的P值都小于0.05，即ARMA-ARCH模型的參數(shù)擬合是顯著的。綜上，得到完整的擬合模型為

表5 ARMA-ARCH模型參數(shù)估計(jì)及其檢驗(yàn)結(jié)果

(13)

式中:ut為序列的殘差項(xiàng);εt為獨(dú)立同分布，且εt～N(0,1);ht為殘差的條件方差。

(4)使用ARMA(2，3)-ARCH(1)對(duì)序列進(jìn)行預(yù)測(cè)，并將原序列畫出來作為對(duì)比，結(jié)果如圖7所示，其中誤差如表6所示。

表6 ARMA-ARCH模型的預(yù)測(cè)誤差

圖7 ARMA-ARCH模型的預(yù)測(cè)結(jié)果

從圖6可以看出：該模型擬合較好，而且預(yù)測(cè)到數(shù)據(jù)在短期內(nèi)呈上升趨勢(shì)。

從表6中可以看出:模型的均方根誤差小于1，模型擬合較好。

4 結(jié)論

(1)在Fluent軟件中使用大渦模擬湍流模型對(duì)燃油計(jì)量系統(tǒng)進(jìn)行了流場數(shù)值模擬計(jì)算，并通過在計(jì)量活門閥芯處設(shè)置monitor的方法獲得了計(jì)量活門閥芯處的壓力時(shí)間序列。

(2)使用時(shí)間序列分析方法對(duì)該序列進(jìn)行分析，建立了該序列的自回歸移動(dòng)平均模型；對(duì)該自回歸移動(dòng)平均模型殘差的平方進(jìn)行了異方差檢驗(yàn)(ARCH效應(yīng)檢驗(yàn))，建立了該模型殘差的異方差模型。結(jié)合該序列的自回歸模型和異方差模型，建立了原始序列的ARMA(2,3)-ARCH(1)模型，使用該模型預(yù)測(cè)了5個(gè)值并與原數(shù)據(jù)進(jìn)行比較，結(jié)果表明預(yù)測(cè)效果較好。

(3)研究結(jié)果為相關(guān)研究中對(duì)壓力變量的處理提供了參考。