賈濤 牛華東 解景濤 盧瑞軍 蘇茂輝

摘要：CAE在汽車和發(fā)動機設計中發(fā)揮越來越重要的作用，然而CAE的本質(zhì)是求解各種常微分方程組或偏微分方程組;常微分方程又稱動力系統(tǒng)方程，通常用來求解動力系統(tǒng)問題。偏微分方程則廣泛應用于電磁學、流體力學、結(jié)構力學等多個領域;對于偏微分方程通常很難求得其解析解，需要借助數(shù)值計算來獲取其方程的近似解。拉普拉斯方程（Laplace）是一種橢圓形二階偏微分方程，并且可以求取其解析解。本文通過解析法以及數(shù)值法對拉普拉斯方程求解，并對比不同求解方法的效率和精度;結(jié)論顯示解析法雖然精度較高，但是需要很大的計算量，并且大多數(shù)偏微分方程沒有解析解。因此，在汽車和發(fā)動機等工程應用中應該根據(jù)精度需求選擇最優(yōu)的途徑求解偏微分方程問題。

關鍵詞：Laplace方程;偏微分方程;數(shù)值計算;CAE

1? 理論分析

拉普拉斯方程的數(shù)學表達式如下，

2? 問題描述

研究某矩形二維傳熱問題，已知傳熱問題的導熱微分方程式如下[2]：

3? 解析解

求解偏微分方程的核心在于邊界條件，不同的邊界條件甚至會出現(xiàn)完全不同的解;根據(jù)該問題描述，通過邊界條件將上述方程簡化：

由于公式（5）比較繁瑣，并且涉及到n（0～∞）求和計算;方程組規(guī)模小則計算量小，可以快速得到計算結(jié)果;方程組規(guī)模大則計算量很大，通常需要很久才能完成計算。本例中算法1取n=110，對5000個數(shù)據(jù)編程計算，則耗時近1小時;而采用算法2則只需要數(shù)秒時間。

算法1：

clc;

clear;

sum=zeros（100，50）;

syms x;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=int（100*sin（k*pi*x/50），0，50）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end

end;

算法2：

clc;

clear;

sum=zeros（100，50）;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=5000*（1-cos（k*pi））/（k*pi）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end;

end;

4? 數(shù)值解法

通過數(shù)值方法求解偏微分方程是常用的求解方法，如下通過有限元方法以及EXCEL迭代、軟件編程等不同方法進行數(shù)值計算，對比不同計算方法的精度和計算效率。

4.1 有限元方法

借助有限元軟件，建立平面模型并搭建穩(wěn)態(tài)熱力學模型，生成5151個節(jié)點、10000個單元;設置邊界條件如上所述，計算得到的溫度場分布如圖2所示：該方法更偏向于工程的應用。在求解過程中，模型搭建過程花費較多時間;但是求解計算效率較高，普通電腦計算只需要20sCPU時間。

4.2 Excel迭代方法

設置步長dx=dy=1，采用五點差分法求解;如圖4所示，使用EXCEL建立數(shù)值計算表格，啟用迭代計算，設置誤差以及最大迭代次數(shù)進行求解;該方法計算效率較高，只需要1～2s CPU時間;得到計算結(jié)果如圖5[3]。

4.3 軟件編程計算

借助軟件編寫五點差分程序，使用高斯賽德爾迭代法，同樣可以計算出最終的結(jié)果;編程耗時較長，計算耗時60s左右的CPU時間[4]。

計算程序：

clear

clc

u=1：1：50;

v=1：1：50;

[x，y]=meshgrid（u，v）;%通過坐標軸點在平面畫網(wǎng)格

p=[x，y];

plot（p）;

spy（p）;

for i=2：50

for j=2：100

p（i，j）=（j-1）+98*（i-2）;

end;

end;

for i=1;

for j=1：100

p（i，j）=0;

end;

end;

for i=50;

for j=1：100

p（i，j）=0;

end;

end;

for j=1;

for i=1：50

p（i，j）=0;

end;

end;

for j=100;

for i=1：50

p（i，j）=0;

end;

end;

PP=delsq（p）;

load（'my.dat'）;

b=spconvert（my）;

N=length（b）;? ?%Gauss – Seidel迭代法

u=inv（PP）*b;

u=zeros（N，1）;%迭代初始值

D=diag（diag（PP））;

E=-tril（PP，-1）;%下三角

F=-triu（PP，1）;%上三角

B=inv（D-E）*F;g=inv（D-E）*b;

eps=0.000001;

for k=1：10000 %最大迭代次數(shù)

y=B*u+g;

if norm（u-y）<eps

break;

end

u=y;

end

s=zeros（50，100）;

for i=1：48

for j=1：98;

s（i+1，j+1）=u（98*（i-1）+j，1）;

end

end

for i=1：50

s（i，1）=100;

end

subplot（1，1，1）;

mesh（s）;

5? 誤差分析

對比計算結(jié)果第二列（x=1，u（x，y））數(shù)據(jù)并繪制圖表;如圖8所示，數(shù)值解和解析解的四條曲線重合在一起，說明這些數(shù)值大小相同，整體趨勢一致。圖9誤差分析顯示，EXCEL和編程計算誤差重合在一起，誤差在0～0.01%之間;有限元計算在中間部分誤差接近0，說明其計算誤差相對較小;三種數(shù)值計算方法誤差在工程可以接受的范圍內(nèi)。

6? 總結(jié)

本文通過不同方法求解偏微分方程，對偏微分方程的求解思路有了深入的理解;各種解法中，效率最高的是使用EXCEL迭代求解，操作簡單并且計算速度快、精度較高。借助軟件編程的方法求解效率相對較低，并且需要對程序調(diào)試，精度方面和EXCEL求解差別不大。借助有限元求解則需要搭建有限元模型，并設置邊界條件，計算效率同樣很高，精度方面可以滿足工程需要。而解析法求解拉普拉斯方程雖然精度高，但是對算法要求較高。因此，對于偏微分方程工程方面應用，在滿足精度的前提下，解析解不一定是最優(yōu)的方法;具體需要根據(jù)實際情況選擇相應的求解方法。

參考文獻：

[1]陸金甫，關治.偏微分方程數(shù)值解法[M].清華大學出版社，2004.

[2]楊世銘，陶文銓.傳熱學[M].四版.高等教育出版社，2006.

[3]賈新民，嚴文.有限差分法求解拉普拉斯方程[J].昌吉學院學報，2009.

[4]劉大衛(wèi)，高明.正方形環(huán)域Laplace方程的簡明數(shù)值解法[J].沈陽師范大學學報，2006.