how to apply subdomain boundary conditions

Hello,

I am trying to do magnetostatic analysis of 2d coaxial cable. I created the geometry here,

%%%%%

emagmodel = createpde(electromagnetic="magnetostatic");

R1 = [3,4,-10,10,10,-10,-10,-10,10,10]';

C1 = [1,0,0,1]';

C1 = [C1;zeros(length(R1) - length(C1),1)];

gm = [R1,C1];

sf = '(R1+C1)';

ns = char('R1','C1');

ns = ns';

g = decsg(gm,sf,ns);

geometryFromEdges(emagmodel,g);

pdegplot(emagmodel,EdgeLabels="on",FaceLabels="on")

axis equal

%%%%%

I selected arbitrary outer radius for amperian loop.

I need to apply boundary conditions on inner boundaries. In order to do that I need to subtract C1 from R1 but if I do it I cannot specify current density inside of the C1 since matlab doesn't assign any face on C1. I am open to any suggestions.

Answer 1

[Divyam]

Hi @mete ustundag,

The code uses the set formula "sf" to decompose the geometry and hence results in a figure with a single face. If you wish to add current density to the face C1, then you can decompose the geometry without specifying the set formula and add the respective properties to C1.

Since this does not align with your required workflow, the alternative solution here is to define a function handle that adds the "Current Density" to the location of the circle C1 when the solver is called. Here is a sample code

emagmodel = createpde(electromagnetic="magnetostatic");

emagmodel.VacuumPermeability = 1.2566370614E-6;

% Define the outer rectangle
R1 = [3,4,-10,10,10,-10,-10,-10,10,10]';

% Define the inner circle
C1 = [1,0,0,1,0,0,0,0,0,0]';

% Combine the geometries
gm = [R1,C1];
sf = 'R1-C1';
ns = char('R1','C1');
ns = ns';

% Create the geometry
g = decsg(gm,sf,ns);
geometryFromEdges(emagmodel,g);
pdegplot(emagmodel,EdgeLabels="on",FaceLabels="on");
axis equal

% Adding material properties
electromagneticProperties(emagmodel, 'Face', 1, 'RelativePermeability', 1);

% Adding current density to the circle
sc1 = electromagneticSource(emagmodel, "Face", 1, "CurrentDensity", @myfunCurrent);

% Example boundary condition
electromagneticBC(emagmodel, 'Edge', 1:4, 'MagneticPotential', 0); 

% Generate mesh
generateMesh(emagmodel);

% Solve the model
result = solve(emagmodel);

% Visualize the results
figure;
pdeplot(emagmodel, 'XYData', result.MagneticPotential, 'Contour', 'on');
title('Magnetic Potential');


function J = myfunCurrent(location, state)
    % Alter this based on the radius of the circle
    radius = 1;
    
    % Extract the x and y coordinates of the location
    x = location.x;
    y = location.y;
    % Initialize the current density array
    J = zeros(size(x));
    
    % Determine which points are inside the circle
    insideCircle = (x).^2 + (y).^2 <= radius^2;
    
    % Assign a current density value to points inside the circle
    % Example: Set a constant current density 
    J(insideCircle) = 1;
end

For more information regarding the use of function handles to define properties, refer to this documentation: https://www.mathworks.com/help/pde/ug/pde.electromagneticmodel.electromagneticbc.html#mw_6509b973-160a-4ae5-a9cf-9931c778f618:~:text=use%20function%20handles