楊 嵩， 朱先勇， 王 輝， 于 萍

(吉林大學 機械科學與工程學院，長春 130025)

射流沖擊是自然界中廣泛存在的一種自然現(xiàn)象，波浪對礁石及船體的沖擊等都可以理解為射流沖擊[1]。射流沖擊同時也具備廣泛的工程應用背景，從19世紀首次應用于礦物開采到如今廣泛應用于航空、建筑、化工、機械加工、醫(yī)藥等領域，如利用液體射流沖擊進行表面清洗、金屬切割、巖石破碎等[2]；射流沖擊在給工業(yè)界帶來便利的同時也引發(fā)出一系列的問題，如液滴沖擊高速旋轉的渦輪葉片引發(fā)葉片表面磨損及損傷，雨滴對高速飛行器表面蒙皮產(chǎn)生的破壞[3]。因此有必要將射流沖擊尤其是水射流沖擊的研究引入到流體力學的日常教學中。動量定律是流體力學基本定律，作為射流沖擊現(xiàn)象的基本原理，有廣泛而深刻的理論及實踐應用意義[4]。利用水射流沖擊固體方法對動量定律進行驗證也是流體力學實驗教學中常基本方法。

水射流沖擊固體是一個涉及水錘力、沖擊力(滯止壓力)、材料變形、損傷等諸多因素的非線性動力學過程，具有大變形、高應變率、瞬態(tài)過程難于觀察和檢測等特點[5]。國內(nèi)外學者分別對水射流沖擊過程和動量定律實驗裝置進行了理論和應用研究。Field從理論、實驗和應用三個方面對液體沖擊引發(fā)的侵蝕損傷進行了研究，Obara等利用高速相機對高速水射流沖擊PMMA(Polymethyl Methacrylate)樣件過程中射流沖擊波的形成、界面反射及PMMA的損傷進行了研究，重慶大學煤礦災害動力學與控制國家重點實驗室對使用脈沖水射流破巖問題進行了深入的研究[5-6]；廖華林等[7-9]對高壓射流沖擊巖石過程中涉及的流-固耦合作用、巖石破碎機理及在石油工程中的應用進行了研究。國內(nèi)對動量定律實驗的研究多從動量方程本身入手，通過理論計算及實測等方法討論動量系數(shù)等對測量精度的影響并提出解決方法，如丁彤[10]對影響動量定律實驗精度的參數(shù)進行了靈敏度分析，提出解決辦法并開發(fā)了驗證樣機；劉銀慶等[11]對動量定律驗證儀測量部分進行改進，有效減小了紊流、彈簧異常波動對精度的影響；胡衛(wèi)紅等[12]和陳國玉[13]分別設計了一種新式動量定律驗證儀。

作為以射流沖擊現(xiàn)象為基礎的動量定律實驗裝置，射流沖擊的引入有助于提升動量定律實驗裝置工作特性分析的可信度、準確性，也便于在日常的動量實驗教學中更好的向學生闡明射流沖擊原理、過程及實驗誤差產(chǎn)生的原因。

綜上，本文將以吉林大學實驗教學所用射流沖擊型動量定律實驗裝置為研究對象，本裝置為平板射流沖擊式，射流沖擊到平板上，通過測定平板所承受的沖擊力來驗證動量定律。本文擬采用SPH-FEM(Smoothing Particle Hydrodynamics, SPH )耦合算法模擬射流沖擊平板過程，分析平板、傳感器結構在射流沖擊載荷作用下的響應特性，結合分析結果展開誤差討論。

1 射流沖擊型動量定律實驗裝置結構

射流沖擊型動量定律實驗裝置由如下部分組成：①水箱；②動力裝置(泵-電機模塊)；③流量計量裝置；④射流發(fā)生器；⑤圓形平板；⑥底座；⑦檢測裝置(傳感器和數(shù)顯表)，具體如圖1(a)所示。射流沖擊過程所涉及的結構如圖1(b)所示，其中包含圓形平板、電阻式應變傳感器、底座及若干連接件；圓形平板通過中心孔與傳感器固連，傳感器通過螺栓與底座固連，底座通過粘膠與實驗裝置基體固連。圓形平板和傳感器是核心組件，圓形平板用于改變射流動量形成射流沖擊力，傳感器用于傳遞、測量射流沖擊力的大小。

