//The scilab code is developed by //B. Sai Ram, Research Scholar, Dept. of Electrical Engineering, IIT Bombay //for the MOOC course: Electrical Equipment and Machines: Finite Element Analysis //by Prof. S. V. Kulkarni, IIT Bombay. clc; clear; n=51; // number of nodes along x-axis m=51; // number of nodes along y-axis n_nodes=(n*m); // total number of nodes in the domain disp('intializing global coefficent matrix'); C=zeros(n_nodes,n_nodes); BJ=zeros(n_nodes,1); x=0:0.1/(n-1):0.1; // nodal points on x-axis y=0:0.1/(m-1):0.1; // nodal points on y- axis // this forms a square with 0.1 X 0.1 dimensions in MKS system disp('forming p matrix'); // formation of p matrix // Given domain is discritized into squares and then each square is // divided into 2 trainguler elements using diagonals of square h=1; // h indicates node number (column index of p matrix) for i=1:m // index for y for j=1:n // index for x // for each y 51 x coordinates are formed p(1,h)=x(j); p(2,h)=y(i); h=h+1; end end disp('----------------------------------------------------------------'); disp('forming t matrix'); // formation of t matrix, 3 gobal node numbers for each element and sub-domain number h=1;// here, h indicates element number (column index of t matrix) // forming t matix, of an element with an edge parallel to x-axis, for i=1:m-1 for j=1:n-1 t(2,h)=j+((i-1)*n); // for i = 1, j =1 , h =1, n=51, t(1,1) = 1 t(3,h)=j+1+((i-1)*n);// t(2,1) = 2 t(4,h)=j+1+((i)*n);// t(3,1) = 53 // x coordinates of each node of the element Xc(1) = p(1,j+((i-1)*n)); Xc(2) = p(1,j+1+((i-1)*n)); Xc(3) = p(1,j+1+((i)*n)); // y coordinates of each node of the element Yc(1) = p(2,j+((i-1)*n)); Yc(2) = p(2,j+1+((i-1)*n)); Yc(3) = p(2,j+1+((i)*n)); // calculating the centroid of the element cX = sum(Xc)/3; cY = sum(Yc)/3; if (cX<=0.07 && cX>=0.03) && (cY<=0.07 && cY>=0.03) then t(1,h)=2; // same material everywhere else t(1,h)=1; // same material everywhere end h=h+1;// incrementing the element number end end for i=1:m-1 for j=1:n-1 t(2,h)=j+((i-1)*n);// for i = 1, j =1 , h =1, n=51, t(1,1) = 1 t(3,h)=j+1+((i)*n);// t(2,1) = 53 t(4,h)=j+((i)*n);// t(3,1) = 52 // x coordinates of each node of the element Xc(1) = p(1,j+((i-1)*n)); Xc(2) = p(1,j+1+((i)*n)); Xc(3) = p(1,j+((i)*n)); // y coordinates of each node of the element Yc(1) = p(1,j+((i-1)*n)); Yc(2) = p(1,j+1+((i)*n)); Yc(3) = p(1,j+((i)*n)); // calculating the centroid of the element cX = sum(Xc)/3; cY = sum(Yc)/3; if (cX<=0.07 && cX>=0.03) && (cY<=0.07 && cY>=0.03) then t(1,h)=2; // same material everywhere else t(1,h)=1; // same material everywhere end h=h+1;// incrementing the element number end end disp('----------------------------------------------------------'); [A3 n_elements]=size(t);// A3=4, //size of t matrix is 4 X no. of elements and 4 gets assigned to A3 [A1 n_nodes]=size(p);// A1=2 //size of t matrix is 2 X no. of nodes and 2 gets assigned to A1 disp('number of elements'); disp(n_elements); //for n = 51, m = 51, n_elements = 5000 disp('number of nodes'); disp(n_nodes); //for n = 51, m = 51, n_nodes = 2601 //Forming element coeffcient matrices disp('forming element coefficient matrix'); for element=1:n_elements nodes=t(2:4,element);// 1st three rows --> 3 global node numbers of each element //'nodes' is a 3 X 1 matrix containing global // node numbers of element under consideration Xc=p(1,nodes'); //x-coordinates of nodes, Xc is a 3 X 1 matrix Yc=p(2,nodes'); //y-coordinates of nodes, Yc is a 3 X 1 matrix P=zeros(3,1); Q=zeros(3,1); P(1)=Yc(2)-Yc(3); P(2)=Yc(3)-Yc(1); P(3)=Yc(1)-Yc(2); Q(1)=Xc(3)-Xc(2); Q(2)=Xc(1)-Xc(3); Q(3)=Xc(2)-Xc(1); //finding x and y coordinates //of the centroid //see if the element lies within the //rectangular conductor cX = sum(Xc)/3; cY = sum(Yc)/3; if t(1,element) == 2 then Mu(element) = 4*%pi*1e-7; bj = 1e3; // J--> current density in the conductor else Mu(element) = 4*%pi*1e-7; bj = 0; end //eq. (15.21c) of //M. N. O. Sadiku and S. V. Kulkarni, Principles of Electromagnetics, //Sixth Edition, Oxford University Press, India, 2015 // area of triangular element delta= 0.5*abs((P(2)*Q(3))-(P(3)*Q(2)));//Absolute value is taken since //the three nodes may not have been numbered in the anticlockwise for i=1:3 for j=1:3 // 3 X 3 element coefficient matrix c(element,i,j)=((P(i)*P(j))+(Q(i)*Q(j)))/(4*delta*Mu(element)); end end // 3 X 1 element level source matrix be(element,:) = bj*(delta/3)*[1;1;1]; end disp('forming global coefficient matrix'); for element=1:n_elements nodes=t(2:4,element); // Formation of the global coefficient matrix for i=1:3 for j=1:3 C(nodes(i),nodes(j))=C(nodes(i),nodes(j))+c(element,i,j); end BJ(nodes(i),1) = be(element,i)+ BJ(nodes(i),1); end end disp('applying boundary conditions'); // e matrix is not formed hence boundary conditions are applied manually // 1st edge, bottom edge for i=1:n C(i,:)=zeros(1,n_nodes); C(i,i)=1; BJ(i,1)=0; end // 2nd edge, vertical left side edge for i=n+1:n:(n*m)-2*n+1 C(i,:)=zeros(1,n_nodes); C(i,i)=1; BJ(i,1)=0; end // 3rd edge, vertical right side edge for i=2*n:n:(n*m) C(i,:)=zeros(1,n_nodes); C(i,i)=1; BJ(i,1)=0; end // 4th edge, top edge for i=(n*(m-1))+1:(n*m)-1 C(i,:)=zeros(1,n_nodes); C(i,i)=1; BJ(i,1)=0; end disp('computing nodal potentials'); A=C\BJ; //Here, A is a column vector - it is converted into n X m matrix, //similarly, two rows of p matrix (x and y coordinates) are converted //into n X m matrix Vn=matrix(A(1:n*m),n,m);//potentials at each node of n X m grid XX=matrix(p(1,1:n*m),n,m);// x coordinates of n X m grid YY=matrix(p(2,1:n*m),n,m);// y coordinates of n X m grid figure(1) surf(XX,YY,Vn') figure(2) contour(x,y,Vn,25)