基于CEL算法的發(fā)動(dòng)機(jī)潤(rùn)滑液瞬態(tài)振蕩過(guò)程可視化研究

馬子遠(yuǎn)， 俞小莉， 黃鈺期， 劉震濤， 黃 瑞

(浙江大學(xué) 能源工程學(xué)院， 杭州 310000)

在內(nèi)燃機(jī)中，油底殼主要作為一個(gè)儲(chǔ)存潤(rùn)滑油的容器，并通過(guò)其吸入系統(tǒng)將潤(rùn)滑油供應(yīng)到發(fā)動(dòng)機(jī)部件上，并直接影響到發(fā)動(dòng)機(jī)的工作性能。油底殼中的總油量是根據(jù)發(fā)動(dòng)機(jī)中潤(rùn)滑旋轉(zhuǎn)部件所需的油流量及滿足整個(gè)潤(rùn)滑系統(tǒng)循環(huán)所需的油流量來(lái)確定[1]。一個(gè)通用的發(fā)動(dòng)機(jī)油底殼-吸油-過(guò)濾器系統(tǒng)，如圖1所示。

油底殼設(shè)計(jì)的重要目標(biāo)是既要滿足基于發(fā)動(dòng)機(jī)整體設(shè)計(jì)標(biāo)準(zhǔn)的NVH性能，同時(shí)還要滿足車(chē)輛設(shè)定工況(操作工況)下的潤(rùn)滑性能，這都與油底殼內(nèi)的潤(rùn)滑油振蕩情況有關(guān)[2]。油底殼中潤(rùn)滑油隨機(jī)體振蕩是一個(gè)涉及多體動(dòng)力學(xué)、彈塑性力學(xué)及流體力學(xué)等多學(xué)科的問(wèn)題，也是一個(gè)典型的流固耦合振蕩問(wèn)題。目前隨著數(shù)值計(jì)算技術(shù)的飛速發(fā)展，一些新型的數(shù)值計(jì)算方法(如SPH法、ALE法、混合CFD法等)為解決復(fù)雜的流固耦合問(wèn)題提供了強(qiáng)大的技術(shù)支持，尤其是為振蕩CFD問(wèn)題的深入研究開(kāi)辟了新的方法[3]。

圖1 油底殼吸油-過(guò)濾-殼布局

拉格朗日算法(Lagrangian)以物質(zhì)的坐標(biāo)為基礎(chǔ)，有限元網(wǎng)格的節(jié)點(diǎn)即為物體的質(zhì)點(diǎn)，可以跟蹤質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡，從而能準(zhǔn)確地描述物體邊界的運(yùn)動(dòng);但是在涉及大變形的問(wèn)題時(shí)，單元網(wǎng)格將會(huì)出現(xiàn)嚴(yán)重的畸變現(xiàn)象，其可能最終會(huì)導(dǎo)致無(wú)法繼續(xù)計(jì)算。歐拉方法(Eulerian) 常用于流體力學(xué)計(jì)算分析，計(jì)算中所采用的網(wǎng)格以空間坐標(biāo)為基礎(chǔ)，有限元節(jié)點(diǎn)即為空間點(diǎn)，因此能夠處理物質(zhì)的扭曲及一些大變形問(wèn)題，但由于歐拉法在捕捉物體邊界信息上較為困難，不能夠精確描述物質(zhì)的邊界，并且該方法固存的數(shù)值耗散問(wèn)題導(dǎo)致其對(duì)計(jì)算資源和時(shí)間要求較高[4]。單純的拉格朗日和單純的歐拉算法都有各自的缺陷和不足，但是又有著各自的優(yōu)勢(shì)。如果將兩者有機(jī)地結(jié)合起來(lái)，可以解決一些只用單一方法所不能解決的問(wèn)題。耦合的歐拉-拉格朗日(Coupled Lagrangian Eulerian，CEL) 方法就是基于這一目的最早由Noh[5]提出的，并采用有限差分法求解了帶有移動(dòng)邊界的二維流體動(dòng)力學(xué)問(wèn)題。在Noh的研究中，網(wǎng)格點(diǎn)可以隨物質(zhì)點(diǎn)一起運(yùn)動(dòng)，但也可以在空間中固定不動(dòng)，甚至網(wǎng)格點(diǎn)可以在一個(gè)方向上固定，而在另一個(gè)方向上隨物質(zhì)一起運(yùn)動(dòng)。因此，CEL算法在解決物體的大位移時(shí)，比如碰撞、流體動(dòng)力學(xué)及流體-固體之間的相互作用時(shí)有強(qiáng)大的優(yōu)勢(shì)。圖2為L(zhǎng)agrangian算法、Eulerian算法及CEL算法在相同時(shí)間Δt上的物體形態(tài)與空間網(wǎng)格變形狀態(tài)示意圖。

