初学matlab，用有限元算法解决一个电势分布问题

最新推荐文章于 2024-09-14 22:38:49 发布

blouwern

最新推荐文章于 2024-09-14 22:38:49 发布

阅读量894

点赞数 15

文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/2301_79878816/article/details/138620786

版权

先确定坐标系的位置，这里取小正方形的突起点中心作为原点，横轴平行于电极板，横轴取值范围有限，取长度10代替无限长

1.构造点集矩阵dots0

finextention表示网格中相邻点的间距，可自行调整，建议能被所有边界线段的长度整除

finextention=0.5;
k=1;
dotnum=(20/finextention+1)*(20/finextention+1);%大致估计点的个数，不能比实际值小
dots0=zeros(dotnum,2);
for i=-10:finextention:10
    if abs(i)>=0.5%根据横坐标分类，突起处以外的部分
        for j=-0.5:finextention:19.5
            dots0(k,:)=[i,j];
            k=k+1;
        end
    else          %突起处以上的部分       
        for j=0.5:finextention:19.5
            dots0(k,:)=[i,j];
            k=k+1;
        end
    end
end
dots0(k:dotnum,:)=[];%删除估计中多余的点（突起处以下的部分）
dotnum=k-1;%记录点的个数
disp('$ dots0 constructed');

2.三角剖分

trigs=delaunay(dots0(:,1),dots0(:,2));

disp('$ trigs constructed');

剖分的图像：

这里三角剖分出现了一点小问题，突起处以内还有三角形出现，需要删掉

for i=1:size(trigs,1)

    if i>size(trigs,1)

        break   %这里的判断语句似乎是多余的，但是在实际运行过程中不知为何i会循环到size(trigs,1)+1超出索引范围

    end

    if (abs((dots0(trigs(i,1),1)+dots0(trigs(i,2),1)+dots0(trigs(i,3),1))/3)<0.5)&&(abs((dots0(trigs(i,1),2)+dots0(trigs(i,2),2)+dots0(trigs(i,3),2))/3)<0.5)

        trigs(i,:)=[];

    end

end

%以上循环删去了位于突起处以内的点，通过三角形重心进行判断

这个是根据最初结果进行修改，在其他的问题中，如果要一步到位的话，应该也可以用类似这样的方法

最后的结果如下

这里三角剖分直接使用了函数delaunay，实际上简单的图形应该可以自行划分，我自己原来写了一个，但是需要花费很长一段运行时间

%此次的问题示例中图形形状简单，可以直接先用矩形网格剖分，每个小矩形用对角线剖分,

%并且用矩阵T储存每个三角形顶点坐标

T=[]

l=1

for i=1:dotnum

   for j=1:dotnum

       if dots0(j,1)-dots0(i,1)==finextn

           T(l,1)=i;

           T(l,2)=j;

           for k=1:dotnum

               if dots0(k,2)-dots0(i,2)==finextn

                   T(l,3)=k;

                   l=l+1;

                end

            end

        end

    end

end

3.电势条件初始化

进行的同时用orderstorage记录已知点的序号，为后续通过索引将矩阵分块做准备

len1=0;

orderstorage=zeros(dotnum,1);

U=zeros(dotnum,1);

for k=1:dotnum

    %上下极板平整处

    if dots0(k,2)==19.5

        U(k)=1;

        len1=len1+1;

        orderstorage(len1)=k;

    end

    if dots0(k,2)==-0.5

        U(k)=0;

        len1=len1+1;

        orderstorage(len1)=k;

    end

    %下级板突起处

    if (dots0(k,1)==-0.5 ||dots0(k,1)==0.5)&&(-0.5<=dots0(k,2)&&dots0(k,2)<=0.5)

        U(k)=0;

        len1=len1+1;

        orderstorage(len1)=k;

    end

    if (dots0(k,2)==0.5) && (-0.5<=dots0(k,1)&&dots0(k,1)<=0.5)

        U(k)=0;

        len1=len1+1;

        orderstorage(len1)=k;

    end

    %电极板左右两边电势可视为匀强电场

    if (dots0(k,1)==-10 ||dots0(k,1)==10)

        U(k)=(dots0(k,2)+0.5)/20;

        len1=len1+1;

        orderstorage(len1)=k;

    end

end

orderstorage=unique(orderstorage,'stable');%删掉for循环中重复判断的点

orderstorage(find(orderstorage==0))=[];%删掉多余的点,也可以用any函数

len1=size(orderstorage,1);%len1确定为所有不重复的已知点的个数

disp('$ electric potential initialized');

4.构建刚度矩阵K和力矩阵F

epsilon=1;

trigarea=zeros(size(trigs,1),1);

K0=zeros(dotnum);

b=zeros(dotnum,size(trigarea,1));

c=zeros(dotnum,size(trigarea,1));

for i=1:size(trigarea,1)

    for k=1:3           %k,l,m轮换

        l=k+1;if l>3

                 l=l-3;

              end

        m=k+2;if m>3

                 m=m-3;

              end

        b(trigs(i,k),i)=dots0(trigs(i,l),2)-dots0(trigs(i,m),2);

        c(trigs(i,k),i)=dots0(trigs(i,l),1)-dots0(trigs(i,m),1);

    end

    trigarea(i)=abs((b(trigs(i,1),i)*c(trigs(i,2),i)-b(trigs(i,2),i)*c(trigs(i,1),i))/2);

    for n=1:3

        for m=1:3

            K0(trigs(i,n),trigs(i,m))=K0(trigs(i,n),trigs(i,m))+epsilon*(b(trigs(i,n),i)*b(trigs(i,m),i)+c(trigs(i,n),i)*c(trigs(i,m),i))/4/trigarea(i);

        end

    end

