设有一长直矩形(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

2、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 - 1 for j = 2 : px - 1 v2(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) %误差判断 m

3、axt = t; end end end v1 = v2;endaxis(0, 20, 0, 10, 0, 10);contour(v1, 30);end计算的结果如下：010101010101010101010101010101010101010004.94516.8677.74258.1988.45898.61748.7168.7768.80838.81868.80838.7768.7168.61748.45898.1987.74256.8674.9451002.91344.78055.90516.59077.02017.29477.47077.57957.63887.65777.6388

4、7.57957.47077.29477.02016.59075.90514.78052.9134001.92813.43654.50665.23955.73626.07086.29246.43256.50986.53456.50986.43256.29246.07085.73625.23954.50663.43651.9281001.36232.5313.44514.12454.61444.95975.19585.34825.43345.46085.43345.34825.19584.95974.61444.12453.44512.5311.3623000.990181.882.61863.1

5、9883.63733.95784.18294.3314.41484.44194.41484.3314.18293.95783.63733.19882.61861.880.99018000.718361.38041.95022.4152.77833.05133.24723.3783.45283.47723.45283.3783.24723.05132.77832.4151.95021.38040.71836000.502890.9731.38681.73282.00952.22182.37652.48122.54142.56112.54142.48122.37652.22182.00951.73

6、281.38680.9730.50289000.320190.621940.891121.11971.30541.44971.55611.62861.67061.68441.67061.62861.55611.44971.30541.11970.891120.621940.32019000.155920.303480.436050.549590.64260.715450.769520.806550.828070.835130.828070.806550.769520.715450.64260.549590.436050.303480.155920000000000000000000000划出的

7、等位线图如下：最后的迭代次数为384次。二 有限单元法首先进行网格划分。大致图形如下：1234三角元的划分如图，这样有利于以三角元的编号来计算出顶点的编号，给出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

8、 + 1; %纵边上的节点数NON = nx * ny; %总节点数%计算各点在矩阵中的坐标再添加到X和Y数组中for n = 1 : nx * ny for j = 1 : nx for i = 1 : ny if (num(i, j) = n) X(n) = j; Y(n) = i; %输入节点坐标 end end endendK = zeros(NON, NON);%系数矩阵for e = 1 : NOE I = A(e); J = B(e); M = C(e); xi = (X(I) - 1) * (2 * a / gx); xj = (X(J) - 1) * (2 * a / gx)

9、; 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) =

10、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)

11、;endP = zeros(NON, 1);%方程右边的矩阵%强制边界条件，修改系数矩阵和方程右边的矩阵for n = 1 : nx ubn = bound(1, n); P = P - 10 * K(: , ubn);end%针对上下边界进行修改for n = 1 : nx ubn = 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(

12、dbn, :) = zeros(1, NON); K(:, dbn) = zeros(NON, 1); K(dbn, dbn) = 1;end%针对左右边界进行修改for n = 2 : ny - 1 lbn = 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) = ze

13、ros(NON, 1); K(rbn, rbn) = 1;endv1 = zeros(NON, 1);v1 = inv(K) * P;v = zeros(ny, nx);m = 0;for j = 1 : nx for i = 1 : ny m = 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

14、 = gy + 1; %nx,ny分别为横纵方向的结点数num(ny, nx) = 0; %用来存储节点号的数组if(ny > nx)%用来存储边界节点的数组 bound(4, ny) = 0; else bound(4, nx) = 0;end%将各个结点的编号存入num数组中for j = 1 : nx for i = 1 : ny if (j = 1) num(i, j) = i; %第一列的结点编号 else num(i, j) = num(i, j - 1) + ny; %非第一列的节点编号 end endend%将边界结点的编号存入bound数组中，其中1，2，3，4行分别对

15、应为上边界下边界左边界和右边界for j = 1 : nx bound(1, j) = num(1, j); %上边界 bound(2, j) = num(ny, j);%下边界 endfor i = 2 : ny - 1 bound(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 : NOE for e1

16、 = 3 : 2 :e if (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; else A(e1) = A(e1 - 2) + 1; B(e1) = B(e1 - 2) + 1; C(e1) = C(e1 - 2) + 1; end endend%计算偶数号三角元的顶点编号for e = 4 :2 : NOE for e1 = 4 : 2 :e if (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; else A(e1) = A(e1 - 2) + 1; B(e1) = B(e1 - 2) + 1; C(e1) = C(