(a)t時(shí)刻 (b)t+Δt時(shí)刻

圖2 3種算法對(duì)應(yīng)的網(wǎng)格變形特點(diǎn)

Fig.2 Characteristics of deformed meshes of Lagrangian, Eulerian and CEL algorithms

耦合歐拉-拉格朗日(CEL)方法對(duì)流固耦合振動(dòng)數(shù)值仿真研究提供了一種新的思路，其本質(zhì)上是以歐拉材料域在受拉格朗日部件運(yùn)動(dòng)影響下的流動(dòng)特性作為研究對(duì)象的，不受網(wǎng)格運(yùn)動(dòng)的限制[6]，因此在振蕩流計(jì)算方面具有獨(dú)特優(yōu)勢(shì)。

鑒于CEL方法在處理此類問(wèn)題方面的優(yōu)勢(shì)，目前國(guó)內(nèi)外已有不少學(xué)者開(kāi)展了耦合歐拉-拉格朗日方法(CEL)的研究，F(xiàn)oucard等[7]采用了CEL算法對(duì)大變形或極端變形情況下的超彈性材料問(wèn)題進(jìn)行了計(jì)算分析，采用擴(kuò)展有限元(XFEM)方法來(lái)對(duì)系統(tǒng)力學(xué)平衡方程和變形梯度張量方程進(jìn)行了離散化，并提出了一種適應(yīng)于材料屬性介于流體和固體之間的物體的運(yùn)動(dòng)的計(jì)算方程；Rostami等[8]采用CEL方法來(lái)模擬兩相納米流(流體中包含AL2O3粒子)，對(duì)其在微通道波紋管內(nèi)的共軛流動(dòng)換熱問(wèn)題進(jìn)行了研究；Diggs等[9]對(duì)采用CEL方法研究多相流問(wèn)題中的體積分?jǐn)?shù)計(jì)算提出了評(píng)價(jià)方法，并對(duì)基于網(wǎng)格法和基于粒子法兩種方法進(jìn)行了對(duì)比研究； Bahremand等[10]采用CEL方法對(duì)定常壁面熱流條件下的螺紋管內(nèi)的納米流湍流流動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬和試驗(yàn)測(cè)試；Vujanovi等[11]采用CEL法對(duì)柴油機(jī)的噴霧、燃燒及污染物的生成過(guò)程進(jìn)行了數(shù)值模擬，并與已有的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比，研究顯示CEL方法在解決此類問(wèn)題上有明顯優(yōu)勢(shì)，同時(shí)具有高的計(jì)算效率和計(jì)算精度；Jarauta等[12]采用嵌入式CEL方法對(duì)高分子聚合物燃料電池電解質(zhì)氣道內(nèi)的水滴動(dòng)態(tài)特性進(jìn)行了數(shù)值研究，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析。

在國(guó)內(nèi)，葉林征等[13]為探究超聲珩磨中不同角度空化微射流沖擊下壁面特性的變化，采用CEL算法進(jìn)行了模型的模擬和分析；卿啟湘等[14]針對(duì)盾構(gòu)機(jī)切刀切削土體仿真收斂困難，仿真結(jié)果震蕩嚴(yán)重等復(fù)雜非線性問(wèn)題，利用CEL法對(duì)模型進(jìn)行了求解，并與理論公式計(jì)算結(jié)果進(jìn)行了對(duì)比，發(fā)現(xiàn)CEL方法在處理此類問(wèn)題時(shí)收斂快，震蕩小，計(jì)算效率高，對(duì)網(wǎng)格大小敏感性低，結(jié)果更接近理論解；姚小虎等[15]采用CEL算法分析了復(fù)雜的水陸兩棲飛機(jī)水上降落時(shí)，機(jī)頭入水瞬時(shí)的結(jié)構(gòu)動(dòng)力響應(yīng)；王建華等[16]根據(jù)模擬鉆井船在黏土層中插樁對(duì)鄰近樁影響的離心模型試驗(yàn)結(jié)果，驗(yàn)證了通過(guò)CEL方法結(jié)合非線性地基梁有限元計(jì)算方法來(lái)分析此類問(wèn)題的可行性。

