基于質(zhì)量集中法和迭代法的平面網(wǎng)衣變形模擬

楊靜宇，李 靖

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點實驗室，上海200240）

0 引言

隨著經(jīng)濟的發(fā)展，人們的生活水平大幅提升，對水產(chǎn)品的需求與日俱增[1]，漁業(yè)的健康發(fā)展已成為全世界眾多國家所共同關(guān)心的問題。我國從1970年代開始發(fā)展淺水網(wǎng)箱，到2018年全國傳統(tǒng)近岸網(wǎng)箱數(shù)量已達到149.3萬只，2019年深水網(wǎng)箱的養(yǎng)殖體積已經(jīng)達到19 358 969 m3[2]。相比于試驗，網(wǎng)衣的數(shù)值模擬要更加簡便和經(jīng)濟[3]。近年來，隨著科技的迅速發(fā)展，計算機的計算能力也突飛猛進，有限元方法、集中質(zhì)量法和質(zhì)量加彈簧的模擬方法被廣泛地應(yīng)用到網(wǎng)衣數(shù)學(xué)建模和動態(tài)模擬上。NIEDZWIED和HOPP[4]以網(wǎng)目為基本單位預(yù)測了網(wǎng)衣的形狀和載荷，并計算出中層拖網(wǎng)的三維形狀；PRIOUR等[5]利用有限元的方法計算出網(wǎng)衣下垂的形狀，并采用有限元方法研究剛度對網(wǎng)衣的影響；LEE等[6]將網(wǎng)衣離散成彈簧和質(zhì)量點的結(jié)合體，模擬出拖網(wǎng)在運動時的狀態(tài)；TAKAGI等[7]考慮了漁網(wǎng)的彈性對網(wǎng)衣變形的影響，建立漁網(wǎng)的數(shù)學(xué)模型，對漁網(wǎng)的運動狀態(tài)進行動態(tài)模擬，在此基礎(chǔ)上開發(fā)出軟件NALA來計算網(wǎng)衣的變形；AARASNES等[8]把柔性的網(wǎng)衣簡化為一些平面網(wǎng)片鉸接在一起，將這些小網(wǎng)片視為由純二維平面組成，忽略了這些網(wǎng)片的三維特性。在國內(nèi)，李玉成等[9]建立了各集中質(zhì)量點的運動方程，通過Runge-Kutta-Verner方法求解運動方程，模擬結(jié)果與試驗相比吻合較好；袁軍亭[10]利用有限元理論，建立網(wǎng)衣的三維力學(xué)數(shù)學(xué)模型（3DNNM），計算出平面網(wǎng)衣和圓柱形網(wǎng)片的張力分布和總阻力，與試驗結(jié)果相比誤差較??；趙云鵬等[3，11]將網(wǎng)衣離散成質(zhì)量點和彈簧的組合體，求解得到網(wǎng)衣在波浪中各個時刻和重力式網(wǎng)箱在來流中的受力和變形；孫宵峰等[12]采用隱式的Newmarkβ方法對網(wǎng)衣的力學(xué)模型進行求解，并建立網(wǎng)衣的可視化工具。

本研究將柔性網(wǎng)衣離散為彈簧和質(zhì)量點的組合體，通過Morison公式計算網(wǎng)線所受水動力，并建立質(zhì)量點的運動方程，基于迭代方法編程求解網(wǎng)衣質(zhì)量點的坐標，最終得到網(wǎng)衣的變形和所受阻力。

1 運動方程的建立

1.1 網(wǎng)繩物理模型

網(wǎng)衣的基本單位為網(wǎng)繩，因此，首先建立網(wǎng)繩的運動方程，并將其拓展以模擬計算三維網(wǎng)衣。在實際應(yīng)用中，組成漁網(wǎng)的網(wǎng)繩大部分都是柔性體，在水流的沖擊下形狀會發(fā)生改變，形狀的改變又會引起作用在網(wǎng)繩上的受力發(fā)生變化，形狀和受力之間密切相關(guān)，相互耦合。為了解決這一耦合問題，將網(wǎng)繩離散成無質(zhì)量的圓柱形彈簧和無體積的質(zhì)量點的結(jié)合體。圓柱體在水中的受力不僅與水流的相對速度也與圓柱體的方位有關(guān)，這里對該模型做出如下假設(shè)：（1）網(wǎng)繩只能承受拉力，不能承受壓力；（2）張力在整個橫截面上為常數(shù)；（3）網(wǎng)繩在拉伸過程中橫截面積保持不變；（4）兩個質(zhì)量點之間的圓柱呈直線，不能發(fā)生彎曲。

