设有一长直矩形(a2a)接地金属槽,上盖板对地绝缘ф1

设有一长直矩形(a*2a)接地金属槽，上盖板对地绝缘ф1=10，侧壁，底ф2=0，求槽内的电位分布。

一．有限差分法取步长h为1，a为10，则长边与短边上的节点数分别为21和11。

具体采用简单迭代法，等位线取30条，代码和计算结果如下：function [n, v1] = FDM(a, h)%h为步长px = 2 * a / h + 1; %长边上的点数py = a / h + 1; %短边上的点数v1 = ones(py, px); %设置一10*20的矩阵v1(1, :) = ones(1, px) * 10; %初始化上边界v1(py, :) = zeros(1, px); %初始化下边界for i = 1 : py %初始化左右边界v1(i, 1) = 0;v1(i, px) = 0;endv2 = v1;maxt = 1;t = 0;n = 0;while(maxt > 1e-6)maxt = 0;n = n + 1;for i = 2 : py - 1for j = 2 : px - 1v2(i, j) = 0.25 * (v1(i, j + 1) + v1(i + 1, j) + v1(i - 1, j) + v1(i,j - 1)); %计算边界内的值t = abs(v2(i, j) - v1(i, j));if(t > maxt) %误差判断maxt = t;endendendv1 = v2;endaxis([0, 20, 0, 10, 0, 10]);contour(v1, 30);end计算的结果如下：0 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 1010 00 4.9451 6.867 7.7425 8.198 8.4589 8.6174 8.716 8.776 8.80838.8186 8.8083 8.776 8.716 8.6174 8.4589 8.198 7.7425 6.8674.9451 00 2.9134 4.7805 5.9051 6.5907 7.0201 7.2947 7.4707 7.5795 7.63887.6577 7.6388 7.5795 7.4707 7.2947 7.0201 6.5907 5.9051 4.78052.9134 00 1.9281 3.4365 4.5066 5.2395 5.7362 6.0708 6.2924 6.4325 6.50986.5345 6.5098 6.4325 6.2924 6.0708 5.7362 5.2395 4.5066 3.43651.9281 00 1.3623 2.531 3.4451 4.1245 4.6144 4.9597 5.1958 5.3482 5.43345.4608 5.4334 5.3482 5.1958 4.9597 4.6144 4.1245 3.4451 2.5311.3623 00 0.99018 1.88 2.6186 3.1988 3.6373 3.9578 4.1829 4.331 4.41484.4419 4.4148 4.331 4.1829 3.9578 3.6373 3.1988 2.6186 1.880.99018 00 0.71836 1.3804 1.9502 2.415 2.7783 3.0513 3.2472 3.378 3.45283.4772 3.4528 3.378 3.2472 3.0513 2.7783 2.415 1.9502 1.38040.71836 00 0.50289 0.973 1.3868 1.7328 2.0095 2.2218 2.3765 2.4812 2.54142.5611 2.5414 2.4812 2.3765 2.2218 2.0095 1.7328 1.3868 0.9730.50289 00 0.32019 0.62194 0.89112 1.1197 1.3054 1.4497 1.5561 1.6286 1.67061.6844 1.6706 1.6286 1.5561 1.4497 1.3054 1.1197 0.89112 0.621940.32019 00 0.15592 0.30348 0.43605 0.54959 0.6426 0.71545 0.76952 0.80655 0.828070.83513 0.82807 0.80655 0.76952 0.71545 0.6426 0.54959 0.43605 0.303480.15592 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0划出的等位线图如下：最后的迭代次数为384次。

二．有限单元法首先进行网格划分。

大致图形如下：三角元的划分如图，这样有利于以三角元的编号来计算出顶点的编号，给出1，2号三角元的顶点编号之后即可得到其他三角元的顶点编号。

然后开始计算系数矩阵K，这里没有使用迭代法。