1 數(shù)值分析模型

1.1 系統(tǒng)描述

分析系統(tǒng)為新設(shè)計(jì)開(kāi)發(fā)的深井后置式4缸柴油發(fā)動(dòng)機(jī)的油底殼系統(tǒng)。系統(tǒng)橫截面，如圖3所示。

圖3 發(fā)動(dòng)機(jī)的橫截面圖

分析的系統(tǒng)模型參數(shù)如表1所示。

油底殼作為拉格朗日部件，潤(rùn)滑油區(qū)域劃分為歐拉網(wǎng)格。采取CEL方法計(jì)算模型時(shí)，歐拉網(wǎng)格立方體模型與油底殼相比要確保足夠大，以保證歐拉材料的自由流動(dòng)和在潤(rùn)滑液的振蕩高度。計(jì)算的材料參數(shù)如表2所示。

表1 算例系統(tǒng)計(jì)算參數(shù)

表2 材料計(jì)算參數(shù)

除上述材料參數(shù)外，還需使用材料的狀態(tài)方程進(jìn)行描述。選用 Mie-Gruneisen 方程，如下所示

p-ph=Γρ(Em-Eh)

(1)

(2)

Eh=(phη)/(2ρ0)

(3)

式中：p為液體壓強(qiáng)；ph為Hugoniot壓強(qiáng)；ρ0為液體初始密度；c0為液體聲速；η為名義體積壓縮應(yīng)變,有η=1-ρ0/ρ,s為待定常數(shù)。Γ=Γ0(ρ0/ρ)為Gruneisen率，Γ0為材料常數(shù),Em為單位質(zhì)量?jī)?nèi)能。

1.2 有限元模型建立

建立的數(shù)值模型，如圖4所示。油底殼劃分為二階四面體單元，C3D10M。歐拉立方體模型劃分為八節(jié)點(diǎn)歐拉單元，EC3D8R。

圖4 在網(wǎng)格狀態(tài)下的裝配剖視圖

在CEL算法中，歐拉材料通過(guò)歐拉體積分?jǐn)?shù)(EVF)來(lái)跟蹤其經(jīng)過(guò)網(wǎng)格的狀態(tài)，所有歐拉單元需通過(guò)指定值來(lái)代表其充滿歐拉材料的比例。針對(duì)這樣的問(wèn)題，在初始條件里必須指定歐拉區(qū)域空間里面被流體材料所占據(jù)的初始體積。體積分?jǐn)?shù)工具執(zhí)行實(shí)體模型和歐拉網(wǎng)格的布爾運(yùn)算。這將創(chuàng)建一組節(jié)點(diǎn)，定義了油底殼中的初始潤(rùn)滑油量狀態(tài)，圖5所示。拉格朗日材料和歐拉材料間的接觸通過(guò)基于罰函數(shù)接觸算法的一般接觸分析來(lái)計(jì)算，當(dāng)歐拉單元中體積分?jǐn)?shù)0時(shí)，代表該歐拉單元中沒(méi)有歐拉材料，拉格朗日單元能沒(méi)有任何阻礙地通過(guò)歐拉單元。

1.3 計(jì)算工況

為獲得對(duì)實(shí)際振蕩情況的估計(jì)，又采用比較低的時(shí)間成本，采用時(shí)間縮放的計(jì)算策略，以減少計(jì)算時(shí)間，同樣也是之所以選擇雙精度ABAQUS/Explicit執(zhí)行計(jì)算的原因。表3為計(jì)算工況。

表3 數(shù)值模擬計(jì)算工況

2 計(jì)算結(jié)果與分析

2.1 啟動(dòng)工況計(jì)算結(jié)果分析

