實(shí)驗(yàn)三二維電流場有限元分析實(shí)驗(yàn)

《生物醫(yī)學(xué)電磁場數(shù)值分析》實(shí)驗(yàn)報告實(shí)驗(yàn)題目實(shí)驗(yàn)三二維電流場有限元分析實(shí)驗(yàn)姓名郭立杰學(xué)號085958班級生醫(yī)C081組別6時間2011-5-23地點(diǎn)7A-705指導(dǎo)教師李穎同組人王佳一、 實(shí)驗(yàn)?zāi)康恼莆斩S電流場有限元分析的方法，編制相應(yīng)程序，包括有限元系數(shù)矩陣的生成，邊界條件的處理,方程組的求解等，培養(yǎng)解決實(shí)際電磁場問題的能力。二、 實(shí)驗(yàn)條件硬件：微型計算機(jī)。軟件：MATLAB6.0軟件。三、 實(shí)驗(yàn)內(nèi)容1、 針對二維電流場，根據(jù)變分原理或加權(quán)余量法推導(dǎo)其有限元方程形式。2、 編制其有限元分析程序，并進(jìn)行求解。四、 實(shí)驗(yàn)步驟1、 預(yù)習(xí)：推導(dǎo)生物醫(yī)學(xué)電磁場中二維穩(wěn)態(tài)電流場的有限元離散方程。2、 根據(jù)場域剖分的結(jié)果，利用三角形單元的形狀函數(shù)求解系數(shù)矩陣。3、 處理第二類和第一類邊界條件。4、 求解有限元線性方程組。主程序：clearall;nsx=3;nsy=1;nsx1=[3,4,5]; %將x軸分3部分，第一部分分3小份,第二部分分成4份，第三部分分成5份nsy1=8;x1=[0,5,10,15];%x軸的范圍為[0,12],所分的3部分區(qū)間分別為[0,3][3,8][8,12]y1=[0,8]；num_nodex=sum(nsx1)+1;num_nodey=sum(nsy1)+1;total_node=num_nodex*num_nodey;total_element=(num_nodex-1)*(num_nodey-1)*2;stepy=(y1(2)-y1(1))/nsy1;p=1;fori=1:nsxstepx(i)=(x1(i+1)-x1(i))/nsx1(i);%x軸上的步長分3種情況form=1:nsx1(i)xx=x1(i)+stepx(i)*(m-1);forj=1:nsy1+1node(p,1)=xx;node(p,2)=y1(1)+stepy*(j-1);%第p個節(jié)點(diǎn)的橫、縱坐標(biāo)p=p+1;endendendxxx=xl(i+l);forn=1:num_nodeynode(p,1)=xxx; %最后一列節(jié)點(diǎn)的橫坐標(biāo)node(p,2)=y1(1)+stepy*(n-1);%最后一列節(jié)點(diǎn)的縱坐標(biāo)p=p+1；endp=1;fori=1:num_nodex-1forj=1:num_nodey-1%a類單元的K,M,N節(jié)點(diǎn)編號element(p,1)=num_nodey*(i-1)+j;element(p,2)=element(p,1)+1+num_nodey;element(p,3)=element(p,1)+1;p=p+1;endforj=1:num_nodey-1 %b類單元的K,M,N節(jié)點(diǎn)編號element(p,1)=num_nodey*(i-1)+j;element(p,2)=element(p,1)+num_nodey;element(p,3)=element(p,2)+1;p=p+1;endendrou(1:total_element)=1;S=qxishuzhen(node,element,rou); %求稀疏矩陣F(total_node,1)=0;F(6,1)=1;F(111,1)=-1; %二維電流場場域內(nèi)電極的節(jié)點(diǎn)位置refe=60;S(refe,:)=0;S(:,refe)=0;S(refe,refe)=1;%邊界條件的處理U=inv(S)*F;%求節(jié)點(diǎn)電位figure;plot(U,'b.-');%繪制節(jié)點(diǎn)電位變化曲線figure;fori=1:total_elementp1=element(i,1);p2=element(i,2);p3=element(i,3);x1=node(p1,1);x2=node(p2,1);x3=node(p3,1);y1=node(p1,2);y2=node(p2,2);y3=node(p3,2);xx=[x1,x2,x3];yy=[y1,y2,y3];zz=(U(p1,1)+U(p2,1)+U(p3,1))/3;%利用線性插值方法求出單元中的電位分布fill(xx,yy,zz); %繪制整個場域的電位分布holdon;endcolormap(gray); %繪制電位分布圖為灰度圖像colorbar;功能函數(shù)：functions=qxishuzhen(node,element,rou)n=size(node,1);

