李 濤 劉 浩 何 綱
1.蘇州科技學院,蘇州,215009 2.南京航空航天大學,南京,2100163.河海大學,常州,213022
保持參數分布狀態的有限元變形網格B樣條曲面重建算法
李濤1劉浩2何綱3
1.蘇州科技學院,蘇州,2150092.南京航空航天大學,南京,2100163.河海大學,常州,213022
針對變形前后的有限元網格模型及原始曲面模型,提出一種保持參數分布狀態的B樣條曲面重建算法。在裁剪區域,通過對原始曲面進行平移、延伸、截取及重新參數化等操作,構造初始擬合曲面;在非裁剪區域,在垂直于大曲率邊界的參數方向上插入截面線并重新參數化兩條大曲率邊界,用雙向蒙皮的方法重建擬合曲面。然后進入重新參數化網格點和重新擬合曲面的迭代過程,直至滿足終止準則。在曲面迭代修改過程中,通過計算基函數的極值點給出一種更精確的欠約束區域判定方法;利用插值四邊/三角細分算法細化粗糙網格模型,補充約束條件,提高擬合曲面的光順性;借助擬合曲面補充邊界約束條件,結合網格點插值和形狀保持約束,改進裁剪曲面的擬合精度。實驗結果驗證了算法的有效性。
參數分布;欠約束區域;雙向蒙皮;插值四邊/三角細分
在板料成形計算機輔助工程(CAE)分析過程中,通常要把計算機輔助設計(CAD)得到的曲面模型劃分成有限元網格單元,然后用CAE軟件模擬沖壓過程進行成形性評估,再修改曲面模型并返回至CAE系統繼續分析[1]。這是一個循環反復的過程。由于經過分析以后的網格模型往往發生變形,原來的曲面模型不能直接輸入到CAD/CAM系統中,必須把修改后的網格轉化成曲面模型,因此,設計高效、穩定的曲面重建算法對縮短產品設計周期、節省開發成本具有極其重要的作用。
曲面重構是逆向工程、計算機圖形學和CAD等領域的一個經典問題,目前已有大量的研究成果[2-8]。按照待擬合數據是否規則,傳統的自由曲面重構算法大致可以分為兩大類:一類是針對有序的截面數據,利用蒙皮的方法構造單張或多張曲面[2-3];一類是針對大規模的散亂數據,通過擬合的方法構造滿足一定精度和光順性要求的逼近曲面[4-6]。由于對原始幾何信息未知,在曲面重建時不可能考慮參數分布,而不良的參數分布可能導致曲面在隨后的網格剖分時產生大量的非結構三角網格,不利于有限元分析。另一方面,有限元網格不同于逆向工程中的測量數據,網格通常劃分得比較稀疏,且自適應采樣的過程也往往使網格點分布不夠均勻。因此,完全借助逆向工程軟件或按照傳統的曲面擬合方法對有限元網格進行曲面反求,效果并不理想。由于有限元網格模型是由CAD曲面模型轉化而來的,變形網格對應的原始幾何信息是已知的,且這些信息正是需要計算出來進行分析和比較的,在曲面重構時應該予以考慮。
關于有限元網格模型的曲面重建,文獻[6]假設已知變形前后網格點的對應關系,采用類似逆向工程里的方法直接進行B樣條曲面的重建,沒有考慮網格分布不均時的異常處理和擬合曲面的參數分布;文獻[7]把對原始曲面的迭代修改和對變形網格的曲面重構結合在一起,但是曲面修改的方法難以有效地控制邊界誤差;文獻[8]把雙向蒙皮、曲面掃掠、延伸等正向設計的思想用于有限元網格模型的曲面重建,取得了較好的效果,但是仍有一些問題沒有解決:①重建裁剪曲面時沒有考慮原始曲面參數分布;②對于非裁剪區域,雙向蒙皮的方法可以構造出逼近程度較高的初始擬合曲面,但是當邊界曲率變化較大,特別是當邊界呈“U”形或“S”形時,蒙皮曲面的參數分布不夠理想;③對擬合曲面進行迭代修改時,按照文獻[4]的方法在被裁剪區域或稀疏網格點處補充形狀保持約束條件,雖然可以保證約束方程組的可解性,但是當大量欠約束區域存在時,該方法構造的擬合曲面光順性較差。
本文主要解決上述幾個關鍵問題,根據輸入的原始曲面和網格及變形網格模型,重建新的曲面模型。本文假設變形前后的網格并無確定的對應關系,但是變形網格已經劃分為若干個與原始曲面模型相對應的區域。
1.1非裁剪域的判定

