鞋涉水問題的多體動力學與流-固耦合系統(tǒng)建模及仿真分析

吳曉建，祁祺?，王斌，周兵

（1.南昌大學機電工程學院，江西，南昌，330031；2.湖南大學汽車車身先進設(shè)計制造國家重點實驗室，湖南長沙 410082）

CAE（Computer Aided Engineering，CAE）技術(shù)在工程分析中發(fā)揮著重要的作用，它為工程問題的發(fā)現(xiàn)、解決提供了關(guān)鍵性的驗證和指導.然而，復雜系統(tǒng)往往同時涉及多體動力學與流-固耦合等不同技術(shù)方向，有效而準確的聯(lián)合仿真已成為工程界共性技術(shù)難點.馬增輝等[1]采用ALE（Arbitrary Lagrange-Euler，ALE）法進行了水陸兩用飛機降落耐波性的研究；江旭東等[2]采用耦合歐拉-拉格朗日CEL（Coupled Euler-Lagrangian，CEL）方法研究了壓差驅(qū)動式管道機器人的動力特性；雷娟棉等[3]采用CFD（Computational Fluid Dynamics，CFD）的6 自由度動網(wǎng)格法分析了導彈掛載分離過程；楊曉彬等[4]采用CFD 的重疊網(wǎng)格法分析了飛機水上降落的水動力砰擊載荷.以上典型復雜工程應用問題側(cè)重于流-固耦合系統(tǒng)建模與仿真，未涉及與多體動力學的聯(lián)合.本文所研究的行人在雨天行走時的鞋涉水問題則是一個多體多力學與流-固耦合仿真的問題，不同的步態(tài)和鞋在涉水過程具有差異化的“涉水-出水”過程，準確地將此過程予以模擬呈現(xiàn)對工程中復雜耦合問題的仿真分析具有良好的參考價值.

鞋的涉水過程指在積水路面行走時，鞋在人體下肢運動的驅(qū)動下和水進行接觸的過程.整個過程同時涉及結(jié)構(gòu)體、空氣和水三相的相互作用，包括流動介質(zhì)突變引起的載荷突變、自由液面變化、姿態(tài)變化等眾多的復雜問題.

同時，鞋的運動往往與人體的下肢運動相關(guān)，其姿態(tài)以及運動形式往往是復雜多變的，進而導致涉水運動成為復雜的多體動力學-流-固耦合問題.

本文將行人行走的鞋涉水問題分析分為兩個基本方面，首先建立人體下肢模型，采用多體動力學分析，得到鞋體運動的邊界條件；然后建立鞋涉水多體運動的流-固耦合模型.其中，耦合歐拉-拉格朗日CEL 方法，采用基于體積分數(shù)的流-固耦合邊界追蹤算法，自動在結(jié)構(gòu)和流體域間進行載荷、位移、速度等信息的傳遞，能夠解決空化射流沖擊[5]、存在大變形的切削問題[6]以及震蕩過程的流固耦合問題[7].因而本文將采用CEL 方法對鞋“涉水-出水”過程的流-固耦合問題進行模擬.

最后，基于以上研究方法，對比分析運動鞋和平底鞋的“涉水-出水”過程，通過涉水過程時間和所需能量等，評價兩種鞋涉水性能的強弱.本文以生活中常見的行人在濕地行走，其鞋涉水現(xiàn)象為研究對象，開展復雜系統(tǒng)多體動力學-流固耦合仿真技術(shù)研究，盡管不是直接的工程問題，但對類似工程問題的建模及分析具有良好的借鑒作用.

1 鞋涉水過程的分析方法

1.1 基于Adams 的人體步態(tài)運動

采用Adams 多體動力學仿真軟件模擬人體行走步態(tài)，進行足部動力學和運動學分析.通過把上身視為一個質(zhì)量塊，下肢部分分成大腿、小腿和足部，以及在各關(guān)節(jié)之間設(shè)立運動副和驅(qū)動，建立人體下肢簡化模型.這種簡化方法可以解決步態(tài)仿真模型不易穩(wěn)定的問題，且較為準確地反映出人體步態(tài)運動.

1.2 耦合歐拉-拉格朗日法