s=sparse(n,n,O);m=size(element,l);fori=1:ms1=sparse(n,n,0);con=1/rou(i); %O=1/pK=element(i,1);M=element(i,2);N=element(i,3);%第i個單元三個節(jié)點(diǎn)的標(biāo)號依次為K,M,Nrk=node(K,1:2);rm=node(M,1:2);rn=node(N,1:2);x1=rk(1,1);y1=rk(1,2);x2=rm(1,1);y2=rm(1,2);x3=rn(1,1);y3=rn(1,2);P1=x2*y3-x3*y2;p2=x3*y1-x1*y3;p3=x1*y2-x2*y1;q1=y2_y3;q2=y3_y1;q3=y1_y2;r1=x3-x2;r2=x1-x3;r3=x2-x1;mianji=((y2-y3)*(x1-x3)-(y3-y1)*(x3-x2))/2;mianji1=4*mianji;s1(K,K)=(q1A2+r1A2);s1(M,M)=(q2A2+r2A2);s1(N,N)=(q3A2+r3A2);s1(K,M)=(q1*q2+r1*r2);s1(K,N)=(q1*q3+r1*r3);s1(M,N)=(q2*q3+r3*r2);s1(M,K)=s1(K,M);s1(N,K)=s1(K,N);s1(N,M)=s1(M,N);s1=con*s1/mianji1;s=s+s1;end五、實(shí)驗(yàn)結(jié)果1、二維電流場場域內(nèi)節(jié)點(diǎn)電位圖2、二維電流場場域內(nèi)電位分布灰度圖0.5-0.5-10.5-0.5-1六、分析與討論在本次實(shí)驗(yàn)中通過對實(shí)驗(yàn)程序的編寫和修改，使我徹底理解了總體合成過程的步驟和意義，同時還讓我體會到使用MATLAB編程來實(shí)現(xiàn)有限元分析與理論知識的差別之處，也讓我第一次真正地體會了到有限元分析方法的現(xiàn)實(shí)使用價值。七、思考題1、 有限元系數(shù)矩陣如何進(jìn)行單元分析和總體合成？答：單元分析：對節(jié)點(diǎn)和單元進(jìn)行編號離散化，利用線性插值得到基函數(shù)，然后得到單元特征式?？傮w合成：求得各單元的導(dǎo)數(shù)式，將它們集合起來，即可裝配成總體代數(shù)方程組。2、 為什么說第二類邊界條件是自然邊界條件，而第一類邊界條件是強(qiáng)加邊界條件？答：第二類邊界條件上的節(jié)點(diǎn)不需要做任何特殊的處理，它決定于邊界上介質(zhì)的情況；而第一類邊界條件是直接規(guī)定的，所以說第二類邊界條件是自然邊界條件，而第一類邊界條件是強(qiáng)加邊界條件。3、如何處理有限元方程的第一類邊界條件?答：主要有四種方法：（1）將離散的強(qiáng)加邊界條件役U0b答：主要有四種方法：（1）將離散的強(qiáng)加邊界條件役U0b加到矩陣方程式u—F_L —.11ruF-(/12—2uFn n上，相當(dāng)TOC\o"1-5"\h\zu n于方程組中一部分未知量b中有已知值，因此只要從上式中刪去編號i=b的b個方程，而在余下的方程7/—7/ —中，將b_0b的已知數(shù)值代入并將相應(yīng)項移到右端，再去解余下的 b個方程組即可。（2）對系數(shù)矩S陣S和右端項陣列F的第b行作一些修改，將主對角線上的元素bb改為一個足夠大的正數(shù)，相應(yīng)的將右FEu u u端的b改為 0b，這樣第b個方程中除b外的其他未知量的系數(shù)與b的系數(shù)相比微不足道，第b個方程Eu—Eu S近似于b~ 0b。（3）將系數(shù)矩陣住對角線元素bb置1,其余第b行