How to set the boundary conditions of 3D Poisson Equation
15 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Anthony Owusu-Mensah
le 2 Mar 2020
Commenté : David Martín Oliva Zavala
le 12 Nov 2020
I am trying to compute the electric potential at point (x,y,z) by solving the 3D Poisson equation below using finite difference method.
I have the charge densities at various positions of (x,y,z).
Below are the boundary condtions;
- Dirichlet boundary condition are applied at the top and bottom of the planes of the rectangular grid.
- Electric potential is to be incorporated by setting and , where h is the height of the simulation box.
- Neumann boundary conditions are also enforced at the remaining box interfaces by setting at faces with constant x , at faces with constant y, and at faces with constant z.
I did the implementation of the boundary conditions with this code and I would want to ascertain if
the implementation of the boundary condtions is right.
x1 = linspace(0,10,20);
y1 = linspace(0,10,20);
z1 = linspace(0,10,20);
V = zeros(length(x1),length(y1),length(z1));
%Dirichlet Boundary Conditions
%Top plane
V(:,end,end) = 0;
V(end,:,end) = 0;
V(end,end,:) = 0;
% Bottom plane
V(:,1,1) = 0;
V(1,:,1) = 0;
V(1,1,:) = 0;
%Incoporated electric potential
V(:,:,1) = 0;
V(:,:,end) = -40*z1(end);
%Neumann Boundary Condition
i = 2:length(x1)-1;
j = 2:length(y1)-1;
k = 2:length(z1)-1;
V(i+1,j,k) = V(i-1,j,k);
V(i,j+1,k) = V(i,j-1,k);
V(i,j,k+1) = V(i,j,k-1);
1 commentaire
David Martín Oliva Zavala
le 12 Nov 2020
Hey Anthony, have you finished solving poisson eqn in 3D? I´m just starting this project and I may need some help in the future jeje:)
Réponse acceptée
Dinesh Yadav
le 5 Mar 2020
Modifié(e) : Dinesh Yadav
le 5 Mar 2020
Hi Anthony I think you have wrongly defined the top and bottom planes.A plane is a 2-D sheet structure however the line of code
V(:,end,end);
represent a single row vector, similar to a line and not a plane. As V is a 20x20x20 3-D matrix, lets assume the first plane facing us is top plane i.e (assuming x direction increases no of columns and y direction increases no of rows, i.e origin is at (20,1,1))
V(:,:,1); % top plane face with constant z
V(:,:,end); % bottom plane face with constant z
V(:,1,:); % face with constant x
V(:,end,:) % face with constant x
V(1,:,:); % face with constant y
V(end,:,:) % face with constant y
dx = 2:length(x1)-1;
dy = 2:length(y1)-1;
dz = 2:length(z1)-1;
V(dy,dx+1,dz) = V(dy,dx-1,dz); % dv/dx = 0
V(dy+1,dx,dz) = V(dy-1,dx,dz); % dv/dy = 0
V(dy,dx,dz+1) = V(dy,dx,dz-1); % dv/dz = 0
3 commentaires
Dinesh Yadav
le 5 Mar 2020
Yes, replace i with dy j with dx and k with dz. Just to make it more readable. I have corrected the above code.
Anthony Owusu-Mensah
le 5 Mar 2020
Modifié(e) : Anthony Owusu-Mensah
le 5 Mar 2020
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Graphics Object Programming dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!