代码如下：主函数FEM(a, gx, gy):function v = FEM(a, gx, gy) %gx gy分别为横纵方向的网格数NOE = 2 * gx * gy; %NOE为三角元总数[A, B, C, num, bound] = NaD(gx, gy, NOE);%网格划分nx = gx + 1; %横边上的节点数ny = gy + 1; %纵边上的节点数NON = nx * ny; %总节点数%计算各点在矩阵中的坐标再添加到X和Y数组中for n = 1 : nx * nyfor j = 1 : nxfor i = 1 : nyif (num(i, j) == n)X(n) = j;Y(n) = i; %输入节点坐标endendendK = zeros(NON, NON);%系数矩阵for e = 1 : NOEI = A(e);J = B(e);M = C(e);xi = (X(I) - 1) * (2 * a / gx);xj = (X(J) - 1) * (2 * a / gx);xm = (X(M) - 1) * (2 * a / gx);yi = - (Y(I) - 1) * (a / gy);yj = - (Y(J) - 1) * (a / gy);ym = - (Y(M) - 1) * (a / gy);bi = yj - ym;bj = ym - yi;bm = yi - yj;ci = xm - xj;cj = xi - xm;cm = xj - xi;s = 0.5 * (bi * cj - bj * ci);K(I, I) = K(I, I) + (bi * bi + ci * ci) / (4 * s);K(J, J) = K(J, J) + (bj * bj + cj * cj) / (4 * s);K(M, M) = K(M, M) + (bm * bm + cm * cm) / (4 * s);K(I, J) = K(I, J) + (bi * bj + ci * cj) / (4 * s);K(J, I) = K(I, J);K(I, M) = K(I, M) + (bi * bm + ci * cm) / (4 * s);K(M, I) = K(I, M);K(J, M) = K(J, M) + (bj * bm + cj * cm) / (4 * s);K(M, J) = K(J, M);endP = zeros(NON, 1);%方程右边的矩阵%强制边界条件，修改系数矩阵和方程右边的矩阵for n = 1 : nxubn = bound(1, n);P = P - 10 * K(: , ubn);end%针对上下边界进行修改for n = 1 : nxubn = bound(1, n);dbn = bound(2, n);%ubn和dbn分别为上下边界点的编号P(ubn, 1) = 10;P(dbn, 1) = 0;K(ubn, :) = zeros(1, NON);K(:, ubn) = zeros(NON, 1);K(ubn, ubn) = 1;K(dbn, :) = zeros(1, NON);K(:, dbn) = zeros(NON, 1);K(dbn, dbn) = 1;end%针对左右边界进行修改for n = 2 : ny - 1lbn = bound(3, n);rbn = bound(4, n);%lbn和rbn分别为左右边界点的编号P(lbn, 1) = 0;P(rbn, 1) = 0;K(lbn, :) = zeros(1, NON);K(:, lbn) = zeros(NON, 1);K(lbn, lbn) = 1;K(rbn, :) = zeros(1, NON);K(:, rbn) = zeros(NON, 1);K(rbn, rbn) = 1;endv1 = zeros(NON, 1);v1 = inv(K) * P;v = zeros(ny, nx);m = 0;for j = 1 : nxfor i = 1 : nym = m + 1;v(i, j) = v1(m, 1);endendcontour(v, 30);end划分网格的函数NaD(gx, gy, NOE)：function [A, B, C, num, bound] = NaD(gx, gy, NOE) %NaD函数用来划分网格gx与gy分别为横纵方向的网格数,NOE为三角元总数nx = gx + 1;ny = gy + 1; %nx,ny分别为横纵方向的结点数num(ny, nx) = 0; %用来存储节点号的数组if(ny > nx)%用来存储边界节点的数组bound(4, ny) = 0;elsebound(4, nx) = 0;end%将各个结点的编号存入num数组中for j = 1 : nxfor i = 1 : nyif (j == 1)num(i, j) = i; %第一列的结点编号elsenum(i, j) = num(i, j - 1) + ny; %非第一列的节点编号endendend%将边界结点的编号存入bound数组中，其中1，2，3，4行分别对应为上边界下边界左边界和右边界for j = 1 : nxbound(1, j) = num(1, j); %上边界bound(2, j) = num(ny, j);%下边界endfor i = 2 : ny - 1bound(3, i) = num(i, 1); %左边界bound(4, i) = num(i, nx);%右边界endA(1) = 1; B(1) = 2; C(1) = ny + 2; A(2) = 1; B(2) = ny + 2; C(2) = ny + 1;epc = 2 * gy;%计算奇数号三角元的顶点编号for e = 3 :2 : NOEfor e1 = 3 : 2 :eif (mod(A(e1 - 2), ny) == ny - 1)A(e1) = A(e1 - 2) + 2;B(e1) = B(e1 - 2) + 2;C(e1) = C(e1 - 2) + 2;elseA(e1) = A(e1 - 2) + 1;B(e1) = B(e1 - 2) + 1;C(e1) = C(e1 - 2) + 1;endendend%计算偶数号三角元的顶点编号for e = 4 :2 : NOEfor e1 = 4 : 2 :eif (mod(A(e1 - 2), ny) == ny - 1)A(e1) = A(e1 - 2) + 2;B(e1) = B(e1 - 2) + 2;C(e1) = C(e1 - 2) + 2;elseA(e1) = A(e1 - 2) + 1;B(e1) = B(e1 - 2) + 1;C(e1) = C(e1 - 2) + 1;endendendend选择a = 10，gx = 9，gy = 4，最后的计算结果如下：10 10 10 10 10 10 10 10 10 10 0 4.3692 6.0704 6.767 7.024 7.024 6.767 6.0704 4.3692 0 0 2.1151 3.4123 4.0902 4.3733 4.3733 4.0902 3.4123 2.1151 0 0 0.89632 1.5378 1.9134 2.0807 2.0807 1.9134 1.5378 0.89632 0 0 0 0 0 0 0 0 0 0 0等位线图：。