陳曉燕，張小霞，張新華

（1.遼寧省清河水庫管理局，遼寧鐵嶺 112003；2.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川成都 610065）

1 研究背景

隨著計算機科學(xué)的不斷進步，很多理論上解難以解決的問題，都可以通過數(shù)值模擬的方式得到滿足工程實際需要的近似解。常用的數(shù)學(xué)模型多采用C語言或Fortran語言編寫的，在解決工程問題時，通常需要針對具體問題建立對應(yīng)的幾何或地形模型并劃分網(wǎng)格，若能用數(shù)學(xué)函數(shù)準確地反映實際地形將會極大地提高模擬結(jié)果的可靠性，也能更好地解釋一些工程現(xiàn)象及其產(chǎn)生的原因。

但是，實際可獲取的地形數(shù)據(jù)總是有限的，而且不一定是規(guī)整的、結(jié)構(gòu)化的。在利用C語言或Fortran語言編制數(shù)值模擬模型中添加精細化的地形模塊比較耗時，不能同步可視化并適時調(diào)整模型。根據(jù)已有數(shù)據(jù)及其具體地形特征，選取適宜的數(shù)值分析方法，并應(yīng)用MATLAB實現(xiàn)數(shù)據(jù)插值或數(shù)據(jù)擬合，求解插值函數(shù)則更為方便、快捷，并且易于核查分析結(jié)果的準確性。

2 理論基礎(chǔ)

2.1 數(shù)值分析方法

應(yīng)用有限的地形數(shù)據(jù)獲取較為精確的物理模型屬于數(shù)據(jù)擬合問題，在數(shù)值分析中，常被稱為數(shù)值逼近，是為離散數(shù)據(jù)（不論何種途徑得到）建立簡單連續(xù)模型的方法。數(shù)值逼近方法有插值法和擬合法兩種[1]，兩種方法又可根據(jù)處理方式或選取的函數(shù)類型不同分為許多具體方法，詳細方法可參考有關(guān)文獻［2～4］。

2.2 MATLAB簡介

MATLAB（Matrix Laboratory）是矩陣實驗室的簡稱，和Mathematica,Maple并稱為三大數(shù)學(xué)軟件。MATLAB主要應(yīng)用于工程計算、控制設(shè)計、信號處理與通訊、圖像處理、信號檢測、金融建模設(shè)計與分析等領(lǐng)域，具有六大主要特性和優(yōu)點：友好的工作平臺和編程環(huán)境[10]，簡單易用的程序語言，強大的科學(xué)計算機數(shù)據(jù)處理能力，出色的圖形處理功能，應(yīng)用廣泛的模塊集合工具箱，實用的程序接口和發(fā)布平臺。

目前，MATLAB已經(jīng)把工具箱延伸到了幾乎所有科學(xué)研究和工程應(yīng)用領(lǐng)域。

3 案列分析

3.1 研究工況

張小霞[12]在對自編二維的淺水模型(SWE)驗證時，采用了蔡榮[13]的丁壩定床水槽實驗數(shù)據(jù)。試驗水槽長7.2 m、寬2 m、高約30 cm，水槽底部鋪設(shè)平整，模型河床的縱比降為0.16%。實驗中丁壩布置在水槽中部，結(jié)構(gòu)如圖1所示。丁壩壩體為不透水丁壩，丁壩橫斷面采用梯形斷面，迎水壩坡為1.0∶1.5，背水壩坡為1∶2，壩頭坡面為光滑曲面，壩頭坡1.0∶3.5。流量為15 L/s，尾水水深3.36 cm。

圖1 實驗丁壩結(jié)構(gòu)示意圖（單位：cm）

3.2 數(shù)據(jù)分析及模型建立

為驗證SWE模型中水動力模塊的可靠性，將丁壩結(jié)構(gòu)作進一步分析。以丁壩左上角底部為坐標原點，建立坐標系如圖2所示，其中分別以壩寬、壩長、壩高方向為x，y，z軸。根據(jù)丁壩結(jié)構(gòu)示意圖可以將丁壩分為4個面，其中面1-3為平面，所有參數(shù)可以準確設(shè)置；面4為曲面，曲面4上各點的坐標確定需要進一步分析，由其形狀可以推知曲面4為一個二次曲面，可由Z=f（x，y）表示出來，式中x和y的最高次項都不超過2次。曲面4上有5條特征線分別為圖2中L1-L5，同時有5個邊界點，可以直接由丁壩結(jié)構(gòu)示意圖推求得到，為提高數(shù)值逼近精度又分別增加特征線L1-L4的4個中點，一共構(gòu)成9個特征點，根據(jù)所建坐標模型推求出各點坐標并列于表1。

