張琴林 杜海 孔文杰 劉文靜 譚周杭 王磊
摘要:無操縱面飛行器是近年來航空領(lǐng)域研究熱點(diǎn),其主要實(shí)現(xiàn)方式為吹氣射流形成環(huán)量,產(chǎn)生所需控制力矩。為研究環(huán)量控制技術(shù)應(yīng)用于無舵面飛行器進(jìn)行偏航姿態(tài)控制的效果,本文對(duì)環(huán)量控制機(jī)翼在不同迎角下,開展了吹氣系數(shù)分別為0,0.005,0.010,0.020,0.030的風(fēng)洞試驗(yàn)及數(shù)值模擬研究。風(fēng)洞試驗(yàn)結(jié)果表明,環(huán)量控制器在較小的吹氣系數(shù)下仍能產(chǎn)生較大的偏航力矩,同時(shí)隨著吹氣系數(shù)的增加,阻力系數(shù)減小,偏航力矩隨之增大;數(shù)值模擬結(jié)果顯示隨著吹氣系數(shù)的增加,向后的射流動(dòng)量增大,產(chǎn)生推力效果,從而產(chǎn)生偏航力矩。
關(guān)鍵詞:環(huán)量控制;飛行器;偏航;風(fēng)洞試驗(yàn);數(shù)值模擬
中圖分類號(hào):V211.7文獻(xiàn)標(biāo)識(shí)碼:ADOI:10.19452/j.issn1007-5453.2020.05.010
基金項(xiàng)目:國家自然科學(xué)基金(51806181)
飛行器的飛行控制,主要通過升降舵、副翼和方向舵等舵面來分別控制飛行器的俯仰、滾轉(zhuǎn)及偏航[1]。而這些活動(dòng)舵面的存在也帶來了一系列的問題,如復(fù)雜繁重的機(jī)械操縱系統(tǒng)、頻繁的檢查維修、隱身性能的降低等,這些都會(huì)影響飛行器的整體性能。若采用無舵面的飛行控制方式,如利用環(huán)量控制[2]技術(shù)來實(shí)現(xiàn)對(duì)飛行器的俯仰、滾轉(zhuǎn)以及偏航姿態(tài)的控制,將會(huì)使飛行器獲得一些優(yōu)勢,如提高飛行器的氣動(dòng)性能、短距起降性能、隱身性能等[3]。
從20世紀(jì)60年代以來,國內(nèi)外研究人員利用風(fēng)洞試驗(yàn)及數(shù)值模擬方法對(duì)環(huán)量控制技術(shù)開展了大量的研究,并嘗試在飛行器上運(yùn)用[4-5]。2005年,環(huán)量控制第一次應(yīng)用于固定翼無人機(jī)Irvine Tutor 40,實(shí)現(xiàn)了滾轉(zhuǎn)姿態(tài)控制[6]。2010年英國BAE系統(tǒng)公司(BAE System)聯(lián)合多所高校開發(fā)了Demon無人驗(yàn)證機(jī),其在試飛期間的滾轉(zhuǎn)姿態(tài)操縱可以通過環(huán)量控制技術(shù)實(shí)現(xiàn)[1]。大量資料顯示[7-11],將環(huán)量控制技術(shù)應(yīng)用于飛行器的滾轉(zhuǎn)姿態(tài)控制已取得突破性進(jìn)展,而將其應(yīng)用于偏航姿態(tài)控制上的研究還很少。
本文將對(duì)后緣采用環(huán)量控制技術(shù)的機(jī)翼進(jìn)行偏航特性研究,在20m/s風(fēng)速、不同迎角下,對(duì)不同吹氣系數(shù)下的吹氣環(huán)量控制機(jī)翼的升阻力系數(shù)以及偏航力矩進(jìn)行對(duì)比分析,同時(shí)利用數(shù)值模擬方法來研究機(jī)翼后緣射流對(duì)偏航力矩的控制機(jī)理。
1環(huán)量裝置
1.1環(huán)量噴口參數(shù)
環(huán)量控制技術(shù)一般由固定的幾何升力面和圓形后緣組成,在升力面和圓形后緣之間(即接近后緣的位置)開口,在翼型內(nèi)部空腔形成高壓氣流并通過開口沿著物面切向產(chǎn)生射流,射流與外流混合后沿著彎曲的圓形后緣表面形成科恩達(dá)效應(yīng)(Coanda Effect)[12],在表面處,氣流受到的壓力與離心力達(dá)到平衡,如圖1所示。環(huán)量控制可以有效延遲邊界層的分離,增加附體流動(dòng)的距離,改變翼型前后緣的駐點(diǎn)位置和環(huán)量,提高升力[13],同時(shí)射流與外流混合,帶動(dòng)外流使流線發(fā)生較大的偏折,產(chǎn)生類似于襟副翼的效果。吹氣射流的動(dòng)量系數(shù)[14]定義為:
影響環(huán)量控制裝置性能的因素有很多,主要有射流出口高度h、后緣半徑RCC、射流出口展向速度均勻性等。為保證環(huán)量控制的高效性,在綜合了國內(nèi)外相關(guān)研究[14-16]的基礎(chǔ)上,本文設(shè)計(jì)了一套環(huán)量控制裝置,參數(shù)如下:射流出口高度h=0.3mm,科恩達(dá)(Coanda)尾緣半徑RCC=4.5mm。
1.2吹氣裝置
本文所設(shè)計(jì)的環(huán)量控制裝置由導(dǎo)流裝置、射流噴口和科恩達(dá)后緣組成。其中試驗(yàn)?zāi)P蜌庠床捎?2MPa的高壓氣體,從模型外部引氣,提供高壓高速氣流。
為了實(shí)現(xiàn)與舵面一樣的功能,使氣流發(fā)生上下偏轉(zhuǎn),將噴管從進(jìn)氣口開始設(shè)計(jì)成上下獨(dú)立的兩部分,由中部的一塊平板隔開,上下兩部分內(nèi)腔體積相等,如圖2所示。當(dāng)上翼面噴口吹氣時(shí),氣流將沿著科恩達(dá)后緣向下偏轉(zhuǎn),形成正環(huán)量;當(dāng)下翼面噴口吹氣時(shí),氣流將沿著科恩達(dá)后緣向上偏轉(zhuǎn),形成負(fù)環(huán)量,從而產(chǎn)生飛行器姿態(tài)控制所需的氣動(dòng)力和氣動(dòng)力矩。為了使噴出的氣流均勻,導(dǎo)流裝置做了氣動(dòng)優(yōu)化,裝置內(nèi)部設(shè)置了導(dǎo)流片,導(dǎo)流片呈曲線狀將內(nèi)腔分為6個(gè)區(qū)域,導(dǎo)流片的厚度為1mm,在接近進(jìn)氣口處設(shè)有倒圓角。整個(gè)機(jī)翼內(nèi)部導(dǎo)流裝置如圖3所示。
2風(fēng)洞試驗(yàn)
2.1風(fēng)洞試驗(yàn)設(shè)備
本試驗(yàn)在西華大學(xué)XHWT風(fēng)洞進(jìn)行,風(fēng)洞為低速回流式風(fēng)洞,試驗(yàn)段尺寸為3m×1.2m×1.2m(長×寬×高),可控制風(fēng)速為0.5~60m/s,收縮比為7.1,湍流強(qiáng)度ε≤0.5%,氣流偏角|Δα|≤0.50,|Δβ|≤0.50。本次試驗(yàn)所采用的模型具體尺寸為:環(huán)量控制機(jī)翼模型展長為725.721mm,后掠角為39.959°,弦長為320.842mm,其中環(huán)量控制噴口出氣寬度為LCC= 488.928mm,環(huán)量控制機(jī)翼實(shí)物圖如圖4所示。
2.2氣動(dòng)測量系統(tǒng)
本試驗(yàn)采用盒式六分量應(yīng)變天平,天平響應(yīng)頻率大于80Hz,總體尺寸為200mm×100mm×60mm的矩形結(jié)構(gòu),該盒式天平六分量量程和檢驗(yàn)精準(zhǔn)度見表1。
整個(gè)天平測力系統(tǒng)由天平本體、高精度直流電源、精密信號(hào)放大器、東華DH8300N動(dòng)態(tài)信號(hào)采集卡,以及采集處理軟件構(gòu)成。高精度直流電源為天平提供直流電源,天平受力產(chǎn)生的微弱測量信號(hào),經(jīng)過信號(hào)放大器放大后輸送到東華采集系統(tǒng),再經(jīng)過采集處理軟件的采集和處理,將各通道電信號(hào)轉(zhuǎn)換為模型各分量的氣動(dòng)力和氣動(dòng)力系數(shù)。
3試驗(yàn)結(jié)果
3.1重復(fù)性試驗(yàn)與結(jié)果
本次研究前,開展了風(fēng)洞標(biāo)模試驗(yàn),其結(jié)果和理論結(jié)果一致,表明測力系統(tǒng)的可靠性。為了測試系統(tǒng)的重復(fù)性,選定來流風(fēng)速為20m/s,不同迎角下測量半模機(jī)翼模型氣動(dòng)力及氣動(dòng)力矩,重復(fù)4次。測量得到,升力系數(shù)的偏差量在±0.002~0.008之內(nèi),標(biāo)準(zhǔn)差為0.016,重復(fù)性良好。滾轉(zhuǎn)力矩系數(shù)的偏差量在±0.006~0.02之內(nèi),標(biāo)準(zhǔn)差為0.04。由升力系數(shù)及滾轉(zhuǎn)力矩系數(shù)可知,該天平測量試驗(yàn)臺(tái)測量的精確性良好,能夠進(jìn)行模型氣動(dòng)力和氣動(dòng)力矩的測量。
3.2偏航控制參數(shù)及結(jié)果
常規(guī)飛行器的偏航運(yùn)動(dòng)是通過方向舵的偏轉(zhuǎn)來實(shí)現(xiàn)的,而基于環(huán)量控制的無舵面飛行器取消了常規(guī)活動(dòng)舵面,則通常利用向后吹氣來實(shí)現(xiàn)偏航控制。
3.2.1偏航控制參數(shù)
本次試驗(yàn)是通過流量計(jì)來測量吹氣裝置的噴氣量,因此需將吹氣系數(shù)轉(zhuǎn)換為吹氣流量。其中,吹氣系數(shù)Cμ與吹氣流量的關(guān)系見式(2):
式中:Q為吹氣流量,V∞為來流速度,c為機(jī)翼弦長,L為機(jī)翼展長,hj為射流出口高度,Lj為環(huán)量控制噴口出氣寬度,Cμ為射流動(dòng)量通量與自由來流動(dòng)量通量的比值(即吹氣系數(shù))。
為了實(shí)現(xiàn)無滾轉(zhuǎn)的偏航控制,固定下噴口Q2吹氣量不變,調(diào)節(jié)上噴口Q1吹氣流量來匹配下吹氣量。本次試驗(yàn)開展了下吹氣系數(shù)Cμ分別為0(即無射流),0.005,0.010,0.020,0.030的偏航控制研究。
上、下噴口對(duì)應(yīng)的吹氣流量如圖5所示(如Cμ= 0.005即表示Q2對(duì)應(yīng)的吹氣系數(shù)為0.005)。
3.2.2偏航控制結(jié)果
本次試驗(yàn)中通過同時(shí)開啟上、下兩側(cè)噴氣裝置(上、下兩側(cè)噴口吹氣流量如圖5所示),保持滾轉(zhuǎn)力矩一致,來測量吹氣裝置的升阻力系數(shù)以及偏航力矩系數(shù)。
圖6是來流風(fēng)速20m/s,副翼位置環(huán)量控制裝置上、下同時(shí)吹氣(上、下噴口吹氣量見圖5)升力系數(shù)曲線,其升力系數(shù)完全重合,實(shí)現(xiàn)了滾轉(zhuǎn)力矩一致的前提。
圖7是來流風(fēng)速20m/s,開啟環(huán)量控制裝置上、下兩個(gè)噴口(上、下噴口吹氣量見圖5),在不同吹氣系數(shù)Cμ下所產(chǎn)生的阻力系數(shù)隨迎角的變化情況。通過曲線可以看出,在不同吹氣系數(shù)Cμ下,阻力系數(shù)曲線趨勢一致。在-4°~32°迎角范圍內(nèi),隨著吹氣系數(shù)Cμ的增大,阻力系數(shù)不斷減?。ㄗ⒁獾阶枇ο禂?shù)隨迎角變化規(guī)律和常規(guī)翼型不一致,其原因是吹氣模型外接吹氣管路的影響)阻力系數(shù)隨吹氣系數(shù)Cμ增大而減小的原因是:吹氣裝置上噴口和下噴口協(xié)同吹氣的射流向后產(chǎn)生了推力。
圖8是來流風(fēng)速20m/s,副翼位置環(huán)量控制裝置上、下同時(shí)吹氣(上、下噴口吹氣量如圖5所示)的偏航力矩系數(shù)曲線??梢钥闯?,在-4°~ 0°迎角范圍內(nèi),偏航力矩系數(shù)變化量較小。在迎角為0°~20°之間時(shí),隨著迎角的增大,偏航力矩系數(shù)隨之增大。隨著吹氣系數(shù)Cμ的增大,偏航力矩也會(huì)隨之增加,當(dāng)下噴口吹氣系數(shù)Cμ=0.03時(shí),偏航力矩控制達(dá)到最優(yōu)。
這一現(xiàn)象充分表明機(jī)翼后緣上、下噴口同時(shí)吹氣可以增加偏航控制效果。其內(nèi)在機(jī)理是:當(dāng)上、下兩個(gè)吹氣縫同時(shí)吹氣時(shí),兩股射流會(huì)發(fā)生混合,混合氣流產(chǎn)生一個(gè)類似于推力的反作用力,并且這種作用力效果隨著吹氣系數(shù)的增加而增加,從而產(chǎn)生偏航力矩。
4偏航力矩控制機(jī)理數(shù)值模擬研究
為了進(jìn)一步探究射流動(dòng)量吹氣系數(shù)對(duì)飛行器偏航特性的影響,從而清晰直觀地揭示出環(huán)量控制裝置使飛行器產(chǎn)生偏航力矩的作用機(jī)理,本次對(duì)2D環(huán)量控制機(jī)翼進(jìn)行了數(shù)值模擬,具體過程如下。
4.1幾何模型
數(shù)值模擬所用翼型為Clark-Y翼型,翼型的尖后緣修形為科恩達(dá)表面形狀,射流出口高度及后緣半徑與試驗(yàn)?zāi)P捅3忠恢?。流域設(shè)計(jì)如圖9所示,二維機(jī)翼位于流場中心,弦長為c,距離上、下邊界均為15c,為保證邊界不影響機(jī)翼周圍流場,流場入口距離機(jī)翼前緣為15c,同時(shí)出口距離機(jī)翼后緣為30c,以保證流場充分發(fā)展。
4.2控制方程
二維黏性不可壓牛頓流體運(yùn)動(dòng)的基本控制方程為連續(xù)性方程和Navier-Stokes方程:
4.3計(jì)算網(wǎng)格劃分
為了減少計(jì)算網(wǎng)格數(shù)量,同時(shí)保證計(jì)算精度,本次采用多域網(wǎng)格的劃分方法,即對(duì)機(jī)翼附近流域進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分,其網(wǎng)格尺寸相對(duì)較小,在外部流域進(jìn)行非結(jié)構(gòu)網(wǎng)格劃分且增大網(wǎng)格尺寸。在ICEM中劃分的計(jì)算網(wǎng)格總數(shù)為205448,網(wǎng)格劃分如圖10所示。
4.4邊界條件及計(jì)算模型
本文用FLUENT數(shù)值模擬軟件對(duì)環(huán)量控制2D機(jī)翼進(jìn)行數(shù)值模擬,采用有限體積法對(duì)控制方程進(jìn)行離散,同時(shí)由于S-A模型對(duì)壁面邊界的空氣流動(dòng)問題以及邊界層中具有逆壓梯度的問題處理結(jié)果較好,本次湍流模型采用S-A模型。
其自由來流入口及射流噴口邊界條件均設(shè)置為壓力入口(對(duì)于噴口而言,即為吹氣系數(shù)分別為0,0.005,0.010,0.020,0.030時(shí)對(duì)應(yīng)的動(dòng)壓)。出口邊界條件為壓力出口,壁面采用無滑移邊界條件。
4.5數(shù)值驗(yàn)證
為了檢驗(yàn)數(shù)值模擬結(jié)果的正確性,本次利用前文中所得到的試驗(yàn)結(jié)果作為數(shù)值模擬的參照對(duì)象。在來流風(fēng)速為20m/s、0°迎角下,后緣吹氣系數(shù)Cμ分別為0、0.005,0.010,0.020,0.030時(shí)(此處的上、下噴口吹氣比例與前文試驗(yàn)保持一致),通過數(shù)值模擬計(jì)算所得到的升阻力系數(shù)與相同條件下試驗(yàn)所得到的升阻力系數(shù)進(jìn)行比較,如圖11、圖12所示。從圖中可以看出,數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)基本吻合。
4.6數(shù)值模擬結(jié)果
圖13是來流風(fēng)速為20m/s、0°迎角下,吹氣系數(shù)Cμ(上、下兩噴口同時(shí)吹氣)分別為0,0.005,0.010,0.020,0.030時(shí)的流線速度云圖,從圖中可以看出,當(dāng)兩個(gè)出氣縫同時(shí)吹氣時(shí),兩股射流會(huì)發(fā)生混合以及抵消。且當(dāng)上、下兩噴口射流動(dòng)量吹氣系數(shù)接近相同時(shí),混合射流不發(fā)生上下的偏轉(zhuǎn),而表現(xiàn)為沿機(jī)翼弦線向后噴射,隨著吹氣系數(shù)的增大,向后的射流動(dòng)量增加,從而對(duì)機(jī)翼產(chǎn)生一個(gè)更強(qiáng)的向前反推力,根據(jù)力的平衡條件,相當(dāng)于減小了機(jī)翼上向后的流動(dòng)阻力。
對(duì)于飛行器而言,若左右兩側(cè)機(jī)翼未同時(shí)向后吹氣或兩側(cè)機(jī)翼吹氣系數(shù)不相同時(shí),在左右兩側(cè)機(jī)翼上便會(huì)出現(xiàn)推力不相等,從而產(chǎn)生繞飛行器重心的偏航力矩,這與風(fēng)洞試驗(yàn)的測力結(jié)果相一致。因此驗(yàn)證了環(huán)量控制器代替?zhèn)鹘y(tǒng)控制舵面的可行性。
5結(jié)論
本文通過風(fēng)洞試驗(yàn)及數(shù)值模擬對(duì)環(huán)量控制機(jī)翼的偏航特性進(jìn)行了探究,得到以下結(jié)論:
(1)在來流風(fēng)速為20m/s,吹氣系數(shù)分別為0,0.005,0.010,0.020,0.030時(shí),環(huán)量控制機(jī)翼的阻力系數(shù)會(huì)隨著吹氣系數(shù)的增加而減小,而偏航力矩系數(shù)隨吹氣系數(shù)的增加而增大。當(dāng)吹氣系數(shù)為0.030時(shí),阻力系數(shù)為最小,偏航力矩最大。
(2)當(dāng)環(huán)量控制器上、下兩個(gè)噴口同時(shí)向后吹氣時(shí),兩股射流會(huì)發(fā)生混合以及抵消,表現(xiàn)為沿機(jī)翼弦線向后噴射,隨著吹氣系數(shù)的增大,向后的射流動(dòng)量增加,產(chǎn)生向前的反推力,從而實(shí)現(xiàn)偏航姿態(tài)控制。
參考文獻(xiàn)
[1]陳黎,常亮.英國惡魔無人驗(yàn)證機(jī)及射流飛控技術(shù)發(fā)展分析[J].飛航導(dǎo)彈,2011(10):56-59. Chen Li, Chang Liang. Development analysis of British Demon unmanned aerial vehicle and jet flight control technology [J].Aircraft Missiles, 2011(10):56-59.(in Chinese)
[2]Liu Y,Sankar L N,Englar R J,et al. Numerical simulations of the steady and unsteady aero-dynamic characteristics of a circulation control airfoil[R]. AIAA-2001-0704,Reston,VA:AIAA,2001:2-4.
[3]張攀峰,燕波,戴晨峰.合成射流環(huán)量控制翼型增升技術(shù)[J].中國科學(xué):技術(shù)科學(xué), 2012, 42(9): 3-5. Zhang Panfeng, Yan Bo, Dai Chenfeng. Lift enhancement method by synthetic jet circulation control[J]. Chinese Science: Technology Science, 2012, 42(9): 3-5.(in Chinese)
[4]Maskell E C,Spence D A. A theory of the jet flap in three dimensions[J].ProceedingsoftheRoyalSocietyA:Mathematical,Physical and Engineering Sciences,1959,51(1266):420-423.
[5]Nielson J N,Biggers J C. Recent progress in circulation control aero-dynamics[R]. AIAA-1987-0001,Reston,VA:AIAA,2002.
[6]Sparks T,Michie S,Gill K,et al. Development of an integrated circulation control/fluidic thrust vectoring flight test demonstrator[C]//1st International Conference on Innovation and Integration in Aerospace Sciences,Belfast,Queens University,2005.
[7]王海洋.基于環(huán)量控制的無人飛行器氣動(dòng)特性研究與飛行試驗(yàn)[D].南京:南京航空航天大學(xué),2014. Wang Haiyang. Aerodynamic characteristics research and flight test of UAV based on circulation control [D]. Nanjing: Nanjing University ofAeronautics andAstronautics, 2014 .(in Chinese)
[8]Yao J,Qiu J,Ji H,et al. Numerical investigation on the unsteady aerodynamic characteristics of flexible variablesweep morphing aircraft[J]. Acta Aerodynamica Sinica,2015,33(3):406-413.
[9]Wetzel D,Griffin J,Cattafesta L N. An experimental study of circulation control on an Elliptic Airfoil[R]. AIAA-2009-4280,2009.
[10]Fielding J P,Lawson C P,Martins Pires R M,et al. Design,build and flight of the Demon demonstrator UAV[C]// 11th AIAAATIO Conference,2011.
[11]Novak C,Cornelius K,Roads R . Experimental investigations of the circular wall jet on a circulation control airfoil [C]// 25th AIAAAerospace Sciences Meeting,1987.
[12]Zhang Yanhua,Zhang Dengcheng,Li Lin. Experimental study on aerodynamic properties of circulation control airfoil with plasma jet[C]//Asia-Pacific International Symposium on Aerospace Technology,Singapore,2018.
[13]Hard S L. Lift augmentation using a NACA 0018 airfoil with circulation control[D]. ProQuest,UMI Dissertations Publishing,2010.
[14]Kanistras K,Saka P C,Valavanis K P,et al. Wind tunnel investigation of a circulation control wing with Dual-Radius flaps[J]. Journal ofAircraft,2018,11(1):1-3.
[15]Swanson R C,Rumsey C L . Computation of circulation control airfoil flows[J]. Computers & Fluids,2009,38(10):1925-1942.
[16]張旺龍.低雷諾數(shù)下附面層分離的抽吸控制及優(yōu)化[D].南京:南京理工大學(xué), 2014.Zhang Wanglong. Suction control and its optimization of boundary layer separation at low Reynolds numbers[D]. Nanjing: Nanjing University of Science and Technology, 2014.(in Chinese)(責(zé)任編輯陳東曉)
作者簡介
張琴林(1997-)女,碩士研究生。主要研究方向:實(shí)驗(yàn)空氣動(dòng)力學(xué)。
Tel:18483678195
E-mail:1015882732@qq.com;
杜海(1985-)男,講師,博士。主要研究方向:空氣動(dòng)力學(xué)及飛行器總體設(shè)計(jì)。
Tel:15196686983
E-mail:duhai@mail.xhu.edu.cn
Study on Control of Yaw Moment by Trailing Edge Blowing
Zhang Qinlin,Du Hai*,Kong WenJie,Liu WenJing,Tan Zhouhang,Wang Lei Xihua University,Chengdu 610039,China
Abstract: In recent years, the aircraft without control surface is a hot research topic in the field of aviation. The main way to achieve this is to form a circulation by blowing jet and generate the required control moment. In order to study the effect of circulation control technology applied to yaw attitude control of flapless aircraft, this paper studies the wind tunnel test and numerical simulation of the circulation control wing with blowing coefficients of 0, 0.005, 0.010, 0.020 and 0.030 at different angles of attack. The wind tunnel test results show that the circulation controller can still produce a large yaw moment with a small blowing coefficient, and with the increase of blowing coefficient, the drag coefficient decreases, and the yaw moment increases. The results of numerical simulation show that with the increase of blowing coefficient, the momentum of backward jet increases, thus resulting in thrust effect and yaw moment.
Key Words: circulation control; aircraft; yaw; wind tunnel experiment; numerical simulation