实验三二维电流场有限元分析实验

1《生物医学电磁场数值分析》实验报告实验题目实验三二维电流场有限元分析实验姓名郭立杰学号085958班级生医C081组别6时间2011-5-23地点7A-705指导教师李颖同组人王佳一、实验目的掌握二维电流场有限元分析的方法，编制相应程序，包括有限元系数矩阵的生成，边界条件的处理，方程组的求解等，培养解决实际电磁场问题的能力。二、实验条件硬件：微型计算机。软件：MATLAB6.0软件。三、实验内容1、针对二维电流场，根据变分原理或加权余量法推导其有限元方程形式。2、编制其有限元分析程序，并进行求解。四、实验步骤1、预习：推导生物医学电磁场中二维稳态电流场的有限元离散方程。2、根据场域剖分的结果，利用三角形单元的形状函数求解系数矩阵。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部分区间分别为[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个节点的横、纵坐标p=p+1;endend2endxxx=x1(i+1);forn=1:num_nodeynode(p,1)=xxx;%最后一列节点的横坐标node(p,2)=y1(1)+stepy*(n-1);%最后一列节点的纵坐标p=p+1;endp=1;fori=1:num_nodex-1forj=1:num_nodey-1%a类单元的K,M,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节点编号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;%二维电流场场域内电极的节点位置refe=60;S(refe,:)=0;S(:,refe)=0;S(refe,refe)=1;%边界条件的处理U=inv(S)*F;%求节点电位figure;plot(U,'b.-');%绘制节点电位变化曲线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;功能函数：functions=qxishuzhen(node,element,rou)n=size(node,1);3s=sparse(n,n,0);m=size(element,1);fori=1:ms1=sparse(n,n,0);con=1/rou(i);%σ=1/ρK=element(i,1);M=element(i,2);N=element(i,3);%第i个单元三个节点的标号依次为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)=(q1^2+r1^2);s1(M,M)=(q2^2+r2^2);s1(N,N)=(q3^2+r3^2);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五、实验结果1、二维电流场场域内节点电位图020406080100120-2-1.5-1-0.500.511.522、二维电流场场域内电位分布灰度图4051015012345678-1-0.500.51六、分析与讨论在本次实验中通过对实验程序的编写和修改，使我彻底理解了总体合成过程的步骤和意义，同时还让我体会到使用MATLAB编程来实现有限元分析与理论知识的差别之处，也让我第一次真正地体会了到有限元分析方法的现实使用价值。七、思考题1、有限元系数矩阵如何进行单元分析和总体合成？答：单元分析：对节点和单元进行编号离散化，利用线性插值得到基函数，然后得到单元特征式。总体合成：求得各单元的导数式，将它们集合起来，即可装配成总体代数方程组。2、为什么说第二类边界条件是自然边界条件，而第一类边界条件是强加边界条件？答：第二类边界条件上的节点不需要做任何特殊的处理，它决定于边界上介质的情况；而第一类边界条件是直接规定的，所以说第二类边界条件是自然边界条件，而第一类边界条件是强加边界条件。3、如何处理有限元方程的第一类边界条件？答：主要有四种方法：（1）将离散的强加边界条件bbuu0加到矩阵方程式nnijFFFuuuS2121上，相当于方程组中一部分未知量bu中有已知值，因此只要从上式中删去编号i=b的bn个方程，而在余下的方程中，将bbuu0的已知数值代入并将相应项移到右端，再去解余下的bnn个方程组即可。（2）对系数矩阵S和右端项阵列F的第b行作一些修改，将主对角线上的元素bbS改为一个足够大的正数，相应的将右端的bF改为bEu0，这样第b个方程中除bu外的其他未知量的系数与bu的系数相比微不足道，第b个方程近似于bbEuEu0。（3）将系数矩阵住对角线元素bbS置1，其余第b行和第b列的元素均置0，同时将5右端项列阵第b行置为bu0，其余元素改为bbibiuSF0。（4）对于简单的情况，边界上部分电位值为0的情况，不把这些边界节点列入方程组中，方程组的阶数仅决定于内部节点的数目，在节点编号时，将边界节点放在最后，由于边界上节点的u值为0，因此在计算系数矩阵式与边界节点有关的所有系数均不必计算。