應用于傳統(tǒng)有限元分析的拉格朗日方法，以物體的坐標為基礎(chǔ)、以有限元網(wǎng)格的節(jié)點為物體的質(zhì)點，可以跟蹤質(zhì)點的運動，從而描述物體的運動，但遇到大變形問題時，會由于網(wǎng)格變形扭曲而失去原有的精度，甚至導致計算無法繼續(xù).歐拉方法常用于流體動力學計算，以網(wǎng)格為空間，以網(wǎng)格節(jié)點為空間點，材料不隨單元變形而在網(wǎng)格中流動，因此可以解決包含大變形和流體流動的問題，但在本文鞋涉水過程流固耦合（Fluid-Structure Interaction，F(xiàn)SI）分析中，歐拉分析雖然可有效進行流體流動分析，卻在捕捉結(jié)構(gòu)流固交界面上存在一定困難；并且歐拉方法固存的數(shù)值耗散問題導致其對計算資源和時間要求較高.由此可知，單純的拉格朗日算法和單純的歐拉算法都有各自的缺陷和不足，也有著各自的優(yōu)勢；如果將兩者有機地結(jié)合起來，可以解決一些只用單一方法所不能解決的問題.耦合歐拉-拉格朗日CEL方法就是基于這一目的最早由Noh[8]提出的，它采用有限差分法求解帶有移動邊界的二維流體動力學問題.在Noh 的研究中，網(wǎng)格點可以隨物質(zhì)點一起運動，但也可以在空間中固定不動，甚至網(wǎng)格點可以在一個方向上固定，而在另一個方向上隨物質(zhì)一起運動.因此，CEL 算法在解決物體的大位移時，比如碰撞、流體動力學及流體-固體之間相互作用時有強大的優(yōu)勢.

2 基于Adams 的人體步態(tài)運動建模

2.1 人體下肢運動的幾何模型

2.1.1 動力學建模數(shù)據(jù)

因本研究只關(guān)注人體下肢的運動，故省略人體上肢建模.下肢結(jié)構(gòu)的主要參數(shù)如表1 所示，其中，腳直接使用鞋的幾何模型.

表1 下肢主要參數(shù)Tab.1 Main parameters of human lower limb

2.1.2 模型平衡控制

人體模型是一種非線性且不易穩(wěn)定的系統(tǒng)，如何調(diào)整各關(guān)節(jié)之角度或速度來保持平衡一直是非常復雜的難點問題.本文提出一種簡化方法解決模型運動過程中的平衡問題，即假設(shè)人體運動過程中質(zhì)心位置變化不大，在此條件下胯部的運動和地面成平行狀態(tài)，因此在胯部設(shè)一平衡桿來控制胯部運動實現(xiàn)平衡，如圖1 所示.

圖1 人體行走幾何模型Fig.1 Geometric model of human walking

2.2 人體下肢運動的動力學模型

2.2.1 模型的運動約束

下肢建模是人體能否實現(xiàn)步行及步行質(zhì)量好壞的關(guān)鍵所在，尤以關(guān)節(jié)結(jié)構(gòu)為集中體現(xiàn).人體膝關(guān)節(jié)由脛骨-股骨和髕骨-股骨關(guān)節(jié)組成，這樣膝關(guān)節(jié)就成為脛骨和股骨相咬合的關(guān)節(jié)[9-10]；同時，膝關(guān)節(jié)和髖關(guān)節(jié)都只有一個自由度，可以用轉(zhuǎn)動副連接.將地面設(shè)置為固定副，胯部和平衡桿之間為移動副.在各轉(zhuǎn)動副設(shè)置驅(qū)動器，設(shè)置地面和腳之間的接觸形式為碰撞.創(chuàng)建出符合實際的步態(tài)運動模型，如圖2 所示.

圖2 步態(tài)運動仿真模型Fig.2 Gait motion simulation model

2.2.2 模型的運動控制

人的行走是一個周期活動過程，一個步態(tài)周期為從腳跟著地到腳跟再次著地的時間.如圖3 所示，步態(tài)周期分為支撐相和擺動相.其中，支撐相是指腳跟著地到腳離地，每只腳的支撐相又分為首次觸地、支撐初期、支撐中期、支撐末期4 個部分，占整個步態(tài)周期的62%左右；擺動相是腳尖離地到腳跟著地，分為擺動早期、擺動中期、擺動末期，大約占整個步態(tài)周期38%[11].本文使用step 函數(shù)對關(guān)節(jié)角度進行控制，以模擬行人的步態(tài)運動.