經(jīng)過(guò)有限元計(jì)算分析，圖6給出了采用CEL算法得到的潤(rùn)滑液在曲軸由靜止到啟動(dòng)過(guò)程中潤(rùn)滑油的瞬時(shí)振蕩過(guò)程，選取了其中幾個(gè)典型時(shí)刻，可以看出，隨著曲軸的轉(zhuǎn)動(dòng)，CEL算法可以清楚地計(jì)算得到潤(rùn)滑油隨曲軸轉(zhuǎn)動(dòng)的運(yùn)動(dòng)情況，啟動(dòng)后隨著曲軸轉(zhuǎn)動(dòng)，潤(rùn)滑油向一側(cè)濺起，沖擊油底殼一側(cè)，并向整個(gè)曲軸-油底殼密封空間飛濺。

為進(jìn)一步體現(xiàn)出CEL算法在流固耦合動(dòng)態(tài)振蕩過(guò)程可視化效果方面的優(yōu)勢(shì)，選取某一時(shí)刻，觀測(cè)計(jì)算模型的材料流動(dòng)情況和網(wǎng)格變化。如圖7所示。

調(diào)取液面濺出網(wǎng)格一節(jié)點(diǎn)的體積分?jǐn)?shù)變化曲線，如圖8所示?？梢钥闯鲈擖c(diǎn)的體積分?jǐn)?shù)由1近似突變?yōu)?，再突變回1，正好驗(yàn)證了該處的潤(rùn)滑油由液體濺出再掉落回液體中的過(guò)程。

調(diào)取潤(rùn)滑油動(dòng)能變化曲線，如圖9所示?？梢钥闯鲈摑?rùn)滑油的動(dòng)能能量變化有許多尖峰，說(shuō)明在振蕩過(guò)程中潤(rùn)滑油有能量的波動(dòng)，可以進(jìn)行優(yōu)化系統(tǒng)的NVH性能。

圖10給出了基于材料體積分?jǐn)?shù)加權(quán)平均的應(yīng)力(SVAVG)分布。這種顯示方式在多種材料混合時(shí)能更直觀的顯示不同液體材料的振蕩情況。

(a) 時(shí)刻1

(b) 時(shí)刻2

(c) 時(shí)刻1

圖8 關(guān)注點(diǎn)處體積分?jǐn)?shù)變化情況

2.2 加速制動(dòng)工況計(jì)算結(jié)果分析

圖11是加速工況下計(jì)算時(shí)刻結(jié)束時(shí)潤(rùn)滑油的振蕩情況及網(wǎng)格的變化情況，可以看出，當(dāng)車(chē)輛突然加速時(shí)，由于油底殼的推動(dòng)作用，潤(rùn)滑油會(huì)被推向前，但又由于慣性的作用，會(huì)導(dǎo)致潤(rùn)滑油在油底殼內(nèi)來(lái)回晃動(dòng)，由網(wǎng)格變化中的波紋可以看出，尾部波紋形狀為兩個(gè)凹峰，伴隨著逐漸向前的傳遞波紋。

圖9 系統(tǒng)動(dòng)能變化情況

圖10 加速工況計(jì)算結(jié)果

調(diào)取井槽尾部中間節(jié)點(diǎn)的應(yīng)力變化情況，如圖12所示，可以看出隨著潤(rùn)滑油在油底殼內(nèi)振蕩，油底殼上的沖擊應(yīng)力也隨之不斷變化，出現(xiàn)多次尖峰，這對(duì)整機(jī)的NVH性能也有很大影響。

圖12 井槽尾部中間區(qū)域節(jié)點(diǎn)應(yīng)力變化曲線

圖13是制動(dòng)工況下計(jì)算時(shí)刻結(jié)束時(shí)潤(rùn)滑油的振蕩情況及網(wǎng)格的變化情況，可以看出，當(dāng)車(chē)輛突然制動(dòng)時(shí)，由于油底殼的推動(dòng)作用，潤(rùn)滑油會(huì)被推向前，但又由于慣性的作用，會(huì)導(dǎo)致潤(rùn)滑油在油底殼內(nèi)來(lái)回晃動(dòng)，由網(wǎng)格變化中的波紋可以看出，尾部波紋形狀為兩個(gè)凸峰，伴隨著逐漸向前的傳遞波紋。

