matlab静电场有限元分析

有限元分析基本流程：几何创建->边界条件设置->物性参数设置->模型求解->结果后处理（可视化）。
问题描述：

matlab代码（非PDE TOOLBOX实现）：

%% Create an electromagnetic model for electrostatic analysis.
close all;
clear all;
emagmodel = createpde('electromagnetic','electrostatic');

%% Create the geometry
% Rectangle is code 3, 4 sides, followed by x-coordinates and then y-coordinates
R1 = [3,4,0,12,12,0,0,0,2,2]'; % air
R2 = [3,4,0,12,12,0,2,2,4,4]'; % dielectric slab
R3 = [3,4,0,12,12,0,4,4,8,8]'; % air
h=0.1;w=1;
R4 = [3,4,4.5,5.5,5.5,4.5,4,4,4.1,4.1]';
R5 = [3,4,6.5,7.5,7.5,6.5,4,4,4.1,4.1]';
% geometry description matrix
geom = [R1,R2,R3,R4,R5];
% Names for the two geometric objects
ns = (char('R1','R2','R3','R4','R5'))';
% Set formula
sf = '(R1+R2+R3)-R4-R5';
% Create geometry
gd = decsg(geom,sf,ns); % geometry description
geometryFromEdges(emagmodel,gd);

% View geometry
figure(1);
pdegplot(emagmodel,'EdgeLabels','on','SubdomainLabels','on')
axis equal

%% Solve problem
mu = 3.0; % relative permittivity of dielectric slab
emagmodel.VacuumPermittivity = 8.8541878128E-12;
electromagneticProperties(emagmodel,'RelativePermittivity',1.0, 'RelativePermeability',1.0,...,
                                    'Face',1);
electromagneticProperties(emagmodel,'RelativePermittivity',1.0, 'RelativePermeability',1.0,...
                                    'Face',2);
electromagneticProperties(emagmodel,'RelativePermittivity',mu, 'RelativePermeability',1.0, ...
                                    'Face',3);
% Apply the voltage boundary conditions on the edges of the square.
electromagneticBC(emagmodel,'Voltage',0.0,'Edge',[1 2 8 9 10 11 12 13]);
electromagneticBC(emagmodel,'Voltage',5.0,'Edge',[3 4 15 19]);
electromagneticBC(emagmodel,'Voltage',-5.0,'Edge',[5 6 17 20]);

% Generate the mesh.
mesh = generateMesh(emagmodel,'Hmax',0.15);
figure(2);pdeplot(mesh);

results = solve(emagmodel);
% plot electric field
figure(3);
pdeplot(emagmodel,'FlowData',[results.ElectricField.Ex ...
                              results.ElectricField.Ey])
axis equal
% plot electric potential
figure(4);
pdeplot(emagmodel,'XYData',results.ElectricPotential,'ZData',results.ElectricPotential,...
    'ColorMap','hot','Contour','on');

% 1.calculate the potential at points P(8,5),Q(4,4),R(6,1),S(6,3),T(2,3)
P_V = interpolateElectricPotential(results,8,5);
P_E = interpolateElectricField(results,8,5);
P_D = interpolateElectricFlux(results,8,5);
disp('Potential, electric field |E| and flux density |D| at point P(8,5) is: ');
disp(P_V);disp(sqrt(P_E.Ex^2+P_E.Ey^2));disp(sqrt(P_D.Dx^2+P_D.Dy^2));

Q_V = interpolateElectricPotential(results,4,4);
Q_E = interpolateElectricField(results,4,4);
Q_D = interpolateElectricFlux(results,4,4);
disp('Potential, electric field |E| and flux density |D| at point Q(4,4) is: ');
disp(Q_V);disp(sqrt(Q_E.Ex^2+Q_E.Ey^2));disp(sqrt(Q_D.Dx^2+Q_D.Dy^2));

R_V = interpolateElectricPotential(results,6,1);
R_E = interpolateElectricField(results,6,1);
R_D = interpolateElectricFlux(results,6,1);
disp('Potential, electric field |E| and flux density |D| at point R(6,1) is: ');
disp(R_V);disp(sqrt(R_E.Ex^2+R_E.Ey^2));disp(sqrt(R_D.Dx^2+R_D.Dy^2));

S_V = interpolateElectricPotential(results,6,3);
S_E = interpolateElectricField(results,6,3);
S_D = interpolateElectricFlux(results,6,3);
disp('Potential, electric field |E| and flux density |D| at point S(6,3) is: ');
disp(S_V);disp(sqrt(S_E.Ex^2+S_E.Ey^2));disp(sqrt(S_D.Dx^2+S_D.Dy^2));

T_V = interpolateElectricPotential(results,2,3);
T_E = interpolateElectricField(results,2,3);
T_D = interpolateElectricFlux(results,2,3);
disp('Potential, electric field |E| and flux density |D| at point T(2,3) is: ');
disp(T_V);disp(sqrt(T_E.Ex^2+T_E.Ey^2));disp(sqrt(T_D.Dx^2+T_D.Dy^2));

% 2.potential cotour plot
figure(5);
phi = results.ElectricPotential;
pdeplot(emagmodel,'XYData',phi,'Contour','on','ColorMap','hot');

%% 3. calculate cap
Dx = results.ElectricFluxDensity.Dx;
Dy = results.ElectricFluxDensity.Dy;

Q = sum(sqrt(Dx.*Dx+Dy.*Dy));
dU = 10; % V
cap = Q/dU; 

结果如下：