圖3 人體步態(tài)周期Fig.3 Human gait cycle

2.2.3 模型的運動仿真

設(shè)置仿真總時間為2 s，仿真步數(shù)為500，即每個仿真步長為0.004 s，仿真結(jié)果如圖4 所示.

圖4 步態(tài)仿真結(jié)果Fig.4 Gait simulation results

對比圖3 和圖4 可以發(fā)現(xiàn)，所建立的人體步態(tài)仿真模型，其模擬結(jié)果與步態(tài)周期吻合度高，能夠較為準確地反映出行走過程中鞋的運動與姿態(tài).

為進一步開展鞋涉水流-固耦合仿真，提取出鞋參考點的速度和角速度作為有限元仿真的邊界條件，參考點選取了髖關(guān)節(jié)和腳跟連線中點的位置，具體如圖5 所示，參考點的速度和角速度值如圖6 所示.

圖5 參考點的位置Fig.5 Position of reference point

圖6 參考點的速度Fig.6 Reference point speed

3 鞋涉水過程流固耦合仿真

3.1 鞋涉水的系統(tǒng)描述

鞋涉水有限元模型主要包括三部分，即歐拉區(qū)域、鞋體和剛性地面.其中歐拉區(qū)域為歐拉體，鞋體和剛性地面為拉格朗日體.考慮到鞋體-地面-流體之間的復雜的接觸關(guān)系，通過罰函數(shù)法的通用接觸算法來描述多體間的運動.

先以運動鞋為例分析其涉水過程.模型如圖7所示，運動鞋單只質(zhì)量為0.29 kg，鞋長250 mm，鞋底最厚處為30 mm、最薄處8 mm.為減少網(wǎng)格數(shù)量和計算消耗，暫不考慮鞋底花紋影響，且對鞋體模型進行適當簡化.

圖7 運動鞋建模Fig.7 Sports shoes model

將鞋體材料簡化為橡膠來模擬鞋體的彈性行為.通過Mooney-Rivlin 模型描述橡膠材料鞋體模型的超彈性本構(gòu)關(guān)系，則有：

式中：W 表示應變能密度；C10和C01為材料常數(shù)；I1和I2為第一、第二Green 不變量.

Ogden[12]已經(jīng)證明，在應變能密度和應力、應變之間存在如下關(guān)系：

式中：t 為真實應力；λ1為主伸長比.對于單軸拉伸實驗，t2=t3=0，三個方向主伸長λ1=λu，λ2=λ3=1/，則應力應變關(guān)系為：

對于雙軸拉伸實驗，t1=t2，t3=0 三個方向主伸長λ3=，λ1=λ2=λB，則應力應變關(guān)系為：

對于平面拉伸實驗，t3=0，三個方向主伸長λ1=λS，λ2=1，λ3=1/λS則應力應變關(guān)系為：

對橡膠材料試件進行單軸拉伸、雙軸拉伸和平面拉伸實驗，得到橡膠材料的實驗數(shù)據(jù).由此，根據(jù)式（5）、（6）、（7）擬合橡膠材料的拉伸實驗數(shù)據(jù)，結(jié)果如圖8 所示（各條曲線表示材料拉伸數(shù)據(jù)和各種本構(gòu)方案下的擬合曲線）.最終獲得Mooney-Rivlin 模型參數(shù)C10=0.181 MPa，C01=0.003 84 MPa.

圖8 擬合的橡膠應變勢能曲線Fig.8 Fitted strain-stress curve

采用CEL 算法時，為防止歐拉區(qū)域和剛體地面接觸失效，同時保證水花濺起高度的捕捉，需要設(shè)置歐拉區(qū)域遠大于拉格朗日區(qū)域.

在Abaqus 中對水的材料屬性進行特殊定義，具體通過材料的狀態(tài)方程實現(xiàn)，在此選用Mie-Gruneisen 方程[13]，如公式（8）所示：

式中：p 為液體壓強，Γ0、s 為材料常數(shù)；η=1 -ρ0/ρ為名義體積壓縮應變；ρ0為材料初始密度；Em為單位質(zhì)量內(nèi)能；C0為材料初始聲速.

3.2 有限元模型的建立