調(diào)取井槽尾部中間節(jié)點(diǎn)的應(yīng)力變化情況，如圖14所示，可以看出隨著潤(rùn)滑油在油底殼內(nèi)振蕩，油底殼上的沖擊應(yīng)力也隨之不斷變化，出現(xiàn)多次尖峰，這對(duì)整機(jī)的NVH性能也有很大影響。

圖14 井槽尾部中間區(qū)域節(jié)點(diǎn)應(yīng)力變化曲線

2.3 轉(zhuǎn)向工況計(jì)算結(jié)果分析

圖15是加速工況下計(jì)算時(shí)刻結(jié)束時(shí)潤(rùn)滑油的振蕩情況及網(wǎng)格的變化情況，可以看出，當(dāng)車(chē)輛轉(zhuǎn)彎時(shí)，由于油底殼的推動(dòng)作用，潤(rùn)滑油會(huì)被推向前轉(zhuǎn)彎方向，但又由于慣性的作用，會(huì)導(dǎo)致潤(rùn)滑油在油底殼內(nèi)來(lái)回晃動(dòng)，由網(wǎng)格變化中的波紋可以看出，波紋形狀為向轉(zhuǎn)向側(cè)濺出，伴隨著逐漸兩側(cè)傳遞波紋。

調(diào)取井槽一側(cè)中部區(qū)域節(jié)點(diǎn)的應(yīng)力變化情況，如圖16所示，可以看出隨著潤(rùn)滑油在油底殼內(nèi)振蕩，油底殼上的沖擊應(yīng)力也隨之不斷變化，出現(xiàn)尖峰，這是油底殼上一側(cè)的節(jié)點(diǎn)調(diào)取結(jié)果。由這曲線尖刺可以看出，潤(rùn)滑液振蕩這對(duì)油底殼以及整機(jī)的NVH性能也有很大影響。

圖16 油底殼一側(cè)中間區(qū)域節(jié)點(diǎn)應(yīng)力變化曲線

3 結(jié) 論

流固耦合振動(dòng)問(wèn)題一直以來(lái)是工程領(lǐng)域研究的熱點(diǎn)和難點(diǎn)，它也是系統(tǒng)NVH性能分析上的一個(gè)不可缺少的研究環(huán)節(jié)，耦合的歐拉-拉格朗日算法(CEL算法)，在解決一些復(fù)雜的流固耦合問(wèn)題方面具有強(qiáng)大的優(yōu)勢(shì)，可以較好地用于潤(rùn)滑油振蕩問(wèn)題的可視化模擬研究，文中建立了發(fā)動(dòng)機(jī)曲軸-油底殼-潤(rùn)滑油流固耦合模型，采用CEL算法分析了潤(rùn)滑油在不同工況下的瞬態(tài)振蕩情況，在加速減速制動(dòng)等工況下，油底殼井槽后部區(qū)域受到了潤(rùn)滑油的不斷的晃動(dòng)和沖擊，應(yīng)力曲線也出現(xiàn)了多處尖峰；在轉(zhuǎn)向工況下，油底殼側(cè)面區(qū)域也受到了潤(rùn)滑油的不斷的晃動(dòng)和沖擊，應(yīng)力曲線也出現(xiàn)了多處尖峰，這對(duì)于發(fā)動(dòng)機(jī)的NVH性能有很大影響；此外，CEL方法對(duì)于后續(xù)噪聲的研究提供了新的思路，也為流固耦合振蕩問(wèn)題的數(shù)值分析探索了一個(gè)新的研究方法。

[1] 袁兆成. 內(nèi)燃機(jī)設(shè)計(jì)[M].北京：機(jī)械工業(yè)出版社，2008.

[2] 陳佐一.流體激振[M].北京：清華大學(xué)出版社，1998.

[3] 徐文杰. 基于CEL算法的滑坡涌浪研究[J]. 工程地質(zhì)學(xué)報(bào),2012,20(3):350-354.

XU Wenjie. CEL algorithm study of reservoir surge induced by land slide [J]. Journal of Engineering Geology, 2012, 20(3):350-354.

