Vector Poisson Equation for a 3D Current Source

I thought given this is a separate issue, it would be worth a separate topic to make it easier to follow where I'm stuck.

In short, the aim of my model is to take a 3D current loop (initially with rectangular cross section - see my other post for a geometry reference) and model the resulting magnetic vector potential in 'all of space', where the containing box is assumed to be large enough that magnetic flux is completely confined within and Dirchlet zero conditions are appropriate for the outer faces. I also hope to calculate the related magnetic field over the same domain for further analysis.

For simplicity, the current density is assumed to have a constant vector magnitude throughout the coil, with vector components assigned to create a current flow around the 3D object. Something like this:



I'm fairly happy that this is possible within this toolbox, however after a significant number of hours effort I'm currently the limiting factor in making this work.

The simplified system of equation being studied is a vector Poisson equation for the magnetic vector potential due to a 3D current source. To implement this through the GUI I thought the inbuilt Poisson equation would be a good place to start to generate syntax giving (after removing the time derivative):

d_poi*(Ax_x + Ay_y + Az_z) = f_poi

From this, assuming mu = mu_o  = 4pi... throughout the domain:

d_poi = mu_o
Ax_x = d/dx Ax ect
f_poi = J.

This is where I come unstuck because I'm not sure how the div -D gradA is defined in this syntax to then convert that to
nabla A = - mu_o J;
d^2/dx^2 Ax = -mu_o * Jx ;
d^2/dy^2 Ay = -mu_o * Jy ;
d^2/dz^2 Az = -mu_o * Jz ;

I think a custom equation with 3 dependent variables is needed, but this scalar expression seems to have vector components too so I'm not actually certain. Also, changing the syntax/expression doesn't seem to alter the expression displayed in the GUI which is making this harder because I'm unsure how exactly the things I type are (or aren't) being interpreted.

I'm fairly certain this is a really simple task (though perhaps assuming the current flows maybe a little more difficult) but I'm struggling with the learning curve to setup something useful. Still working at this so will post an update if I make any progress - apologies for all the very basic questions and thanks for the help in advance.

Unfortunately, I can't quite follow and see the equation you want to define. If you can maybe write it down out paper (ideally all components if its a vector equation), scan it, and post a link I can likely give better suggestions.

So this is an magnetostatics/quasistatics problem in the presence of source currents without any permanent magnetic material, the full equation I'd like to model is:

full_eq.png,

As a set of PDEs then this is:
anisi-pdes.png

where A is the magnetic vector potential, mu = mu(r) & v = 1/mu, is a spatially varying permativity and J is a vector describing current density. Note in this form A can vary by any arbitrary scalar function without physically affecting the physical fields and from reading solving it can be very difficult. Overcoming this then..


Assuming mu(r) = mu_o throughout the entire domain and using the coulomb gauge del . A =0 the following vector identity holds:
vector_ident.PNG


This reduces the equation to:
vector_poisson.PNG

Which forms 3 (separable?) PDEs of the form:
1_comp.PNG,
where Az is the component of vector A along the Z axis, and Jz = the component of current density along the same axis.

(I would have followed your suggestion, but I can barely read my own handwriting xD)

So if I understand correctly you want to solve the system $\frac{\partial^2 Ax_i}{\partial x^2} + \frac{\partial^2 Ax_i}{\partial y^2} + \frac{\partial^2 Ax_i}{\partial z^2} = -\mu_0 Jx_i$ for $x_i = x,\ y,\ z$ in this case the corresponding FEATool equation form would read:

    A1x_x + A1y_y + A1z_z = -mu0*J1
    A2x_x + A2y_y + A2z_z = -mu0*J2
    A3x_x + A3y_y + A3z_z = -mu0*J3

where Ax = A1, Jx = J1, ... etc for clarity. You can set up this system by either using 3 custom equations, or by using a single custom equation and entering A1 A2 A3 in the "Dependent Variable Names" edit field.

custom_equation_3D_magnetostatics.fes

Note that even though you might use A1-3 to represent a vector field, these are technically three uncoupled scalar equations. (P.S. you can now use Latex/MathJax to enter equations.)

Brilliant thanks for the help, I can hopefully do some experimenting soon. Does the toolbox offer ways to impose gauge conditions explicitly or do I need to ensure I'm solving a numerically stable system of equations?

The (hopefully I know I must be getting annoying!) final question I have is with regards to the coefficients J1 J2 J3.

From trying to do FEM in other methods I have a function (attached with its dependent) which, for a 3D coil with rectangular cross section, can map the current density to nodes (well kind of, honestly the function isn't great because it only works if nodes are symmetrically distributed through the coil and the surface nodes at the corners should definitely be assigned differently. The aim is sqrt(J1^2 + J2^2 + J3^2) = Constant everywhere in the coil, with rotational flow in the XY plane like the initial figure on this post in case you need to test this to answer my question).

Currently it is of the form: [Jx,Jy,Jz] = J_coeff(x,y,z,N_wire,I_mag,mu_o).

Given this setup I have two simple questions:

1. Do I need to split this into 3 functions which only output one of J1, J2, J3?
2. In my geometry for the initial problem, I have two sub-domains: an internal coil of the type discussed above, surrounded by a 'universe' box which forms a typical black-box setup. Currently if I understand the GUI correctly, I set Jx as a constant which seems to apply to both sub-domains. The expectation is that this will only apply to the internal sub-domain, but the values in the outer-sub-domain call the internal as a 3D source. How do I approach ensuring this is achieved? I thought perhaps a different equation is needed for each sub-domain, but I am unsure.




J_coeff.mcoilcorner.m

ADSW1243 wrote

Currently it is of the form: [Jx,Jy,Jz] = J_coeff(x,y,z,N_wire,I_mag,mu_o).

1. Do I need to split this into 3 functions which only output one of J1, J2, J3?

Yes, but you can potentially use the "persistent" variables in a function "wrapper" to only calculate Jx-z once, for example:

function [ Ji ] = J_coeff_i( x, y, z, i_component )
  persistent Jx Jy Jz
  if( isempty(Jx) || ~isequal(size(x),size(Jx)) )
    [Jx,Jy,Jz] = J_coeff(x,y,z,N_wire,I_mag,mu_o);
    Jx = reshape(Jx,size(x));
    Jy = ...
  end
  switch i_component
    case 1
      Ji = Jx;
    ...
  end

You would then define J1-3 using the wrapper as "J1 := J_coeff_i(x,y,z,1)" etc ... If FEATool can find the function it will call it with the arguments. Note, that the output must have the same size as the input coordinates (for example using the reshape command).

ADSW1243 wrote

2. In my geometry for the initial problem, I have two sub-domains: an internal coil of the type discussed above, surrounded by a 'universe' box which forms a typical black-box setup. Currently if I understand the GUI correctly, I set Jx as a constant which seems to apply to both sub-domains. The expectation is that this will only apply to the internal sub-domain, but the values in the outer-sub-domain call the internal as a 3D source. How do I approach ensuring this is achieved? I thought perhaps a different equation is needed for each sub-domain, but I am unsure.

Assuming "J1" is defined in your equation, and you define an expression "J1". If the expression J1 is a scalar it will be valid for all subdomains, but if you define it as a space separated list, each entry will be used for the corresponding subdomains, for example "J1 := '1 0'" would set J1 to one in subdomain 1 and zero in subdomain 2.

So combining the two parts gives (for one component): "J1 := 'J_coeff_i(x,y,z,1) 0'", where the function is called only for subdomain 1?