3D poisson for electrostatics using solvepde

These are coordinates of atoms: [1.8150 3.6300 6.2895; 3.6300 1.8150 7.5000; 3.6300 1.8150 5.0789; -2.7225 4.7186 6.2895; 0.9075 6.5336 7.5000; 0.9075 6.5336 5.0789; 2.7225 4.7186 6.2895; 6.3525 6.5336 7.5000; 6.3525 6.5336 5.0789; 0 9.4373 6.2895; 3.6300 11.2523 7.5000; 3.6300 11.2523 5.0789; 5.4450 9.4373 6.2895; 9.0750 11.2523 7.5000; 9.0750 11.2523 5.0789; 2.7225 14.1559 6.2895];

There are 16 atoms. I get a row matrix of charge density(rho) from negf solver of size 1×16. x axis is the transport direction with bias voltage 0.5V. in z axis I apply electric field of 5V. In y axis electrostatic doesn’t change. I want dirichlet boundary condition in x axis. Neumann boundary in z axis. Periodic in y axis. All atoms with lowest & highest x values correspond to boundary in x direction. And atoms with lowest & highest z values correspond to boundary in z direction. Poisson equation to be solve is in form Del^2U(x,y,z)=-rho(x,y,z)/epsilon. Any guide on how to do it using solvepde?

Answer 1

[Karan Singh]

Hi Md Niloy Khan,

To solve the Poisson equation using the “solvepde” function in MATLAB, you need to define the geometry, specify the boundary conditions, and set up the PDE coefficients.

1. Define the geometry: Create a 2D geometry using the “geometryFromEdges” function. Define the boundaries corresponding to the x-axis and z-axis boundaries as separate edge sets.
2. Create a mesh: Generate a mesh using the “generateMesh” function. Specify the geometry and desired mesh size.
3. Define the PDE coefficients: Create a PDE model using the “createpde” function. Set the coefficients of the Poisson equation, including the diffusion coefficient (c), the source term (f), and the electric permittivity (epsilon).
4. Define the boundary conditions: Use the “applyBoundaryCondition” function to specify the boundary conditions. For the x-axis boundaries, use the Dirichlet boundary condition ('dirichlet') with a constant value. For the z-axis boundaries, use the Neumann boundary condition ('neumann') with a zero gradient. For the periodic boundary condition in the y-axis, use the periodic condition ('periodic').
5. Solve the PDE: Call the “solvepde” function, passing the PDE model, mesh, and boundary conditions as arguments. This will solve the Poisson equation and return the solution.

Here is an example code snippet:

% Define the coordinates of the atoms 
coords = [1.8150 3.6300 6.2895; 3.6300 1.8150 7.5000; 3.6300 1.8150 5.0789; -2.7225 4.7186 6.2895; 0.9075 6.5336 7.5000; 0.9075 6.5336 5.0789; 2.7225 4.7186 6.2895; 6.3525 6.5336 7.5000; 6.3525 6.5336 5.0789; 0 9.4373 6.2895; 3.6300 11.2523 7.5000; 3.6300 11.2523 5.0789; 5.4450 9.4373 6.2895; 9.0750 11.2523 7.5000; 9.0750 11.2523 5.0789; 2.7225 14.1559 6.2895]; 

% Create the geometry 
g = geometryFromEdges(edges); 

% Create the mesh 
mesh = generateMesh(g, 'Hmax', 0.1); 

% Create the PDE model 
model = createpde(); 

% Set the PDE coefficients 
c = 1; % Diffusion coefficient 
f = -rho; % Source term 
epsilon = 1; % Electric permittivity 

specifyCoefficients(model, 'm', 0, 'd', 0, 'c', c, 'a', 0, 'f', f, 'epsilon', epsilon); 

% Define the boundary conditions 
applyBoundaryCondition(model, 'dirichlet', 'Edge', xBoundaryEdges, 'u', 0.5); 
applyBoundaryCondition(model, 'neumann', 'Edge', zBoundaryEdges); 
applyBoundaryCondition(model, 'periodic', 'Edge', yBoundaryEdges); 

% Solve the PDE 
results = solvepde(model, mesh); 

% Extract the solution 
U = results.NodalSolution; 

This code provides a starting point for solving the Poisson equation using the “solvepde” function in MATLAB. You may need to adapt and modify it based on your specific requirements and problem setup.

Attached below are some documentation links that you may find helpful:

• Create 2-D geometry from decomposed geometry matrix - MATLAB geometryFromEdges (mathworks.com)
• Create triangular or tetrahedral mesh - MATLAB generateMesh (mathworks.com)
• Create model - MATLAB createpde (mathworks.com)
• Specify coefficients in a PDE model - MATLAB specifyCoefficients (mathworks.com)
• Add boundary condition to PDEModel container - MATLAB applyBoundaryCondition (mathworks.com)

Hope this helps!

[Duncan Carlsmith]

How does one define rho? Please provide a working example.