完全浸沒在水中的網(wǎng)繩受到的重力和浮力的合力FB可以表示為：

其中，L為網(wǎng)線長度，dw為網(wǎng)線的直徑，ρnet和ρ分別指網(wǎng)線的密度和水的密度。

網(wǎng)線的直徑通常不超過1 cm，相比常見海浪的波長是小量，因此可以利用Morison公式計算網(wǎng)線的水動力。Morison公式為半經(jīng)驗半理論公式，將靜止在水中的構(gòu)件的水動力分為拖曳力和慣性力兩部分：

其中，F(xiàn)D為拖曳力，F(xiàn)I為慣性力，CD為拖曳力系數(shù)，A為沿來流方向的投影面積，CM為慣性力系數(shù)，V為來流的速度。

當構(gòu)建在水中運動時，根據(jù)Brebbia和Waler[13]的建議，Morison公式可以寫為：

其中，Cm=CM-1為附加質(zhì)量力系數(shù)，VR為水質(zhì)點相對于構(gòu)件的移動速度，R為構(gòu)件的位移。

對于聚乙烯網(wǎng)衣，利用Wilson[14]推導(dǎo)的關(guān)于網(wǎng)繩伸長量和張力之間的關(guān)系，網(wǎng)繩的拉力為：

其中，T為網(wǎng)線所受到的拉力，L0為網(wǎng)繩的原始長度，L為網(wǎng)繩變形后的長度，彈性系數(shù)C1和C2的選取參考Gerhard Klust[15]在《纖維繩索》中彈力與變形量之間的實驗數(shù)據(jù)，對于聚乙烯材料網(wǎng)線彈性系數(shù)C1=345.37×106，C2=1.012 1。

對于尼龍網(wǎng)衣，采用H.Le Dret[16]在EZ 20 Lloyd拉力機上對直徑為4 mm的尼龍網(wǎng)線進行實驗并根據(jù)實驗結(jié)果擬合得到的拉力與相對伸長量之間的關(guān)系式：

1.2 運動方程的建立

將網(wǎng)繩的質(zhì)量和所受的重力、浮力、水動力以及拉力分配到兩端的質(zhì)量集中點上，如圖1和圖2所示，在每個質(zhì)量點上建立運動方程。圖2中Bi、Wi分別表示分配在第i個質(zhì)量點的重力和浮力，Ti1和Ti2分別表示質(zhì)量點i所受質(zhì)量點i-1和i+1的拉力，F(xiàn)i表示分配在第i個質(zhì)量點的水動力。

圖1 網(wǎng)繩模型示意圖

圖2 模型質(zhì)量點受力分析示意圖

根據(jù)牛頓第二定律，每一段網(wǎng)繩在水中的受力和加速度可以表示為：

教師要根據(jù)高職學(xué)生的實際學(xué)情，指導(dǎo)學(xué)生參與選擇實用性強、富有趣味的傳統(tǒng)文化內(nèi)容，因材施教，靈活運用情境教學(xué)法、討論法、朗讀法、質(zhì)疑探究法、多媒體教學(xué)法等多種教學(xué)方法，營造出傳統(tǒng)文化學(xué)習(xí)的氛圍，帶領(lǐng)高職學(xué)生走進傳統(tǒng)文化情境中，增進學(xué)生對有關(guān)傳統(tǒng)文化的理解，促進學(xué)生人文素養(yǎng)的提高。例如在教學(xué)《鴻門宴》這一文章時，可以用多媒體播放作者的生平介紹及課文中出現(xiàn)的劉邦、項羽、張良、范增等人物生平經(jīng)歷、古代座次文化等知識，激發(fā)學(xué)生參與語文學(xué)習(xí)活動的興趣，讓學(xué)生更深刻地理解和掌握有關(guān)文學(xué)、文化知識，促進了高職學(xué)生人文素養(yǎng)的不斷提升。

其中，M為網(wǎng)繩的質(zhì)量，將式（3）代入（6）中，可以得到：

若將離散了的每一段網(wǎng)繩的水動力、重力和浮力平均分配到兩端的質(zhì)量點上在加上拉力，如圖2所示，得到的第i個質(zhì)量點的運動方程為：

其中，Mi為分配在第i個質(zhì)量點的質(zhì)量，ΔMi為分配在第i個質(zhì)量點的附加質(zhì)量，ai為第i個質(zhì)量點的加速度，下標x、y、z分別表示在大地坐標系下x、y、z三方向的分量。