歐拉區(qū)域的主要尺寸為1 000 mm × 500 mm ×150 mm，為了保證歐拉區(qū)域和剛性地面的成功接觸，設(shè)置接觸區(qū)域尺寸為1 000 mm×500 mm×5 mm.劃分鞋體單元為四結(jié)點線性四面體單元C3D4，地面單元為四結(jié)點三維雙線性剛性四邊形單元R3D4，歐拉區(qū)域為八結(jié)點線性歐拉六面體單元EC3D8R.建立有限元模型如圖9 所示.

圖9 流固耦合過程的有限元模型Fig.9 Finite element model

為防止歐拉材料流出邊界，設(shè)置歐拉區(qū)域邊界的所有法向速度為0.地面為解析剛體，約束其所有自由度.將鞋體表面節(jié)點耦合到鞋的參考點上，將Adams 得到的運動數(shù)據(jù)添加到鞋參考點上，以設(shè)置鞋體的邊界條件，且設(shè)置歐拉區(qū)域、鞋體和剛性地面的接觸為通用接觸.

在CEL 算法中，歐拉區(qū)域及自由表面通過體積分數(shù)描述[14-15].體積分數(shù)1 代表歐拉材料充滿歐拉單元；0 代表歐拉單元中不存在歐拉材料.初始狀態(tài)下歐拉區(qū)域體積分數(shù)都為0，這表示初始狀態(tài)下歐拉區(qū)域沒有歐拉材料分布.為了描述水的初始分布，需設(shè)定歐拉材料在歐拉區(qū)域的初始分布，本文初始分布的體積模型來自于布爾操作，通過體積分數(shù)工具指定其初始材料分布.根據(jù)鞋體的初始位置，得出歐拉模型和初始水域模型如圖10 所示.

圖10 歐拉模型和水域初始節(jié)點分布Fig.10 Euler model and initial node distribution in waters finite element model

3.3 計算結(jié)果分析

經(jīng)過有限元分析計算，圖11 給出了CEL 計算所得鞋涉水過程.選取了幾個比較重要的時刻，可以看出，隨著鞋的運動，CEL 算法可以清楚地計算水的運動過程.

圖11 涉水計算結(jié)果Fig.11 Wading results

進一步地，將鞋“涉水-出水”過程分為5 個階段，如圖11（a）～（e），左側(cè)為水的體積分數(shù)云圖，右側(cè)為鞋體的應力云圖.第一階段鞋體運動使其前端與地面接觸，在這一過程中，發(fā)生了沖擊，使鞋體周圍的水有向外排出的趨勢，為排水階段，如圖11（a）；第二階段鞋體開始向前運動，推鞋前端的水向前運動，為推水階段，如圖11（b）；第三階段，鞋繼續(xù)向前運動，鞋前的水繼續(xù)堆積，在鞋前端產(chǎn)生一個水膜，為水膜生成階段，如圖11（c）；第四階段，鞋繼續(xù)向前運動，并伴隨向上抬起動作，鞋前的水膜開始破裂，為水膜破裂階段，如圖11（d）；第五階段，鞋前的水膜完全破裂，鞋從水中出來，為完全出水階段，如圖11（e）.

讀取水區(qū)域動能變化曲線如圖12 所示.可以看到出現(xiàn)了能量峰值，峰值和陡升起點分別在0.125 s和0.175 s.

圖12 水區(qū)域動能變化圖Fig.12 Water kinetic energy change chart

這兩個時刻的水的體積分數(shù)云圖如圖13（a）所示；水域速度云圖如圖14（a）所示.可以發(fā)現(xiàn)，陡升起點出現(xiàn)在排水階段，峰值出現(xiàn)在水膜發(fā)展階段，這是由于，排水和推水過程使水體發(fā)生大范圍運動，使系統(tǒng)動能陡升.在后續(xù)階段中，因鞋體出水，水和鞋的相對沖擊減小，導致水域動能減小，但仍出了一個峰值，其起始時間和峰值時間分別在0.335 s 和0.365 s.后兩個時刻的體積分數(shù)云圖如圖13（b）所示；水域速度云圖如圖14（b）所示.可以發(fā)現(xiàn)，為了補充因鞋的劃水運動而產(chǎn)生的水坑，水開始向水坑運動；并且這段時間內(nèi)飛濺的水下落到水面上，使系統(tǒng)動能增大.