[4] 魏鵬,史勇杰,徐國(guó)華. 復(fù)雜旋翼流場(chǎng)的耦合歐拉-拉格朗日數(shù)值方法[J]. 航空學(xué)報(bào),2013,34(7):1538-1547.

WEI Peng, SHI Yongjie, XU Guohua. Coupled Eulerian-Lagrangian method for complicated rotor flow field prediction [J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(7): 1538-1547.

[5] NOH W F. CEL:A time-dependent two-space-dimensional coupled Eulerian-lagrangian code[C]∥Methods in Computational Physics, Volume 3,Fundamental Methods in Hydrodynamics. New York,NY: Academic Press, 1964:117-179.

[6] 王懿,賈旭,黃俊,等. 基于CEL的船舶拋錨入泥深度分析[J]. 石油機(jī)械,2014,42(12):44-47.

WANG Yi,JIA Xu,HUANG Jun,et al. Analysis of penetration depth of dropped anchor based on CEL [J].Chian Petroleum Machinery,2014,42(12):44-47.

[7] FOUCARD L, ARYAL A, DUDDU R, et al. A coupled Eulerian-Lagrangian extended finite element formulation for simulating large deformations in hyperelastic media with moving free boundaries[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 283:280-302.

[8] ROSTAMI J, ABBASSI A. Conjugate heat transfer in a wavy microchannel using nanofluid by two-phase Eulerian-Lagrangian method[J]. Advanced Powder Technology, 2016, 27(1): 9-18.

[9] DIGGS A, BALACHANDAR S. Evaluation of methods for calculating volume fraction in Eulerian-Lagrangian multiphase flow simulations[J]. Journal of Computational Physics, 2016, 313: 775-798.

[10] BAHREMAND H, ABBASSI A, SAFFAR-AVVAL M. Experimental and numerical investigation of turbulent nanofluid flow in helically coiled tubes under constant wall heat flux using Eulerian-Lagrangian approach[J]. Powder Technology, 2015, 269: 93-100.

[11] VUJANOVI M, PETRANOVI Z, EDELBAUER W, et al. Modelling spray and combustion processes in diesel engine by using the coupled Eulerian-Eulerian and Eulerian-Lagrangian method[J]. Energy Conversion and Management, 2016,125:15-25.

[12] JARAUTA A, RYZHAKOV P, SECANELL M, et al. Numerical study of droplet dynamics in a polymer electrolyte fuel cell gas channel using an embedded Eulerian-Lagrangian approach[J]. Journal of Power Sources, 2016, 323: 201-212.

[13] 葉林征，祝錫晶，王建青，等.基于 CEL 不同角度超聲空化微射流沖擊的仿真分析[J].振動(dòng)與沖擊，2016,35(16):130-134.

YE Linzheng, ZHU Xijing, WANG Jianqing,et al.Simulations of ultrasonic cavitation micro-jet impact with different angles based on CEL [J]. Journal of Vibration and Shock, 2016,35(16):130-134.

[14] 卿啟湘,楊會(huì),劉杰,等. 基于CEL法的EPB切刀切削土體仿真分析[J]. 工程設(shè)計(jì)學(xué)報(bào),2015,22(3):230-235.

QING Qixiang, YANG Hui, LIU Jie, et al. Simulation analysis on cutting soil for cutter of EPB based on CEL method [J]. Chinese Journal of Engineering Design, 2015, 22(3): 230-235.

[15] 姚小虎,黃愉太,歐智成,等. 基于CEL算法的水陸兩棲飛機(jī)水上降落動(dòng)力特性分析[J]. 華南理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,43(6):110-115.

YAO Xiaohu, HUANG Yutai, OU Zhicheng, et al.CEL algorithm-based analysis of dynamics characteristics of amphibious aircraft landing on water [J]. Journal of South China University of Technology (Natural Science Edition), 2015, 43(6): 110-115.

[16] 王建華,蘭斐. 鉆井船插樁對(duì)鄰近樁影響的耦合歐拉-拉格朗日有限元方法研究[J]. 巖土力學(xué),2016,37(4):1127-1136.

WANG Jianhua, LAN Fei. A coupled Eulerian-Lagrangian FEMmethod for analyzing the effects of spudcan penetration on an adjacent pile [J]. Rock and Soil Mechanics, 2016, 37(4):1127-1136.