205 lines
5.5 KiB
Matlab
205 lines
5.5 KiB
Matlab
% cemLaplace01.m
|
|
|
|
% Solution of Laplace's equation
|
|
clear
|
|
clc
|
|
tic
|
|
|
|
% INPUTS ================================================================
|
|
% Number of XY grid points (integer)
|
|
Nx = 101; Ny = 101;
|
|
% Range for X and Y values [minX maxX minY maxY minZ maxZ]
|
|
XY = [0 4 -1 1];
|
|
% Boundary potential
|
|
|
|
% tolerance for ending iterations
|
|
tol = 1;
|
|
% limits for electric field plot
|
|
minXL = 0; maxXL = XY(2);
|
|
minYL = 0; maxYL = XY(4);
|
|
|
|
% boundary values for the potential
|
|
V0 = 100; V1 = -5; V2 = -10;
|
|
V = zeros(Ny,Nx);
|
|
V(:,1) = V0;
|
|
% V(:,end) = 0;
|
|
% V(1,:) = 0;
|
|
% V(end,:) = 0;
|
|
|
|
% linearly increasing potentials on boundaries
|
|
%enter INPUT details in setup
|
|
|
|
|
|
% SETUP ==================================================================
|
|
|
|
minX = XY(1); maxX = XY(2);
|
|
minY = XY(3); maxY = XY(4);
|
|
x = linspace(minX, maxX,Nx);
|
|
y = linspace(minY, maxY,Ny);
|
|
[xx, yy] = meshgrid(x,y);
|
|
hx = x(2) - x(1); hy = y(2) - y(1);
|
|
Kx = hx^2/(2*(hx^2+hy^2)); Ky = hy^2/(2*(hx^2+hy^2));
|
|
|
|
% V(:,1) = V0 .* cos(2*pi*y/(4*maxY));
|
|
|
|
% linearly increasing potentials on boundaries
|
|
m1 = 20/maxY; m2 = 20/maxX; m3 = 30/maxY; m4 = 10/maxX;
|
|
b1 = 0; b2 = 20; b3 = 10; b4 = 0;
|
|
|
|
V(:,1) = m1 .* y + b1;
|
|
V(:,end) = m3 .* y + b3;
|
|
V(end,:) = m2 .* x + b2;
|
|
V(1,:) = m4 .* x + b4;
|
|
|
|
% CALCULATIONS ===========================================================
|
|
|
|
% dSum difference in sum of squares / n number of iterations
|
|
dSum = 10; n = 0;
|
|
|
|
while dSum > tol
|
|
sum1 = sum(sum(V.^2));
|
|
|
|
for ny = 2: Ny-1
|
|
for nx = 2: Nx-1
|
|
V(ny,nx) = Ky * (V(ny,nx+1) + V(ny,nx-1)) + Kx * (V(ny+1,nx) + V(ny-1,nx));
|
|
end
|
|
% uncomment the next line of code if you need to calculate this column of values
|
|
% V(ny,Nx) = Ky * (V(ny,Nx) + V(ny,Nx-2)) + Kx * (V(ny+1,Nx) + V(ny-1,Nx));
|
|
end
|
|
|
|
sum2 = sum(sum(V.^2));
|
|
dSum = abs(sum2 - sum1);
|
|
n = n+1;
|
|
dSum
|
|
end
|
|
% electric field
|
|
[Exx, Eyy] = gradient(V,hx,hy);
|
|
Exx = -Exx; Eyy = -Eyy;
|
|
E = sqrt(Exx.^2 + Eyy.^2);
|
|
|
|
|
|
% exact solution for boundary conditions
|
|
% left boundary V0, other boundaries V = 0
|
|
% Vexact = (2*V0/pi) .* atan(sin(pi*y(51)/maxX)./sinh(pi*x/maxX));
|
|
% Vexact = (2*V0/pi) .* atan(1./sinh(pi*x/maxX));
|
|
|
|
|
|
% linearly increasing potentials on boundaries
|
|
m1 = 20/maxY; m2 = 20/maxX; m3 = 30/maxY; m4 = 10/maxX;
|
|
b1 = 0; b2 = 20; b3 = 10; b4 = 0;
|
|
|
|
V(:,1) = m1 .* y + b1;
|
|
V(:,end) = m3 .* y + b3;
|
|
V(end,:) = m2 .* x + b2;
|
|
V(1,:) = m4 .* x + b4;
|
|
|
|
% GRAPHICS ===============================================================
|
|
|
|
figure(1) % VECTOR FIELD V -----------------------------------------
|
|
set(gcf,'units','normalized','position',[0.02 0.52 0.3 0.32]);
|
|
surf(xx(1:5:end,1:5:end),yy(1:5:end,1:5:end),V(1:5:end,1:5:end));
|
|
xlabel('x [m]'); ylabel('y [m]'); zlabel('V [ V ]');
|
|
title('potential','fontweight','normal');
|
|
set(gca,'fontsize',14);
|
|
rotate3d
|
|
box on
|
|
axis tight
|
|
colorbar
|
|
view(55,49);
|
|
% set(gca,'ZTick',[0 5 10]);
|
|
|
|
%%
|
|
figure(2)
|
|
set(gcf,'units','normalized','position',[0.65 0.52 0.3 0.32]);
|
|
c = 0; yStep = zeros(1,6);
|
|
|
|
for n = 1:10:51
|
|
xP = xx(n,:); yP = V(n,:);
|
|
plot(xP,yP,'b','linewidth',1.2);
|
|
hold on
|
|
c = c+1;
|
|
yStep(c) = yy(n,1);
|
|
|
|
end
|
|
tt = num2str(yStep,2);
|
|
tm1 = 'y = ';
|
|
tm2 = tt;
|
|
tm = [tm1 tm2];
|
|
xlabel('x [m]'); ylabel('V [ V ]');
|
|
set(gca,'fontsize',16)
|
|
h_title = title(tm);
|
|
set(h_title,'fontsize',12,'FontWeight','normal');
|
|
|
|
|
|
% hold on
|
|
% plot(x,Vexact,'r');
|
|
|
|
%%
|
|
figure(3)
|
|
set(gcf,'units','normalized','position',[0.65 0.1 0.3 0.32]);
|
|
|
|
hold on
|
|
sx = x(5);
|
|
for sy = -1 : 0.1 : 1
|
|
h = streamline(xx,yy,Exx,Eyy,sx,sy);
|
|
set(h,'linewidth',2,'color',[1 0 1]);
|
|
end
|
|
|
|
|
|
index1 = 10 : 1: Nx; index2 = 10 : 1 : Ny;
|
|
p1 = xx(index1, index2); p2 = yy(index1, index2);
|
|
p3 = Exx(index1, index2); p4 = Eyy(index1, index2);
|
|
h = quiver(p1,p2,p3,p4);
|
|
set(h,'color',[0 0 1],'linewidth',2)
|
|
xlabel('x [m]'); ylabel('y [m]');
|
|
|
|
title('electric field','fontweight','normal');
|
|
hold on
|
|
axis equal
|
|
set(gca,'xLim',[minX,maxX]);
|
|
set(gca,'yLim',[minY,maxY]);
|
|
set(gca,'xTick',minX:1:maxX);
|
|
set(gca,'yTick',minY:1:maxY);
|
|
set(gca,'fontsize',14)
|
|
box on
|
|
|
|
%%
|
|
figure(4)
|
|
set(gcf,'units','normalized','position',[0.02 0.1 0.3 0.32]);
|
|
surf(xx,yy,E);
|
|
shading interp
|
|
xlabel('x [m]'); ylabel('y [m]'); zlabel('|E| [ V/m ]');
|
|
set(gca,'fontsize',14)
|
|
%set(gca,'zTick',[0 1000 2000]);
|
|
rotate3d
|
|
box on
|
|
axis tight
|
|
colorbar
|
|
view(43,54);
|
|
title('electric field','fontweight','normal');
|
|
|
|
figure(5)
|
|
set(gcf,'units','normalized','position',[0.33 0.52 0.3 0.32]);
|
|
contourf(xx,yy,V,16);
|
|
shading interp
|
|
xlabel('x [m]'); ylabel('y [m]');
|
|
title('potential','fontweight','normal');
|
|
set(gca,'fontsize',14)
|
|
box on
|
|
colorbar
|
|
axis equal
|
|
|
|
figure(6)
|
|
set(gcf,'units','normalized','position',[0.33 0.1 0.3 0.32]);
|
|
contourf(xx,yy,E,32);
|
|
shading interp
|
|
xlabel('x [m]'); ylabel('y [m]');
|
|
title('electric field','fontweight','normal');
|
|
set(gca,'fontsize',14)
|
|
box on
|
|
colorbar
|
|
axis equal
|
|
|
|
|
|
toc
|