圖13 峰值時刻的體積分數(shù)云圖Fig.13 Volume fraction at peak time

圖14 峰值時刻的水速度云圖Fig.14 Water velocity at peak time

4 不同鞋的涉水性能比較

選取鞋長和運動鞋相同的平底鞋作為對比，兩者比較明顯的特征區(qū)別在于平底鞋鞋頭前端沒有上翹.為確?？杀刃裕O(shè)置運動鞋底厚度的最厚和最薄尺寸與運動鞋一致，且平底鞋的前端鞋面高度與運動鞋的大致相當.根據(jù)以上原則選擇圖15 所示平底鞋模型并加以適當簡化.

圖15 平底鞋建模Fig.15 Flat shoes model

除鞋體模型外，其他模型和行人步態(tài)等邊界條件均保持不變，通過計算，比較相同時刻水的體積分數(shù)云圖和速度云圖等來分析兩種鞋涉水性能.

圖16 和圖17 分別為兩種模型在主要時刻下的體積分數(shù)云圖和水域速度云圖.時刻1 中運動鞋模型處于水膜破裂階段，平底鞋處于水膜生成階段，可以發(fā)現(xiàn)平底鞋生成的水膜大于運動鞋生成的水膜.時刻2 中運動鞋處于水膜破裂末期，平底鞋處于水膜破裂初期，平底鞋水膜包裹鞋體的面積大于運動鞋.時刻3 中運動鞋處于完全出水階段，平底鞋處于水膜破裂階段，平底鞋濺起的水花大于運動鞋濺起的水花.

圖18 為運動鞋和平底鞋計算模型的系統(tǒng)動能比較，可以發(fā)現(xiàn)平底鞋的水域系統(tǒng)動能大于運動鞋的水域動能，因為步態(tài)數(shù)據(jù)的一致，動能峰值位置基本相同.

圖16 體積分數(shù)云圖結(jié)果比較Fig.16 Comparison of volume fraction results

圖17 速度云圖結(jié)果比較Fig.17 Comparison of velocity results

圖18 不同鞋水區(qū)域系統(tǒng)動能比較Fig.18 Different water kinetic energy comparison

通過對不同時刻下的體積分數(shù)云圖和水域速度云圖分析，可以發(fā)現(xiàn)，在相同運動步態(tài)下，平底鞋濺起的水花大于運動鞋，因而從涉水到完全出水過程所花費的時間大于運動鞋；通過分析不同模型水系統(tǒng)的動能變化，可以發(fā)現(xiàn)平底鞋系統(tǒng)水區(qū)域動能大于運動鞋，因此，平底鞋的涉水-出水過程要花費更多的能量.以涉水時間和消耗能量大小為涉水性能評價指標，平底鞋的涉水性能弱于運動鞋.

5 結(jié)論與展望

1）建立了行人步態(tài)多體動力學仿真模型和鞋涉水流-固耦合模型，采用耦合歐拉-拉格朗日方法完整分析了鞋“涉水-出水”過程，并對瞬態(tài)過程進行了可視化，從而將鞋“涉水-出水”過程分為排水、推水、水膜生成、水膜破裂和完全出水5 個階段.

2）在相同的行人步態(tài)和鞋底高度情況下，對兩種不同的鞋進行涉水仿真，通過對比各個階段的體積分數(shù)、水域速度和水域動能變化，發(fā)現(xiàn)所研究的平底鞋的涉水性能弱于運動鞋.

3）本研究開展的多體動力學-流-固耦合系統(tǒng)仿真，可為類似復雜工程應用問題仿真分析提供技術(shù)參考，也可為鞋的涉水性能分析及鞋的優(yōu)化設(shè)計提供一定的借鑒.

本文下一步可開展的研究工作主要有：

1）不同的步態(tài)可能導致不同的涉水性能，可以選擇多種步態(tài)運動以分析其在鞋涉水過程中的影響.

2）在流-固耦合仿真當中，可考慮鞋的底部花紋，以及鞋的親水性，由此引起行走時水隨著鞋的運動被澆起從而將鞋頭打濕，探究不同的鞋頭形狀和步態(tài)產(chǎn)生的鞋頭澆水現(xiàn)象的差異化，在此基礎(chǔ)上以進一步對鞋頭設(shè)計給予理論指導.