圖1 射流沖擊測試裝置結構圖Fig.1 The structure graph of jet impact test apparatus

2 SPH-FEM耦合算法

水射流沖擊固體結構的過程是一個典型的流-固耦合問題，其中涉及流體的大變形，如果采用傳統(tǒng)的拉格朗日法計算水射流沖擊，容易出現(xiàn)網(wǎng)格畸形等錯誤從而導致計算困難甚至終止[14]。本文采用SPH-FEM耦合算法對水射流沖擊過程進行模擬，此算法由SPH算法和FEM算法兩部分組成，水射流用SPH算法模擬，固體結構采用FEM方法模擬，充分考慮水射流與被沖擊結構之間的作用關系。SPH算法是一種無網(wǎng)格方法，在處理材料大變形、不連續(xù)介質(zhì)動力學問題上優(yōu)勢明顯，廣泛應用于大變形沖擊動力學問題的求解，但SPH算法在求解材料小變形問題的求解精度及計算效率均遜于傳統(tǒng)有限單元法。因此可以將SPH算法與FEM(Finite Element Method)算法結合，充分利用兩種算法的優(yōu)勢，求解大變形、沖擊、小變形等共存耦合問題[15]。

2.1 SPH計算理論

SPH的全稱為光滑粒子流體動力學方法，SPH方法是一種無網(wǎng)格、純拉格朗日方法，其最初被應用于天體物理領域并取得成功。1994年Monaghan首次將SPH方法應用于自由表面流動模擬，隨后SPH方法被廣泛的應用于水力學、沖擊碰撞、自由表面流動等相關的學科研究中[16-19]。

SPH算法的基本邏輯是用一系列任意分布的粒子來代替連續(xù)介質(zhì)流體，通過粒子集合和插值核函數(shù)來估算N-S方程的空間函數(shù)及其導數(shù)實現(xiàn)基本方程→計算公式的轉化，將原來同時含有時間和空間導數(shù)的偏微分方程轉化成只含有時間導數(shù)的方程。SPH方法的基礎是核函數(shù)，通過對核函數(shù)的積分場函數(shù)近似，其表達式如式(1)所示，在此基礎上應用粒子對核近似方程進行再近似，通過應用局部區(qū)域內(nèi)相鄰粒子對應的值來疊加和取代函數(shù)及其導數(shù)積分形式，稱為粒子近似，如式(2)。從而實現(xiàn)流體力學方程的粒子近似過程[20-21]。

(1)

(2)

式中：f為場函數(shù)；Ω為計算區(qū)域；x為坐標向量;h為光滑長度；W為核函數(shù)；ρ為密度；m為質(zhì)量。

將SPH法應用于N-S方程，具體如下：質(zhì)量守恒方程(式(3))，動量守恒方程(式(4))，能量守恒方程(式(5))。

(3)

(4)

(5)

SPH算法處理液-固沖擊問題的優(yōu)點在于對流項直接通過粒子的運動來模擬，完全消除了自由界面上的數(shù)值發(fā)散問題，保證了自由界面追蹤的清晰準確；無網(wǎng)格特性同時避免了網(wǎng)格扭曲及重構帶來的問題。

2.2 SPH-FEM耦合處理

SPH與FEM的耦合是通過接觸罰函數(shù)的方法將質(zhì)點的力作用于有限元單元表面。接觸類型為點-面接觸，SPH粒子為從節(jié)點，F(xiàn)EM單元為接觸主面，在每個時間步內(nèi)均需要檢查SPH粒子與有限單元表面的穿透狀態(tài)，如果存在穿透則應用罰函數(shù)，如果無穿透則不進行任何處理。

3 計算模型描述

結合實驗臺使用實際，本文擬分析射流射速為17.7 m/s條件下的沖擊過程；射流發(fā)生器產(chǎn)生直徑為8 mm的水射流；射流沖擊位置為平板中心，靶距為2 mm，忽略重力；沖擊作用時間設定為30 ms。