end

F=zeros(dotnum,1);%这里F事实上不是全零的，只有在电极板以内点的序号对应的元素是全零的

disp('$ matrix K&phro constructed');

5.所有矩阵和向量重排顺序

已知未知分成两份，注意dots0也需要进行换序，否则无法一一对应

U1=zeros(len1,1);%U1用来储存已知点的电势

K1=zeros(dotnum,len1);

dots1=zeros(len1,2);

i1=1;

K2=zeros(dotnum,dotnum-len1);

dots2=zeros(dotnum-len1,2);

i2=1;

for i=1:len1

    U1(i)=U(orderstorage(i));

end

for i=1:dotnum

    if ismember(i,orderstorage)

       K1(:,i1)=K0(:,i);

       dots1(i1,:)=dots0(i,:);

       i1=i1+1;

    else

        K2(:,i2)=K0(:,i);

        dots2(i2,:)=dots0(i,:);

        i2=i2+1;

    end

end

%第一次换序，K换列，U（实际上是dots）换行

K0(:,1:len1)=K1;

K0(:,len1+1:dotnum)=K2;  

dots=zeros(dotnum,2);

dots(1:len1,:)=dots1;

dots(len1+1:dotnum,:)=dots2;

%第二次换序，K换行，phro换行

K3=zeros(len1,dotnum);

i1=1;

K4=zeros(dotnum-len1,dotnum);

i2=1;

for i=1:dotnum

    if ismember(i,orderstorage)

       K3(i1,:)=K0(i,:);

       i1=i1+1;

    else

        K4(i2,:)=K0(i,:);

        i2=i2+1;

    end

end

K0(1:len1,:)=K3;

K0(len1+1:dotnum,:)=K4;  

disp('$ reorder completed');

以上的动作相当于把K进行一系列行列对称的变换，同时U和F也进行了变换，只是仅储存了U1，U2是未知的，后续直接通过方程求解，

而对于F，它的边界点序号对应的元素（即P1部分）值不是零，无法确定，但是P2是可以知道全是零

6.K和F分块

求解方程

$(K)(\Phi) = (p)$

$\left. \Longleftrightarrow\binom{\left( K_{11} \right)\left( K_{12} \right)}{\left( K_{12} \right)\left( K_{22} \right)}\binom{\left( \Phi_{1} \right)}{\left( \Phi_{2} \right)} = \binom{\left( p_{1} \right)}{\left( p_{2} \right)} \right.$

矩阵（p）即程序中的F，（Φ）即矩阵中的U，这里我们可以确定P2全为零，事实上矩阵分块之后是有两个矩阵方程的，

$\left( K_{11} \right)\left( \Phi_{1} \right) + \left( K_{12} \right)\left( \Phi_{2} \right) = \left( p_{1} \right)\ldots\ldots{（1）}$

$\left( K_{21} \right)\left( \Phi_{1} \right) + \left( K_{22} \right)\left( \Phi_{2} \right) = \left( p_{2} \right)\ldots\ldots{（2）}$

但含有方程（1）中由于P1是未知的，不为零，故难以求出Φ2，因此目前只能选方程（2），此方程中只有Φ2是未知的，且K22是正定矩阵，也就是它可逆

K=mat2cell(K0,[len1 dotnum-len1],[len1 dotnum-len1]);

P=mat2cell(F,[len1 dotnum-len1]);

U2=linsolve(K{2,2},(P{2}-K{2,1}*U1));

U=[U1;U2];%求得的U2和已知的U1合并

disp('$ all potentail computed');

7．方程解U转化成点点对应的电势矩阵phi

phi=zeros(20/finextention+1,20/finextention+1);

for i=-10:finextention:10

    for j=-0.5:finextention:19.5

        for k=1:dotnum

            if dots(k,:)==[i,j]

                phi((j+0.5)/finextention+1,(i+10)/finextention+1)=U(k);%这里要注意坐标顺序，行对应的是纵坐标，列对应横坐标

            end

        end

    end

end

disp('$ potentail data meshlized');

这一步是整个程序中最花时间的地方，把列向量转化成为方阵，因为把矩阵顺序重排了，所以不得不使用for循环遍历所有点

8.画图

contour(-10:finextention:10, -0.5:finextention:19.5, phi,30);

title('Electric Potential Distribution');%等势线

hold on;

poi1=[-10-0.5i,-0.5-0.5i,-.5+0.5i,0.5+0.5i,0.5-0.5i,10-0.5i];

ang1=[0,1/2,3/2,3/2,1/2,0];

p1=polygon(poi1,ang1);%下极板

poi2=[-10+19.5i,19.5i,10+19.5i];

ang2=[0,1,0];

p2=polygon(poi2,ang2);%上极板

%这里画电极板用的是polygon函数，把两个极板的一维投影看成了多边形，应该还有更好的画极板的方式

plot(p1);

hold on;

plot(p2);

hold on;

%最后画出来的颜色是可以编辑的，可以根据需求自行添加参数

disp('$ process finished');

运行结果：