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

《生物醫(yī)學(xué)電磁場(chǎng)數(shù)值分析》實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)題目實(shí)驗(yàn)三二維電流場(chǎng)有限元分析實(shí)驗(yàn)姓名郭立杰學(xué)號(hào)085958班級(jí)生醫(yī)C081組別6時(shí)間2011-5-23地點(diǎn)7A-705指導(dǎo)教師李穎同組人王佳一、實(shí)驗(yàn)?zāi)康恼莆斩S電流場(chǎng)有限元分析的方法，編制相應(yīng)程序，包括有限元系數(shù)矩陣的生成，邊界條件的處理，方程組的求解等，培養(yǎng)解決實(shí)際電磁場(chǎng)問(wèn)題的能力。二、實(shí)驗(yàn)條件硬件：微型計(jì)算機(jī)。軟件：MATLAB6.0軟件。三、實(shí)驗(yàn)內(nèi)容1、針對(duì)二維電流場(chǎng)，根據(jù)變分原理或加權(quán)余量法推導(dǎo)其有限元方程形式。2、編制其有限元分析程序，并進(jìn)行求解。四、實(shí)驗(yàn)步驟1、預(yù)習(xí)：推導(dǎo)生物醫(yī)學(xué)電磁場(chǎng)中二維穩(wěn)態(tài)電流場(chǎng)的有限元離散方程。2、根據(jù)場(chǎng)域剖分的結(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軸上的步長(zhǎng)分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個(gè)節(jié)點(diǎn)的橫、縱坐標(biāo)p=p+1;endendendxxx=x1(i+1);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)編號(hào)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)編號(hào)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;%二維電流場(chǎng)場(chǎng)域內(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);%繪制整個(gè)場(chǎng)域的電位分布holdon;endcolormap(gray);%繪制電位分布圖為灰度圖像colorbar;功能函數(shù)：functions=qxishuzhen(node,element,rou)n=size(node,1);

s=sparse(n,n,0);m=size(element,1);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個(gè)單元三個(gè)節(jié)點(diǎn)的標(biāo)號(hà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、二維電流場(chǎng)場(chǎng)域內(nèi)節(jié)點(diǎn)電位圖21.510.50-0.5-1-1.5-20204060801001202、二維電流場(chǎng)場(chǎng)域內(nèi)電位分布灰度圖

0.5-0.5-1六、分析與討論0.5-0.5-1在本次實(shí)驗(yàn)中通過(guò)對(duì)實(shí)驗(yàn)程序的編寫和修改，使我徹底理解了總體合成過(guò)程的步驟和意義，同時(shí)還讓我體會(huì)到使用MATLAB編程來(lái)實(shí)現(xiàn)有限元分析與理論知識(shí)的差別之處，也讓我第一次真正地體會(huì)了到有限元分析方法的現(xiàn)實(shí)使用價(jià)值。七、思考題1、有限元系數(shù)矩陣如何進(jìn)行單元分析和總體合成？答：?jiǎn)卧治觯簩?duì)節(jié)點(diǎn)和單元進(jìn)行編號(hào)離散化，利用線性插值得到基函數(shù)，然后得到單元特征式。總體合成：求得各單元的導(dǎo)數(shù)式，將它們集合起來(lái)，即可裝配成總體代數(shù)方程組。2、為什么說(shuō)第二類邊界條件是自然邊界條件，而第一類邊界條件是強(qiáng)加邊界條件？答：第二類邊界條件上的節(jié)點(diǎn)不需要做任何特殊的處理，它決定于邊界上介質(zhì)的情況；而第一類邊界條件是直接規(guī)定的，所以說(shuō)第二類邊界條件是自然邊界條件，而第一類邊界條件是強(qiáng)加邊界條件。3、如何處理有限元方程的第一類邊界條件？答：主要有四種方法：（1）將離散的強(qiáng)加邊界條件‘=U一、.一、0b加到矩陣方程式"]ijU1u2F1F2上，相當(dāng)于方程組中一部分未知量^b中有已知值，因此只要從上式中刪去編號(hào)i=b的nb個(gè)方程，而在余下的方程中，將Ub=U0b的已知數(shù)值代入并將相應(yīng)項(xiàng)移到右端，再去解余下的〃—nb個(gè)方程組即可。（2）對(duì)系數(shù)矩S陣S和右端項(xiàng)陣列F的第b行作一些修改，將主對(duì)角線上的元素bb改為一個(gè)足夠大的正數(shù)，相應(yīng)的將右端的Fb改為EU0b，這樣第b個(gè)方程中除Ub外的其他未知量的系數(shù)與Ub的系數(shù)相比微不足道，第b個(gè)方程近似于EUb=EU0b。（3）將系數(shù)矩陣住對(duì)角線元素Sbb置1，其余第b行和第