1.2截面線的插入
在一些鈑金件模型中,通常會含有一些狹長的區域,這些區域大都含有四條比較規則的邊界線,可以直接用雙向蒙皮方法構造插值四條邊界的初始擬合曲面[8],但是由于缺乏對曲面內部形狀的約束,當邊界曲率比較大時,內部網格點與曲面的偏差會比較大,曲面內部呈現拱起或凹陷的現象,如圖1所示算例,圖1b是對圖1a中的網格直接插值四條邊界得到的蒙皮曲面,內部網格點的最大誤差(如無特別說明,本文所說的誤差均為相對誤差,即絕對誤差與曲面片最小軸向矩形包圍盒對角線長度的比值)為3.7310-2,平均誤差為1.8810-2。另一方面,當曲面邊界曲率較大時,由于網格點分布不均勻(如圖1a中曲率較大的兩條對邊),蒙皮曲面的參數分布往往不夠理想,如圖1d所示。解決上述2個問題的方法是在垂直于大曲率邊界的參數方向上插入一條內部截面線,這樣一方面可以在蒙皮時增加內部截面線約束,另一方面,可以利用此截面線調整與其相交的兩條對邊的參數分布,使之協調一致。

(a)網格模型

(b)蒙皮曲面光照圖

(c)原始曲面的參數分布

(d)曲面與原曲面參數對比圖1 直接插值四條邊界的雙向蒙皮曲面
在兩條大曲率對邊上各選取離中點最近的一點作為截面線候選起始點,然后比較兩條曲線在相應點處的曲率值。為了提高算法的穩定性,取曲率較小的候選點作為新截面線的起始端點,如圖1a中的A點。以邊界線在該點的切線方向為參考方向d1,在變形網格上沿網格邊搜索一條大致垂直于該方向的路徑,直至對邊邊界,如圖1a中的路徑AB。由于有限元網格具有自適應采樣的特點,一般而言,上述方法得到的路徑并不光滑,作為曲面蒙皮和重新參數化至關重要的一條中間截面線,需要進一步調整。方法如下:在這條新的截面線上尋找一條光滑的最長子路徑,如果這段子路徑的頂點數大于2,且其走向d2與前面的參考方向d1大致垂直,則把d2作為投影方向,向兩條對邊邊界重新投影并延伸這段光滑子路徑,得到修正后的截面線。如果網格模型的結構性很差,在初始截面線上找不到上述光滑子路徑,則直接從起始端點開始根據參考方向d1向對邊邊界投影。如果對邊邊界的曲率比較大,這樣直接投影的效果不一定理想,可以根據投影曲線與對邊邊界的夾角情況進行迭代修改,使新的截面線與兩條相交邊界的夾角之差最小。
1.3重新參數化邊界線
在構造插值四條邊界的雙向蒙皮曲面時,由于網格模型具有自適應采樣的特性,對于曲率較大的兩條對邊,傳統的參數化方法有時不能得到協調的參數分布,從而影響蒙皮曲面的參數分布,如圖1d所示。因此,需要重新參數化曲率較大的兩條對邊,以調整蒙皮曲面的參數分布。文獻[9]在計算降階前后兩條Bezier曲線的Hausdorff距離時,設置局部曲率極值點或參數中點為分段點,然后對逼近曲線的參數引入分段線性變換,從而實現曲線的重新參數化。本文借鑒該方法,利用新插入的中間截面線調整兩條大曲率對邊的參數分布。

(1)
經過修改后,中間截面線與該邊界交點參數由m1變為m,其他點的參數也隨之做線性調整,另一條曲線的調整方法類似。圖2是對圖1所示算例插入一條中間截面線并調整參數后得到的蒙皮曲面,內部網格點的最大誤差為9.49×10-3,平均誤差為4.09×10-3。比較可見,插入一條截面線后,蒙皮曲面無論在參數分布還是擬合精度上都有較大的改善。

(a)蒙皮曲面與原曲面參數對比

