無網(wǎng)格數(shù)值求解方法學(xué)習(xí)小結(jié).doc

無網(wǎng)格數(shù)值求解方法 學(xué)習(xí)小結(jié) 一、無網(wǎng)格法的介紹有限元法存在的那些問題都來源于網(wǎng)格，在用有限元方法處理諸如金屬沖壓成型、高速沖擊、動態(tài)裂紋擴(kuò)展、流固耦合等涉及大變形和移動邊界的問題時，由于網(wǎng)格可能發(fā)生嚴(yán)重扭曲，往往需要網(wǎng)格重構(gòu)，不但精度受到了嚴(yán)重影響，計算也大幅度提高，因此有限元方法在這些領(lǐng)域的應(yīng)用遇到了困難。直接在有限元基礎(chǔ)上對其進(jìn)行改進(jìn)，效果自然不會達(dá)到最好，于是研究者把革命的對象鎖定在了網(wǎng)格上。幾經(jīng)嘗試以后，一種基于點(diǎn)集的插值方法被研究者廣泛采用，現(xiàn)今的無網(wǎng)格方法，一般就指的是這一類基于點(diǎn)集的數(shù)值方法。無網(wǎng)格方法的位移函數(shù)是在點(diǎn)的領(lǐng)域內(nèi)構(gòu)造的，并且這些區(qū)域是可以重疊的，因此在處理大變形和移動邊界等問題時，沒有網(wǎng)格的初始劃分和重構(gòu)問題，這不僅有利于這類問題計算精度的提高，還可以減少數(shù)值計算難度。目前已存在十余種無網(wǎng)格方法，它們之間的區(qū)別主要在于試函數(shù)的選擇和微分方程的等效形式。雖然無網(wǎng)格方法對于大變形和移動邊界問題具有優(yōu)勢，但其存在收斂性、數(shù)值穩(wěn)定性和效率等問題，因此無網(wǎng)格方法還只能作為有限元方法的補(bǔ)充。 無網(wǎng)格方法基本思想是將有限元法中的網(wǎng)格結(jié)構(gòu)去除，完全代之以一系列的結(jié)點(diǎn)排列。二、求解方法方法 基于位移最小二乘(MLS)近似方法EFG (Element-free Galerkin Method, Belytschko, 1994)。EFG方法計算穩(wěn)定 ，精度較高，是無網(wǎng)格方法中較為成熟的一種 方法。無網(wǎng)格法就目前來說，仍沒有有限元法發(fā)展得那么快。而且，大規(guī)模地使用無網(wǎng)格法將大大增加計算時間。因此通常只需要在那些不連續(xù)、大變形或應(yīng)力集中區(qū)域使用無網(wǎng)格法進(jìn)行離散，如沖擊區(qū)域、裂紋擴(kuò)展區(qū)域、大變形區(qū)域等，其余區(qū)域仍然可采用其他數(shù)值方法。微分方程組邊界條件等效積分形式等效積分弱形式2.1加權(quán)余量法 求解域中，若場函數(shù)是精確解，則在域中任一點(diǎn)都滿足微分方程，同時在邊界上任一點(diǎn)都滿足邊界條件式，此時等效積分形式或等效積分弱形式必然嚴(yán)格地得到滿足。但是對于復(fù)雜的實(shí)際問題，這樣的精確解往往是很難找到的，因此, 人們需要設(shè)法找到具有一定精度的近似解。設(shè)是一個近似解，即為試函數(shù)，它可以表示成為一組已知函數(shù)或Ritz基函數(shù)的線性組合，即式中為待定系數(shù)或Ritz基坐標(biāo)。將權(quán)函數(shù)代入加權(quán)余量積分式，由于系數(shù)的任意性，有上式給出了個方程。用于求解個待定系數(shù)。如果，則上式是超定的，需要借助于最小二乘法解。對上式進(jìn)行分部積分得到等效積分弱形式的近似形式2.2伽遼金法按照對權(quán)函數(shù)的不同選擇就得到不同的加權(quán)余量的計算方法并賦以不同的名稱。如果取權(quán)函數(shù)與試函數(shù)相同，則稱為Galerkin方法。我們將會看到，在很多情況下，采用伽遼金法得到的求解方程的系數(shù)矩陣是對稱的，這是在用加權(quán)余量法建立有限元格式時幾乎毫不例外地采用伽遼金法的主要原因，而且當(dāng)存在相應(yīng)的泛函數(shù)時，伽遼金法與變分法往往導(dǎo)致同樣的結(jié)果。2.3移動最小二乘近似 構(gòu)造方法：考慮求解域，其中共有個結(jié)點(diǎn)，在各個結(jié)點(diǎn)處有，但。考慮計算點(diǎn)（對于無網(wǎng)格配點(diǎn)法為結(jié)點(diǎn)；對于伽遼金無網(wǎng)格方法為高斯積分點(diǎn)），其鄰域內(nèi)的近似函數(shù)可以寫為式中：為Rits基函數(shù)，為Rits基坐標(biāo)或待求系數(shù)，是計算點(diǎn)鄰域內(nèi)任意點(diǎn)的坐標(biāo)，它包括，是基函數(shù)的個數(shù)。而， 值得注意的是，在經(jīng)典Ritz方法中， Ritz基坐標(biāo)是常數(shù)，并且基函數(shù)要滿足位移邊界條件。在式(1)中，基函數(shù)要滿足如下條件：式中：，表示在域內(nèi)具有直到階連續(xù)導(dǎo)數(shù)的函數(shù)空間。2.4邊界條件無網(wǎng)格方法的結(jié)點(diǎn)形函數(shù)多數(shù)都不滿足關(guān)系，因此位移邊界條件的處理是比較困難的。若采用緊支徑向基函數(shù)來構(gòu)造形函數(shù)，則可以像一般有限元方法那樣來處理位移邊界條件。在MLS近似中，若選奇異函數(shù)為權(quán)函數(shù)，則近似函數(shù)具有插值特性即，因此可以直接施加本質(zhì)邊界條件。對與其他情況，可以借助拉格朗日乘子方法來處理邊界條件。拉格朗日乘子法包括兩種，一種是利用邊界積分中直接引入邊界條件，即3、 具體算例 左端固定的懸臂梁，右端面受拋物線剪切載荷作用主程序：tic clear; Lx = 20; Ly = 10; young = 210; nu=0.3; q = -1;a = 0; nx = 30; ny = 20; ndivl=10; ndivw=6; dmax=2.89; Dmat = (young/(1-nu2)*1 nu 0;nu 1 0;0 0 (1-nu)/2;x,numnod,dm = mesh1(Lx,Ly,nx,ny,dmax);figurehold onplot(x(1,1:(ny+1),x(2,1:(ny+1),k-,linewidth,3);axis equal; plot(x(1,(ny+1):(ny+1):numnod),x(2,(ny+1):(ny+1):numnod),k-,linewidth,2); plot(x(1,numnod:-1:(numnod-ny),x(2,numnod:-1:(numnod-ny),k-,linewidth,2); plot(x(1,1:(ny+1):(numnod-ny),x(2,1:(ny+1):(numnod-ny),k-,linewidth,2);%plot(x(1,:),x(2,:),k.);_axis off;plot(x(1,:),x(2,:),k.);axis equal; axis off; hold offxc,conn,numcell,numq = mesh2(Lx,Ly,ndivl,ndivw);nnu,nnt,numT1,numT2 = mesh3(numq,xc,Lx,Ly,a);% nnu - % nnt - % numT1 - % numT2 - % numq-quado = 4;gauss = gauss2(quado);numq2 = numcell*quado2;gs = zeros(4,numq2); gs = egauss(xc,conn,gauss,numcell); k=kjuzhen(numnod,gs,x,dm,dmax,Dmat);rfa=400e12;ka=kajuzhen(numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa);K=k+ka;f = fjuzhen( numnod,nnt,numT2,xc,gauss,x,dm,dmax,q,Ly);%fa = zeros(2*numnod,1);%fa = fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Ly); fa = fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Lx,Ly)F=f+fa;u=zeros(2*numnod,1);for i=1:numnodu2(1,i) = u(2*i-1); u2(2,i) = u(2*i); end nx1=2; ny1=10;I = Ly3/12;for i=1:(ny1+1) xjm(1,i) = Lx/2; xjm(2,i) = -(Ly)/ny1*(i-1)+Ly; yjm(i) = -(Ly/ny1)*(i-1)+Ly/2; stress11ex(i) = -q*(Lx-xjm(1,i)* yjm(i)/I; stress12ex(i) = q/(2*I)*(Ly2/4 -yjm(i)2 ); endind = 0; enorm=0;for gg=xjm ind = ind+1; gpos = gg(1:2);v = domain(gpos,x,dm,numnod);L = length(v);en = zeros(1,2*L);phi,dphix,dphiy = shape(gpos,dmax,x,v,dm);Bmat=zeros(3,2*L);for j=1:LBmat(1:3,(2*j-1):2*j) = dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j);end for i=1:Len(2*i-1) = 2*v(i)-1;en(2*i) = 2*v(i);end stress(1:3,ind) = Dmat*Bmat*u(en); %stressex(1,ind) = ; % stressex(2,ind) = 0; % stressex(3,ind) = 0; % err = stress(1:3,ind)-stressex(1:3,ind); % err2 = weight*jac*(0.5*(inv(Dmat)*err)*(err);% enorm = enorm + err2;end%uex=zeros(2,numnod);I = Ly3/12;ind4 = 0;for i=1:numnod if(x(2,i)=Ly/2) ind4=ind4+1; uex2(ind4) = q/(6*young*I)*(3*nu*(x(2,i)-Ly/2)2*(Lx-x(1,i)+(4+5*nu)*(Ly/2)2*x(1,i)+(3*Lx-x(1,i)*x(1,i)2 ); endfigure hold onplot(x(1,(ny+1)/2:(ny+1):numnod),u2(2,(ny+1)/2:(ny+1):numnod),r.);plot(x(1,(ny+1)/2:(ny+1):numnod),uex2,-);%plot(xz,u2jy,o);xlabel(x/m,fontweight,bold);ylabel(ux/m,fontweight,bold);legend(Uynode,Exact Solution);hold off % figure % hold on% plot( xjm(2,1:(ny1+1),stress(1,1:ind),r*); % plot( xjm(2,1:(ny1+1),stress11ex(1,1:ind),.-); % legend(EFG Solution,exact solution);% % xlabel(y/m,fontweight,bold);% ylabel(Stress ,fontweight,bold);% % hold off% % figure % hold on% plot