實驗裝置的計算模型如圖2所示。模型中各部分的材料、單元屬性如表1所示。材料PMMA和STEEL的基本屬性[22]如表2所示，射流(WATER)狀態(tài)方程EOS選擇US-UP，具體參數(shù)如表2所示。依據(jù)實驗條件，經(jīng)計算可知，射流沖擊力小于20 N，結構所屬材料PMMA和STEEL均處于彈性變形區(qū)，無塑性變形。

圖2 射流沖擊計算模型圖Fig.2 The numerical simulation graph of jet impact model

編號名稱材料單元類型單元尺寸/mm1底座PMMASOLID22傳感器STEELSOLID13圓形平板PMMASOLID54射流WaterParticle1

表2 PMMA，STEEL，WATER材料屬性Tab.2 The properties of material PMMA STEEL and WATER

4 結構位移響應特性分析

結構響應特性分析主要分析傳感器及平板結構在射流沖擊作用下的位移響應，目的在于還原射流沖擊過程中平板及傳感器的運動過程，為誤差分析提供參考。

4.1 傳感器位移響應特性

計算模型中可以將傳感器結構簡化為懸臂梁模型，懸臂梁自由端承載。為簡化描述，結合傳感器應變片粘貼位置，在傳感器結構側面取5個測量點A，B，C，D，E，利用測量點間接描述傳感器結構的響應，A～E位于一條直線上，具體位置如圖3所示。

圖3 傳感器測量點位置分布圖Fig.3 The distribution graph of observation points on sensor

A～E測量點在X，Y，Z方向的位移變化如圖4～圖6所示。由圖4和圖5可知，在X方向(射流沖擊方向)和Y方向上，A～E點發(fā)生了移動且運動趨勢保持一致。A點位移幅值SA最大，E點位移幅值SE最小，B，C，D點位移幅值介于A，E之間，這表明傳感器在X，Y方向上受力并產(chǎn)生了彎曲變形；但傳感器在X，Y方向的具體形變過程及變形幅度略有不同。

由圖4可知，傳感器在X方向受到射流沖擊力作用產(chǎn)生移動和變形，A～E點的位移曲線經(jīng)過振蕩調(diào)整后在t=8 ms時趨于穩(wěn)定，這表明傳感器在X方向的位移和變形在經(jīng)歷振蕩調(diào)整后于t=8 ms時趨于穩(wěn)定，傳感器在X方向的受力Fsx趨于穩(wěn)定。

圖4 傳感器測量點在X方向位移變化圖Fig.4 The displacement graph of sensor observation points with time on direction of X

由圖5可知，傳感器在Y方向受到外力作用產(chǎn)生移動和變形，A～E在Y方向的位移曲線持續(xù)振蕩無穩(wěn)定趨勢，但呈一定周期性。這表明傳感器在Y方向的位移和變形持續(xù)變化，傳感器在Y方向受力Fsy非恒定且保持振蕩。對比X，Y方向上A～E點的位移差值可知，X方向上A～E的位移差值大于Y方向上A～E的位移差值，即傳感器在X方向的變形大于Y方向變形，這表明傳感器在X方向所受外力Fsx大于Y方向所受外力Fsy。

圖5 傳感器測量點在Y方向位移變化圖Fig.5 The displacement graph of sensor observation points with time on direction of Y

由圖6可知，A～E點在Z方向上發(fā)生運動，位移曲線持續(xù)振蕩，A～E點Z方向的位移變化規(guī)律與X,Y方向略有不同，B,E兩點的位移幅值最大、方向相反，C,D兩點次之、方向相反，A點位移最小，造成此種情況的原因有如下兩個：①傳感器在X,Y方向的變形間接引起Z方向位移變化;②傳感器受到了不規(guī)則的彎、扭作用,Z方向的位移變化值遠小于其他兩個方向。

圖6 傳感器測量點在Z方向位移變化圖Fig.6 The displacement graph of sensor observation points with time on direction of Z