圖2 丁壩特征點/線/面示意圖

3.3 Matlab模型求解

以上分析表明：需要求解的問題——曲面4的表達式，已知條件——曲面4上的幾個特征邊界點線及曲面4大致符合的函數(shù)類型，所求函數(shù)在已知數(shù)據(jù)空間內(nèi)，并且知道函數(shù)為二次曲面函數(shù)，結(jié)合插值法和擬合法兩種數(shù)值逼近方法的優(yōu)缺點及適用性，可知該問題更適合用擬合法求解。

表1 特征點坐標m

因此，采用MATLAB曲線擬合工具cftool(Curve Fitting Tool)實現(xiàn)對1—9點的二次曲面擬合。通過選用不同的擬合增強方式并調(diào)整x,y的最高項次數(shù)（2次或1次），得到的最佳擬合結(jié)果如圖3所示。

圖3 點1-9二次曲面擬合結(jié)果

此次擬合回歸平方和SSE=1.45×10-4，決定系數(shù)R2=0.990 3，擬合得到的曲面4的表達式：

比較圖3和丁壩示意圖圖1，可以發(fā)現(xiàn)模擬結(jié)果形式符合，但存在一定偏差，尤其是壩頭底部的曲線明顯偏向右側(cè)，和原實驗中丁壩壩頭底部曲線不符。為調(diào)整擬合圖形，可在丁壩壩頭底部曲線L5上增加2個控制點。為此需要先求出L5的表達式，通過點3，4，5三個點進行二次曲線擬合，最佳擬合結(jié)果如圖4所示，擬合曲線L5的表達式為式（2）。

圖4 丁壩壩頭底部曲線擬合結(jié)果

由式（2），推求出曲線L5上兩點10（0.08，0.79，0.00）和11（0.25，0.83，0.00），加上前面表2中列出的9個點，用此11個點對曲面4進行二次擬合，調(diào)整不同的擬合增強方式并改變x,y的最高項次數(shù)，得到的最佳擬合結(jié)果如圖5所示，其回歸平方和SSE=1.399×10-4，決定系數(shù)R2=0.992 7，擬合得到的曲面4的表達式為式（3）。比較圖5和丁壩示意圖圖1，發(fā)現(xiàn)模擬結(jié)果符合較好。

圖5 點1-11二次曲面擬合結(jié)果

3.4 擬合結(jié)果應(yīng)用情況

在SWE模型中采用擬合式3求取壩頭各網(wǎng)格點上的高程，得到初始河道地形圖如圖6所示。對丁壩細節(jié)應(yīng)用數(shù)值擬合，并采用MATLAB可視化分析得到了較為準確細致的丁壩結(jié)構(gòu)，將圖6與圖1相比較，上述方法具有較高的真實性，應(yīng)用模擬函數(shù)可以構(gòu)建出與上下游壩面以及河底均能平滑相接的地形模型。

圖6 SWE模型水動力驗證模擬試驗地形圖

為檢驗丁壩模型對水動力模擬結(jié)果的影響，分別應(yīng)用上述擬合丁壩結(jié)構(gòu)和簡化為立方體結(jié)構(gòu)的丁壩進行河槽水動力計算。除丁壩結(jié)構(gòu)外，其余參數(shù)均按照蔡榮實驗進行設(shè)置。計算區(qū)域長7.2 m、寬2 m，采用非均勻網(wǎng)格，在丁壩位置處網(wǎng)格較密，離丁壩越遠網(wǎng)格尺寸越大。對于擬合丁壩X，Y方向的最小網(wǎng)格均為0.01 m，總網(wǎng)格數(shù)為342×82個，時間步長取為0.001 s/步。對于簡化丁壩采用長0.67 m、寬0.35 m、高0.1 m的立方體，在計算區(qū)域內(nèi)X，Y方向的最小網(wǎng)格均為0.02 m，總網(wǎng)格數(shù)為222×62個，時間步長取為0.01 s/步。