(b)蒙皮曲面光照圖圖2 插入一條中間截面線后構造的蒙皮曲面
初始擬合曲面的曲面品質和逼近精度對于最終的擬合結果都有至關重要的影響。作為參數化投影基曲面,采用B樣條擬合曲面比最小二乘擬合平面可以得到更好的參數化結果[3]。
關于裁剪曲面的擬合,由于被裁剪區域幾何信息的缺失,初始基曲面的構造是一個難題。因為裁剪曲面的外環邊界形狀不定,難以直接根據裁剪邊界確定基曲面的參數走向。因此,用曲面重建的方法構造的擬合曲面與原始曲面的參數分布往往有很大的差異[8]。文獻[4]用最小二乘擬合平面作為投影基曲面參數化擬合數據。這種方法不適合大曲率曲面,因為投影參數化方法可能導致參數分布很不理想。
本文通過平移并修改原始基曲面的方法構造與之有相同參數分布的裁剪曲面。對每一張曲面片,根據變形前后網格模型估算一個大致的偏置矢量,把原始曲面加上偏置矢量后得到的新曲面作為裁剪曲面的初始近似。以新曲面為投影面計算變形網格點的參數,同時得到網格點到其投影點的距離,即擬合誤差。如果所有網格點的誤差滿足給定閾值,即變形曲面剛好是原曲面的一個近似平移,則平移后的新曲面即可作為最終的擬合曲面;否則,檢查是否有邊界點沒有投影到基曲面上,若出現上述情況,則根據最大偏離點的位置延伸基曲面,從而使所有網格點都投影到基曲面上。
按照上述方法構造的初始擬合曲面,特別是裁剪曲面的初始基曲面,往往不能滿足用戶的誤差要求,需要進一步修改,以降低擬合誤差。主要思想是用最新得到的擬合曲面作為參數化投影曲面,重新計算網格點參數和誤差,如果達到誤差需求或最大迭代次數,則停止計算;否則重新擬合。下面著重介紹曲面迭代修改過程中的一些關鍵技術細節。3.1基曲面的預處理
原始曲面模型本身可能存在一些設計缺陷,如裁剪曲面的基曲面過大,參數分布不佳等。直接以這種曲面的偏置面為初始曲面進行迭代修改,往往難以保證擬合曲面的精度和光順性,因此,對這種有缺陷的基曲面需要加以修改。
若裁剪曲面的基曲面過大,如圖3a、圖3b所示,基曲面的大部分區域屬于“無用”區域被裁剪掉,且這部分區域在曲面重新擬合過程中沒有任何作用,同時保留區域的約束作用被弱化,故采用曲面截取的方式[10]去掉部分“無用”區域。圖3c是截取后的裁剪曲面及其基曲面。

(a)原始裁剪曲面

(b)原始裁剪曲面及其基曲面

(c)截取后的裁剪曲面及其基曲面