通過對傳感器位移響應特性的分析可知，在射流沖擊力的作用下，傳感器在射流方向(X方向)發(fā)生移動，經(jīng)振蕩調(diào)整后趨于穩(wěn)定，傳感器在Y，Z方向由于受到非恒定外力的作用而持續(xù)振動，Y方向的振動類似于彎曲，Z方向振動為不規(guī)則振動，Y方向的振蕩幅度大于Z方向，Z方向的振蕩頻率大于Y方向。Y，Z兩方向上的不規(guī)則振動會對沖擊力的測量帶來誤差。

4.2 平板位移響應特性

計算模型中，平板在射流沖擊作用下的位移響應云圖如圖7所示。在射流沖擊瞬間，平板中心區(qū)域在沖擊力的作用下沿射流沖擊方向發(fā)生移動，由于慣性作用平板的位移分布同幾何結構一致，呈同心圓形式，中心位置位移量最大，最外側位移最小，徑向位置相同區(qū)域的位移保持一致，此時平板發(fā)生凹變形，如圖7(a)所示。隨著沖擊的持續(xù)，平板各點位移量均增加，但分布規(guī)律保持不變，仍為同心圓，位移分布關于Y，Z軸對稱，平板發(fā)生凹變形，如圖7(b)所示。

圖7 平板位移變化圖Fig.7 The displacement graph of plate with time

隨著沖擊的進一步作用，平板各點位移量繼續(xù)增加，但分布規(guī)律發(fā)生了變化，由同心圓形式轉化成僅沿Z軸對稱，最大位移區(qū)域由平板的中心位置沿Z+方向逐漸上移直至頂端，位移分布由圓形變?yōu)閹l狀，以中心孔為界，Z-方向的位移小于等徑向位置Z+方向的位移，這表明平板發(fā)生傾斜并產(chǎn)生彎曲變形，如圖7(c)～圖7(g)所示。t=10 ms以后，平板的位移趨于穩(wěn)定但依舊存在小幅振蕩波動，此時平板的位移分布規(guī)律與之前不同，最大位移區(qū)位于平板的中上段，呈帶條狀分布，平板下半部分的最小位移區(qū)與最大位移區(qū)關于中心孔呈對稱分布，平板位移量的徑向不一致表明平板發(fā)生較嚴重的不規(guī)則彎曲變形。

通過對射流沖擊過程中平板位移響應特性進行分析可知，在射流沖擊作用下，平板中心首先沿射流方向發(fā)生移動并帶動其他區(qū)域運動，隨著射流沖擊的持續(xù)作用，平板的運動由沿射流方向的運動和不規(guī)則的擺動兩部分組成， 沿射流方向的運動經(jīng)振蕩調(diào)整后趨于穩(wěn)定，而不規(guī)則的擺動會持續(xù)存在。平板的不規(guī)則擺動作為一個附加載荷作用在與之相連接的傳感器上，進而產(chǎn)生一個沖擊力測量誤差。

5 傳感器應力特性分析

傳感器應力特性主要分析傳感器結構的應力響應特性，尤其是A～E測量點的應力響應特性。

當射流與圓形平板接觸的瞬間，射流動量發(fā)生突變，形成一個類脈沖沖擊力，平板中心部分沿射流方向發(fā)生移動，平板與傳感器的接觸，產(chǎn)生一個接觸應力，在射流沖擊力作用下，傳感器發(fā)生彎曲產(chǎn)生彎曲應力，如圖8(a)所示。隨著射流沖擊的持續(xù)，傳感器彎曲變形增加，傳感器內(nèi)部彎曲應力增加，傳感器中心的薄壁結構處為一個彎曲應力集中區(qū)，中心薄壁結構抗彎截面系數(shù)小是產(chǎn)生應力集中的主要原因，如圖8(b)和圖8(c)所示。

圖8 沖擊過程中傳感器Mises應力分布圖Fig.8 The Mises contour of sensor in impacting process with time

對比圖8(a)和圖8(d)可發(fā)現(xiàn)，傳感器結構存在另外一個應力集中區(qū)域，即傳感器與平板的連接處。由上文分析可知，在射流沖擊作用下，平板存在不規(guī)則擺動，平板和傳感器之間會產(chǎn)生接觸應力，且最大接觸應力點非固定。