圖7為模擬流速與蔡榮實驗測得的流速分布圖，其中圖7a為蔡榮實測數(shù)據(jù)。對擬合丁壩水動力模擬和簡化丁壩水動力模擬結(jié)果提取相關(guān)斷面的流速數(shù)據(jù)，應(yīng)用tecplot繪制流速矢量分別如圖7b和圖7c所示。整體上SWE模型計算結(jié)果都很好地模擬出了丁壩作用下水流流速分布，與實測保持較高的一致性。在丁壩附近，擬合丁壩的計算流速與實驗更加接近。

模擬結(jié)果與蔡榮實驗數(shù)據(jù)誤差分析如表2所示。如圖7（a）所示，測量了9個斷面的流速，從左到右依次編號。在有測量數(shù)據(jù)的點（速度為零的點不包括在內(nèi)）位上，分別比較擬合丁壩和簡化丁壩計算結(jié)果的相對誤差和絕對誤差，然后取平均求得斷面平均的相對誤差和絕對誤差分別列于表2中。

圖7 應(yīng)用結(jié)果比較圖

表2 特征點坐標

簡化丁壩整體平均流速絕對誤差為0.056 m/s，比擬合丁壩整體平均流速絕對誤差（0.047 m/s）高出近0.01 m/s；相應(yīng)的整體平均相對誤差高出3%。由此可見，應(yīng)用上述擬合所得地形模型進行數(shù)值模擬計算，得到水動力模擬結(jié)果與蔡榮實驗數(shù)據(jù)吻合更好，可降低或避免由于丁壩設(shè)置不妥造成的模擬誤差。

4 結(jié)語

針對大型數(shù)值計算中常遇到的物理模型二次曲面擬合問題，以丁壩壩頭曲面為例，應(yīng)用MATLAB中的數(shù)值分析工具成功的擬合了丁壩壩頭曲面，提高了數(shù)值計算的精度和效果。當(dāng)然它也存在一些不利影響，為了保證模擬計算穩(wěn)定，提高物理模型精度的同時必然要求更小的網(wǎng)格尺寸及更短的時間步長，在以上案列中擬合丁壩模擬耗時遠遠高于簡化丁壩模擬耗時。另外，為了取得較好的擬合精度，應(yīng)特別注意，獲取盡量準確的原始數(shù)據(jù)以及體現(xiàn)模型特征的數(shù)據(jù)點。

［1］張文娟，張國華，楊繼業(yè).基于廣義延拓的數(shù)值逼近方法［J］.紡織高校基礎(chǔ)科學(xué)學(xué)報，2007，20（3）：281—285.

［2］奚梅成.數(shù)值分析方法［M］.合肥：中國科學(xué)技術(shù)大學(xué)出版社，2004.

［3］南京大學(xué)數(shù)學(xué)系.數(shù)值逼近方法［M］.北京：科學(xué)出版社，1978.

［4］程正興.數(shù)據(jù)擬合［M］.西安：西安交通大學(xué)出版社，1986.

［5］吳亞娟，劉曉鋒，趙興強.MATLAB在《計算方法》教學(xué)與實驗中的應(yīng)用［J］.通化師范學(xué)院學(xué)報，2006，27（2）：33—36.

［6］韓艷麗，李鳳萍.插值與曲線擬合［J］.中國西部科技，2009，8（11）：91—92.

［7］李慶揚，王能超，易大義.數(shù)值分析［M］.武漢：華中理工大學(xué)出版社，1982.

［8］邱斌，樂科軍，朱建軍，等.衛(wèi)星精密星歷的曲線移位插值方法［J］.工程勘察，2008（11）：64—68.

［9］張雪峰，李文林.Matlab在高等數(shù)學(xué)中的幾點應(yīng)用［J］.安陽師范學(xué)院學(xué)報，2009，2009（5）：121—123.

［10］姬慧玲，孔群璐，馬紅星.基于matlab的電力系統(tǒng)潮流計算方法：CN，CN 103441496 A［P］.2013.

［11］吳云天，史啟朝.三維模擬四維的Matlab實現(xiàn)［J］.陜西科技大學(xué)學(xué)報，2009，27（4）：146—150.

［12］張小霞.丁壩對改善水生生物棲息地作用的數(shù)值模擬研究［D］.成都：四川大學(xué)，2016.

［13］蔡榮.瀾滄江丁壩局部沖刷試驗及其防護研究［D］.重慶：重慶交通大學(xué)，2012.