2 編程計算過程

Matlab軟件是一款數(shù)學(xué)軟件，有著強大的矩陣運算功能和簡單的語法結(jié)構(gòu)，使用Matlab軟件編程來模擬繩子的受力變形過程，過程如下：

（1）輸入繩子的參數(shù)：繩總長la、繩子直徑dw，分段數(shù)N，網(wǎng)線材料彈性系數(shù)C1，C2，網(wǎng)線材料密度ρnet，流速V、密度ρ、阻力系數(shù)CD和附加質(zhì)量力系數(shù)Cm，計算步長δt，最終迭代進度tend。

（2）將整個繩子劃分為單獨的點，并將這些點的位置、速率、加速度進行初始化。同時，建立時間迭代步的位置、速度和加速度矩陣，并進行初始化。

（3）依照網(wǎng)繩質(zhì)量點i當前的位置，計算出第i個點在第j次迭代中，在大地坐標系下的受力大小。

（4）由步驟（3）得到的受力大小進一步計算出節(jié)點加速度aix、aiy、aiz；根據(jù)加速度和速度計算出在該時間段內(nèi)繩子的質(zhì)點的運動速度Uix、Uiy、Uiz；根據(jù)加速度和速度計算出網(wǎng)繩質(zhì)量點的位移Rix、Riy、Riz，并更新質(zhì)量點i的位置。

（5）網(wǎng)衣的計算結(jié)果收斂時質(zhì)量點受力達到平衡狀態(tài)，加速度趨于0，中間迭代過程的加速度并不影響最終的計算結(jié)果，所以兩次迭代計算的加速度的相對誤差小于0.1時，可認為當前位置已經(jīng)足夠準確，進入步驟（6）進行下一節(jié)點的迭代計算，否則，令j=j+1，進入（3）繼續(xù)循環(huán)迭代。

（6）若i是最后一個點質(zhì)量點則進入下一整體迭代步（t=t+δ）的計算，否則令i=i+1，j=0進入（c）中繼續(xù)進行下一節(jié)點的迭代計算。

（7）重復(fù)（3）至（6）的步驟直至計算到預(yù)先設(shè)定的迭代進度tend。

在步驟（4）中采用NewMark-β方法以增加計算的穩(wěn)定性。NewMark-β方法是對物體運動的加速度做出線性變化的假定，采用該方法后速度Ui和位移Ri的基本方程分別為：

其中，Ri、Ui，ai分別表示質(zhì)量點i的位移、速度和加速度。令α=0.5，β=0.25這相當于加速度在δt內(nèi)為常量，其值為兩個迭代步（t，t+δt）兩端加速度的平均值。

網(wǎng)衣是由縱橫的網(wǎng)線構(gòu)成，三維網(wǎng)衣的模型與網(wǎng)繩離散不同的是，質(zhì)量點（i，j）周圍有4個與之相連的圓柱桿件（i，j-1），（i，j+1），（i+1，j），（i-1，j）。將每個圓柱桿件的質(zhì)量平均分配到兩端的質(zhì)量點上，并且質(zhì)量點承擔(dān)了圓柱桿件所受的水動力、浮力、重力以及拉力，質(zhì)量點所在位置即為網(wǎng)衣目腳交點的所在位置，如圖3所示。網(wǎng)衣質(zhì)量點的運動方程與式（7）相同，編程求解過程與網(wǎng)線一致。

圖3 網(wǎng)衣離散示意圖

當流速大于0.5 m/s時，網(wǎng)繩軸線方向的阻力系數(shù)CDτ和垂直于軸線方向的阻力系數(shù)CDn的選取參考Casarella[17]的實驗結(jié)果，分別取0.8和1.8。沿網(wǎng)線軸線方向的附加質(zhì)量力系數(shù)Cmτ和垂直于軸線方向的阻力系數(shù)Cmn選取，參考趙云鵬[3，18]的研究成果，分別取0.0和1.0。

3 網(wǎng)衣的模擬

3.1 四周固定網(wǎng)衣的模擬

建立單位面積（1 m×1 m）無結(jié)尼龍網(wǎng)衣，命名為S0，S0四邊固定，網(wǎng)目邊長為0.0625 m，網(wǎng)目數(shù)為16×16，網(wǎng)線直徑為0.002 m，密度為1.14 g/cm3。袁軍亭[11]在其博士論文中將網(wǎng)衣離散成集中質(zhì)量和彈簧的組合體，建立平面網(wǎng)衣數(shù)學(xué)力學(xué)模型3DNNM，采用Runge-Kutta方法求解運動方程，對網(wǎng)衣S0模擬計算。本小節(jié)將采用迭代求解的方法，在不同的流速下，計算出網(wǎng)衣S0的變形和受到的水動力以及網(wǎng)線中的拉力，并對比兩種計算方法的結(jié)果。