傳感器測量點A～E的Mises應力變化曲線如圖9所示。由圖9可知，A～E測量點的Mises應力持續(xù)振蕩變化，Mises應力的變化趨勢與位移變化趨勢相對應；C點的Mises應力最大，因為C點位于抗彎截面系數(shù)較小的中心薄壁位置，A，B兩點的Mises應力較小，因為A，B兩點更靠近懸臂梁結構的自由端，D，E兩點的Mises應力介于中間，因為D，E兩點位于抗彎截面系數(shù)較大的固定端。

圖9 傳感器測量點Mises應力變化圖Fig.9 The Mises stress graph of sensor observation points with time

選取測量點C為研究對象，對其σx,σy,σz,σxy,σxz,σyz變化情況進行研究，具體如圖10所示。

圖10 傳感器測量點C應力變化圖Fig.10 The stress graph of sensor observation point C with time

由圖10(a)和圖10(b)可知，傳感器C測量點同時存在拉/壓應力和剪應力，即傳感器同時承受彎曲和扭轉作用，且彎曲和扭轉作用是振蕩變化的，這與上文位移響應分析結果相對應。

綜上可知，射流沖擊過程中，射流對平板和傳感器的作用力為非恒定的振蕩變化量；在此外力作用下，平板和傳感器發(fā)生不規(guī)則的擺振運動。在射流沖擊力作用下，傳感器承受彎、扭聯(lián)合作用，傳感器內(nèi)部同時存在彎曲和剪切兩種應力。

6 傳感器測量誤差分析

由上文分析可知，射流沖擊載荷作用下，平板和傳感器會產(chǎn)生不規(guī)則振動。電阻式應變傳感器作為測量元件，其不規(guī)則振動會導致傳感器產(chǎn)生橫向變形，依據(jù)應變傳感器的橫向效應理論，橫向變形會引起測量誤差[23]。由電阻式應變傳感器的測量原理可知，傳感器應變片的彈性變形量s與沖擊力F存在映射關系，即F=f(s)，因此可通過傳感器應變片的彈性形變量s來間接表征射流沖擊力F。

結合傳感器的布置形式及上文分析，利用A～E測量點相對位移來表示應變片的形變量。E點位于固定端，E點的相對位移為0，其他點的相對位移以E點為參考；Z方向的位移變化較小，此處不予考慮。傳感器在X和Y方向相對位移及變形如圖11～圖13所示。

圖11 傳感器測量點X方向相對位移變化圖Fig.11 The relative displacement graph of sensor observation points with time on direction of X

由圖11可知，在射流沖擊作用下，傳感器在X方向產(chǎn)生運動經(jīng)振蕩調(diào)整后趨于穩(wěn)定，A，B兩測量點的相對位移變化曲線并不平滑，存在突變，此突變由圓形平板運動調(diào)整造成，由于A點更接近自由端，A的變化量更大。

由圖12可知，傳感器在射流沖擊力的作用下，在Y方向進行往復振蕩運動，A～E的振蕩頻率相同、振幅依次逐漸降低，振蕩周期隨沖擊作用時長逐漸增加。

圖12 傳感器測量點Y方向相對位移變化圖Fig.12 The relative displacement graph of sensor observation points with time on direction of Y

以A～E測點相對位移數(shù)據(jù)為基礎對傳感器主體部分在X，Y方向上隨時間的變形進行3D擬合，如圖13所示。

圖13 傳感器變形圖Fig.13 The deformation graph of sensor with time

由以上分析可對射流沖擊作用下傳感器所測沖擊力的數(shù)值做如下推測，即傳感器所測定沖擊力的值應該按照c(t)+Asin(wt+φ)規(guī)律進行變化。c(t)由傳感器應變片在X方向的變形決定，經(jīng)過振蕩調(diào)整后為一個穩(wěn)定值；Asin(wt+φ)主要由應變片在Y方向的變形決定，A，w，φ為變量，將A/c定義為位移響應比。若以測量點C的相對位移來表征傳感器的測量誤差變化，考慮應變傳感器橫向效應引發(fā)的測量誤差。由電阻應變式傳感器的測量原理可知，應變片在二向應變場下的響應為