(d)重新參數化后的裁剪曲面及其基曲面圖3 初始基曲的預處理示例
若基曲面的參數分布不佳,不僅降低擬合算法的效率,而且影響后續的應用,如有限元網格劃分,數控刀軌生成等。引起曲面參數分布不理想的原因主要有兩個:一是B樣條基函數的節點布局不均勻,二是曲面控制頂點的空間分布不均勻。對于前者,通過插入或刪除節點[10]方法使其分布盡量均勻,即計算基函數相鄰互異節點的節點距ki和節點距的平均值ka,若ki>rka,則在兩節點之間插入[ki/(rka)] - 1個節點,其中[·]表示取整操作,r為比例因子,通常可以取3~6之間的值;若ki 3.2約束方程組的建立 根據網格點參數和初始擬合曲面 (2) 可以建立約束線性方程組 AD=B (3) 式中,di,j為曲面控制頂點;Ni,k(u)與Nj,l(v)分別為k次和l次B樣條基函數;A為將網格點參數代入曲面基函數得到的系數矩陣;D為由曲面控制頂點di,j構成的未知量;B為由網格點坐標構成的常矢量。 由于有限元網格點分布的不均勻性以及大量裁剪曲面的存在,方程組可能出現欠約束與過約束同時存在的問題,即系數矩陣A可能行或列不滿秩,或者兩者都不滿秩。這給線性系統(式(3))的求解帶來了困難。為了保證(式(3))解的唯一性與穩定性,本文采取的策略是通過尋找欠約束區域并為之適當補充約束條件,構造一個系數矩陣恒滿足列滿秩的新的線性系統: (4) 3.3欠約束區域的判定 文獻[4]用均勻B樣條基函數節點區間的中間部分作為有效域進行孔洞識別,如果落在有效域內的網格點參數個數少于一定的閾值,則認為該區域為欠約束區域,需要在此補充形狀保持約束條件: -0.25(di-1,j-1+di-1,j+1+di+1,j-1+di+1,j+1) + 0.5(di-1,j+di,j-1+di,j+1+di+1,j)=di,j (5) 從而提高了最小二乘解的穩定性。但是上述有效域的判別方法并不適應于一般的非均勻B樣條曲面,因為當基函數節點分布不均,特別是當基函數含有內部重節點時,最大值點可能出現在節點附近。因此,在此給出一種更加準確的欠約束區域識別方法。 圖4 有效域示意圖 3.4混合網格的插值細分 在逆向工程中,通過對欠約束區域施加頂點形狀保持約束可以提高最小二乘解的穩定性[4]。然而,對于比較稀疏的有限元網格,由于大量欠約束區域的存在,該方法的擬合效果并不理想。在圖5所示的算例中,圖5c是對圖5a中的網格用文獻[4]的方法擬合得到的結果,雖然最大誤差(3.23×10-3)達到一般的工程要求,但是擬合曲面的光順性較差。導致擬合曲面品質較低的主要原因是網格模型過于粗糙,插值網格點位置的“剛性”約束條件不足。解決上述問題的一個有效辦法是增加網格點密度,補充“剛性”約束條件。顯然,線性插補會降低網格模型的光順性,影響擬合曲面的品質,因此考慮網格模型的插值細分。 (a)待擬合的網格模 (b)細分二次后的網格圖 (c)擬合原網格所得曲面 (d)擬合細分網格所得曲面圖5 通過對網格插值細分增加約束條件 由于有限元網格大都是由四邊形和三角形構成的混合網格,本文采用四邊/三角混合網格的插值細分。目前關于混合網格插值細分算法的研究比較多[11-13],針對有邊界的不規則混合網格插值細分,本文采用文獻[13]的方法,把細分算法分解為線性分割和幾何平均兩個過程,根據插值細分與相應逼近型細分的關系,把線性分割得到的新邊點和新面點加上相應偏置量,得到最終的新位置,其中新邊點和新面點的偏置量是以網格頂點相對于相應逼近型細分點的偏置量為基礎計算得到的,權因子為相應的逼近型細分模板。圖5b是對圖5a所示網格細分二次所得到的網格圖,比較圖5a和圖5b可見,細分以后網格點更加稠密,網格模型也更加光順,圖5d是用細分后的網格擬合得到的曲面模型,最大誤差為2.73×10-3,與圖5c相比,光順性有了顯著的提高。 網格細分方法在增加“剛性”約束條件的同時也加大了運算開銷,因此,網格細分與否是以網格點是否過于稀疏或者是否存在大量欠約束區域為判定準則的。 3.5狹長裁剪面的迭代修改 (a)網格模型 (b)擬合曲面圖 (c)圖6b的局部放大 如1.2節所述,鈑金件幾何模型中的狹長區域,大部分可以用雙向蒙皮的方法重構出高精度和高品質的擬合曲面,但是也有部分狹長區域沒有形狀規則的四條邊界,如圖6a所示的網格模型,這時只能作為裁剪曲面來處理。由于基曲面的大部分區域被裁剪掉,因此曲面迭代修改時會出現幾何約束不足的問題。雖然文獻[4]根據文獻[14]給出了形狀保持約束,但是當存在大量欠約束區域,尤其是基曲面的邊界無法得到有效控制時,擬合曲面的光順性難以得到保證,原因是Farin的持久準則[14]實際上是一個離散形式的Euler-LagrangePDE,邊界約束必不可少。本文給出的解決方案是,將前一次擬合曲面的邊界信息和網格點一樣,作為曲面修改的“硬約束”,其他內部欠約束區域作為“軟約束”由形狀保持準則給出。作為“硬約束”,初始擬合曲面的邊界誤差不能太大,如果誤差過大,如絕對誤差大于0.5mm或相對誤差大于0.05,則將該曲面先做變形處理:根據網格點的誤差計算曲面的法向偏置量,如果曲面采樣點參數在裁剪域內,則其偏置量由鄰近網格點的偏置量內插得到;否則,由已知偏置量的鄰近點外插得到,然后根據偏置信息重新擬合曲面。這樣得到的變形曲面比原來的擬合曲面更加“貼近”網格模型,可以用來補充邊界插值約束條件。把邊界點插值約束和內部形狀保持約束(式(5))加入線性系統(式(3))中,得到最終的約束方程式(4)。在圖6所示的算例中,圖6b是用上述方法對圖6a擬合的結果,最大誤差為7.60×10-4,平均誤差為1.16×10-4,最大絕對誤差為0.436mm。筆者同時用文獻[7]和文獻[4]的方法對此算例做了實驗,文獻[7]在去掉光順項約束后,最大絕對誤差仍為0.808mm,而文獻[4]的方法由于擬合曲面嚴重扭曲而導致算法失敗。 (d)擬合裁剪曲面及其基曲面圖6 狹長裁剪面示例 圖7是一個綜合算例,其中圖7a和圖7b分別為原始曲面和擬合曲面的線框圖,圖7c是擬合曲面的光照圖。圖8是另一模型的擬合結果。兩個算例的運行情況見表1,其中最大絕對/相對誤差是整個模型所有擬合曲面的絕對/相對誤差的最大值,平均絕對誤差是整個模型所有擬合曲面絕對誤差的平均值。本文所有算法都用C++語言編程實現,程序運行的軟硬件環境為WindowsXP操作系統,VC++ 6.0編譯器和主頻2.0GHz,內存2.0GB的筆記本電腦。由算例可知,本文算法具有較高的精度和較快的速度,擬合曲面不僅具有較好的光順性,而且具有與原始曲面基本一致甚至更佳的參數分布狀態。 (a)原始曲面線框圖 (b)擬合曲面線框圖 (c)擬合曲面光照圖圖7 綜合算例1 (a)原始曲面線框圖 (b)擬合曲面線框圖圖8 綜合算例2 算例曲面片數網格單元數網格點數最大絕對誤差/最大相對誤差平均絕對誤差運行時間(s)算例119015933158910.245mm/2.94′10-22.18×10-2mm55.4算例240623475217380.460mm/1.83′10-22.40×10-2mm52.1 基于有限元網格與逆向工程中的大規模網格模型的不同,提出了更有針對性的有限元變形網格曲面重建算法。曲面蒙皮、延伸、截取、重新參數化等許多正向設計的方法被用于初始擬合曲面的構造;在曲面的迭代修改過程中,通過計算基函數的極值點給出穩定性更好的有效域判別方法;通過混合網格插值細分算法細化部分稀疏的網格,補充必要的插值約束條件,提高了擬合曲面的光順性;利用最新擬合的基曲面或其變形曲面,補充部分狹長裁剪曲面的邊界插值約束,給出了穩定性更高的裁剪曲面擬合算法。本文方法構造的擬合曲面不僅光順性好,誤差小,運算效率也很高,達到了實時的要求。更為重要的是,擬合曲面與原始曲面具有基本一致甚至更優的參數分布狀態,可以幾乎無需修改而直接用于后續的CAE或CAM。當然,部分狹長區域裁剪曲面的擬合誤差雖然較現有方法有所改進,但有時仍然偏高,這是算法有待進一步提高的地方。 [1]陳文亮. 板料成型CAE分析教程[M], 北京:機械工業出版社, 2005. [2]ParkH.B-splineSurfaceFittingBasedonAdaptiveKnotPlacementUsingDominantColumns[J].Computer-AidedDesign, 2011,43(3): 258-264. [3]MaWeiyin,KruthJ.ParameterizationofRandomlyMeasuredPointsforLeastSquaresFittingofB-splineCurvesandSurfaces[J].Computer-AidedDesign, 1995,27(9):663-675. [4]譚昌柏, 周來水, 張麗艷,等. 裁剪B樣條曲面重建算法研究[J]. 中國機械工程, 2007, 18(19): 2366-2370. TanChangbai,ZhouLaishui,ZhangLiyan,etal.ResearchonTrimmedB-splineEurfaceReconstruction[J].ChinaMechanicalEngineering, 2007, 18(19): 2366-2370. [5]HuangYunbo,QianXiaoping.DynamicB-splineSurfaceReconstruction:ClosingtheSensing-and-modelingLoopin3DDigitization[J].Computer-AidedDesign, 2007,39(11): 987-1002. [6]譚光華, 陳志楊, 張三元,等. 基于原始曲面信息的變形網格的曲面重構[J]. 中國機械工程, 2009, 20(1): 38-43. TanGuanghua,ChenZhiyang,ZhangSanyuan,etal.SurfaceReconstructionforDeformedMeshBasedonOriginalSurfaceInformation[J].ChinaMechanicalEngineering, 2009, 20(1): 38-43. [7]曾建江, 李衛國, 陳文亮. 基于有限元網格變形的飛機外形曲面修改[J]. 南京航空航天大學學報, 2008, 40(4): 534-538. ZengJianjiang,LiWeiguo,ChenWenliang.MappingAircraftSurfacesbyDeformationofFiniteElements[J].JournalofNanjingUniversityofAeronautics&Astronautics, 2008, 40(4): 534-538. [8]LiTao,WangYuguo,ChenWenliang.PiecewiseSmoothSurfacesReconstructionforFiniteElementMixedMeshes[C]//InternationalConferenceonDigitalManufacturing&Automation.Changsha,China, 2010: 90-94. [9]ChenXiaodiao,MaWeiyin,PaulJ.Multi-degreeReductionofBezierCurvesUsingReparameterization[J].Computer-AidedDesign, 2011,43(2):161-169. [10]PieglL,TillerW.TheNURBSBook[M].2Ed.Berlin:SpringerPress, 1996. [11]JiangQingtang,LiBaobin,ZhuWeiwei.InterpolatoryQuad/triangleSubdivisionSchemesforSurfaceDesign[J].ComputerAidedGeometricDesign, 2009, 26(8): 904-922. [12]LinShujing,LuoXiaonan,YouFang,et.al.DeducingInterpolatorySubdivisionSchemesfromApproximatingSubdivisionSchemes[J].ACMTransactiononGraphics, 2008, 27(5):article146(1-7). [13]LinShujing,LuoXiaonan,XuSonghua,etal.ANewInterpolationSubdivisionSchemeforTriangle/QuadMesh[J].GraphicalModels, 2013, 75(5): 247-254. [14]FarinG,HansfordD.DiscreteCoonsPatches[J].ComputerAidedGeometricDesign, 1999, 16(7): 691-700. (編輯王艷麗) Reconstruction of B-spline Surfaces from Finite Element Deformed Meshes with Consistent Parametric Distribution Li Tao1Liu Hao2He Gang3 1.Suzhou University of Science and Technology,Suzhou,Jiangsu,215009 2.Nanjing University of Aeronautics and Astronautics,Nanjing,210016 3.Hohai University,Changzhou,Jiangsu,213022 According to the information of original surfaces, original meshes and deformed meshes, a method for reconstructing B-spline surfaces with consistent parametric distribution was presented. In trimmed regions, initial fitting surfaces were constructed by some operations on original surfaces such as translation, extension, segment and reparametrization, etc. In untrimmed regions, bidirectional skinning method was adopted for the production of initial surfaces after inserting section lines across the highly curved boundaries and revising their parameterization. Then a loop of reparameterizing the nodes and refitting the surfaces proceeded until some ending rules were satisfied. During the process of iteration, a more precise method was presented for determining under-constrained regions by calculating extreme points of basis functions; fairness of the fitting surfaces was improved evidently by splitting coarse meshes to provide enough constraints with interpolatory quad/triangle subdivision scheme. Fitting precision is improved for trimmed surfaces by dint of boundaries and internal nodes interpolation constraints together with shape preservation constraints of under-constrained regions. Test results indicate effectivity of the proposed method. parametric distribution; under-constrained region; bidirectional skinning;interpolatory quad/triangle subdivision 2014-02-18 國家自然科學基金資助項目(51175248,51375141);江蘇省省屬高校自然科學研究資助項目(12KJB460009) TP391< class="emphasis_italic">DOI :10.3969/j.issn.1004-132X.2015.05.019 李濤,男,1978年生。蘇州科技學院數理學院講師。主要研究方向為CAGD。發表論文10余篇。劉浩,男,1972年生。南京航空航天大學機電學院副教授。何綱,男,1975年生。河海大學機電學院副教授。











4 算例與應用






5 結語