設(shè)置計算的時間步長δt=10-4，t=48后結(jié)果收斂，這里的時間t并不是指物理時間，而是代表網(wǎng)衣質(zhì)量點的迭代進度。結(jié)果如表1所示，其中，r為相對誤差，T為網(wǎng)線中的最大拉力，F(xiàn)為網(wǎng)片所受的總水動力，X為節(jié)點的最大偏移量，下標ite和3D分別代表本文迭代法和3DNNM方法的計算結(jié)果。

表1 不同流速下兩種計算結(jié)果的比較

網(wǎng)線之間的最大拉力，網(wǎng)衣所受的水動力以及網(wǎng)衣中的節(jié)點的最大偏移量隨著流速的增加不斷增大，不同流速下網(wǎng)衣的變形如圖4所示，網(wǎng)衣的最大變形量出現(xiàn)在中間節(jié)點。兩種計算方法相比，最大拉力之間的相對誤差在10%以內(nèi)，總水動力之間的相對誤差在5%以內(nèi)，最大偏移量之間的相對誤差在10%以內(nèi)，當流速大于0.8 m/s后相對誤差在5%以內(nèi)，表明通過迭代法求解網(wǎng)衣的變形的結(jié)果是可靠的。

3.2 四角固定網(wǎng)衣的模擬

建立模擬網(wǎng)衣S1，S1四角固定，方形網(wǎng)目大小為4 cm（網(wǎng)目目腳長度為4 cm），網(wǎng)線直徑為1.18mm，網(wǎng)衣邊長為80 cm×32 cm，網(wǎng)目數(shù)量為20×8，材料為聚乙烯。對S1進行模擬，設(shè)置流速分別為1.0 m/s、1.5 m/s和2.0 m/s。計算中迭代步長δt=2×10-5，t=50后質(zhì)量點的坐標收斂，停止計算，計算結(jié)果如表2所示。結(jié)果表明，流速越大，網(wǎng)衣中網(wǎng)線的最大拉力、網(wǎng)衣所受的總水動力以及網(wǎng)衣的最大偏移量越大，不同流速下網(wǎng)衣的變形如圖5所示。不論是網(wǎng)衣S0還是網(wǎng)衣S1，流速越大，網(wǎng)衣的變形量越大，但表1和表2中Fite/V2的數(shù)值均越小。根據(jù)網(wǎng)衣阻力系數(shù)CDn的計算公式（11），網(wǎng)衣的面積An、水流的密度ρ和網(wǎng)線的阻力系數(shù)在模擬中保持恒定，表明影響網(wǎng)衣阻力系數(shù)的是形變量，且網(wǎng)衣的形變量越大其阻力系數(shù)越小。

圖5 不同流速下S1的變形

表2 不同流速下四角固定網(wǎng)衣S1計算結(jié)果

4 結(jié)語

通過將網(wǎng)繩離散成無質(zhì)量的彈簧和質(zhì)量點的模型，利用Morison公式求解網(wǎng)線在來流中的水動力，將網(wǎng)線的受力分配集中到質(zhì)量點上，在質(zhì)量點上建立運動方程。采用迭代方法求解各個的運動方程，并利用NewMark-β算法使迭代過程更加穩(wěn)定。將網(wǎng)線的求解模型拓展到三維網(wǎng)衣，在四周固定網(wǎng)衣S0模擬計算中，與通過Runge-Kutta求解運動方程的3DNNM模型的結(jié)果相比，網(wǎng)線之間的拉力和網(wǎng)衣中最大的變形量的相對誤差均在10%以內(nèi)，網(wǎng)衣的水動力的相對誤差在5%以內(nèi)，表明本文通過迭代法求解網(wǎng)衣的變形的方法是可靠的，具有工程應(yīng)用價值。在四角固定網(wǎng)衣S1的計算結(jié)果中發(fā)現(xiàn)，網(wǎng)衣的網(wǎng)繩最大拉力、總水動力和最大形變量隨著流速的提高都逐漸增大。模擬結(jié)果表明，不論是四周固定網(wǎng)衣S0還是四角固定網(wǎng)衣S1的阻力系數(shù)均隨著網(wǎng)衣變形量的增大而逐漸減小。