(6)

式中：R為電阻值；Sx為縱向靈敏度；Sy為橫向靈敏度；Sα為剪切靈敏度；εx為縱向應變；εy為橫向應變；εα為剪切應變。

一般情況下剪切靈敏度較小，可以忽略。此時式(6)變?yōu)?/p>

(7)

式中：k=Sy/Sx為橫向靈敏度系數(shù)，此處可假定εy=-μ1εx，μ1為傳感器結構材料的泊松比。對式(7)進行變換可得

(8)

式中:Sg為廠家測定的靈敏度系數(shù)，對比式(7)和式(8)可得

Sg=Sx(1-μ1k)

(9)

綜合式(7)～式(9)可得

(10)

經(jīng)變換可得

(11)

如果傳感器僅考慮應變片的靈敏度系數(shù)，則有

(12)

(13)

傳感器橫向效應所引發(fā)的誤差為

(14)

若選取傳感器測量點C在Y方向上相對E點的位移作為橫向效應的邊界值。此時可認為，測量過程中Y和X方向上的位移響應比的絕對值|A/c|與應變比εy/εx相等，即εy/εx∈[0%,7%]。傳感器結構材料的泊松比為0.3，即μ1=0.3，對k取不同值，可估算出傳感器的測量誤差如表3所示。由表3可知，傳感器橫向效應所引發(fā)的測量誤差與橫向靈敏度系數(shù)k、應變比εy/εx有關，傳感器的橫向靈敏度系數(shù)可以認為是傳感器的固有特性參數(shù)，此時傳感器的測量誤差是一個關于應變比εy/εx的一次函數(shù)，依據(jù)式(14)可知，εy/εx越大，傳感器橫向效應所引發(fā)的測量誤差越大。

表3 傳感器誤差估算表Tab.3 Sensor estimated error

依據(jù)仿真參數(shù)對動量定律實驗裝置進行實測，測試數(shù)據(jù)如表4所示，實測沖擊力的誤差范圍小于1%，與依據(jù)式(14)估算的誤差(見表3)相近，此結果表明，文中所采用的測量誤差估算方法針對此問題有效。造成估算誤差與實測誤差之間差異的原因如下：

① 傳感器內(nèi)部存在誤差補償；

② 溫度、濕度等引起的傳感器測量誤差；

③C點的位移響應不能完全反應應變片的真實變化；

④ 實測過程中水泵流量脈動引發(fā)的射流速度變化；

⑤ 計算模型、材料屬性等與真實系統(tǒng)之間的差異。

表4 實驗誤差表Tab.4 Experiment error

7 結 論

本文以射流沖擊型動量定律實驗裝置為研究對象，采用SPH-FEM方法分析了射流沖擊載荷作用下，圓形平板及傳感器結構的位移響應特性和傳感器的應力特性，以此為基礎對傳感器測量誤差產(chǎn)生的原因及誤差范圍進行了分析、計算，并與實測實驗結果進行對比，驗證分析、計算方法的可靠性。綜合以上分析獲得如下結論：

(1)射流沖擊型動量定律實驗裝置在工作過程中，圓形平板及傳感器結構在非射流沖擊方向上會發(fā)生持續(xù)不斷的不規(guī)則振動，此不規(guī)則振動引發(fā)傳感器產(chǎn)生測量誤差。

(2)應變電阻的橫向理論可用于應變傳感器的測量誤差估計。

(3)不規(guī)則振動引發(fā)的測量誤差的大小與傳感器位移響應比(應變比)呈線性函數(shù)關系，位移響應比越大，測量誤差越大。

(4)射流沖擊力的測量值由固定值和變化值兩部分組成，固定值主要由傳感器在射流沖擊方向的變形決定，變化值由傳感器在其他方向上振動引發(fā)的變形決定。

本文所采用的分析方法及所獲得的結論已部分應用于本科生動量方程實驗教學環(huán)節(jié)，為學生理解實驗誤差的產(chǎn)生原因提供幫助。