I am trying to solve a field-circuit coupled problem using the Gridap.jl package.
I want to simulate the current distribution inside a high-voltage cable, as well as the associated magnetic field.
This means I have two equations I would like to solve simultaneously:
A second order differential equation for the magnetic potential in terms of an imposed current density. This equation includes a term for eddy currents (time-harmonic assumption). \nabla\cdot(\nu \nabla A_z) - j\omega\sigma A_z = -J_z
A circuit equation which imposes a certain current I_{coil} on the cable. \iint_{coil} J_z d\Omega = I_{coil}
I have successfully solved this problem using a home-brew FEM implementation, which manually constructs the linear system with a loop over the elements (see this Jupyter notebook, and solution figures for 3 different frequencies). The results from this code match quite well with those from COMSOL Multiphysics.
I am wondering how to solve this problem using Gridap.jl? The first equation is no problem, as this can be solved by just following the first tutorial. However, I am unable to find a suitable starting point for including equation 2 in the documentation, and would greatly appreciate if someone could point me in the right direction.
I have been looking into alternatives, and it seems this can be done quite easily using constraints in Ferrite.jl (Constraints · Ferrite.jl). I have yet to try this, but it looks promising.
Does anyone know if something similar is possible in Gridap.jl?
Did you look at tutorial 18 (18 Topology optimization · Gridap tutorials) by chance? As part of the optimization, they use a radius to filter data (look at the “Topology Optimization” and “Filter and Threshold” sections).
Do you know of a way to rewrite your problem such that you could turn part of it into an optimization problem?
Thanks for the pointer! I’ve (very quickly) read through this tutorial and it seems that the section on the weak form might help in adding an additional equation. I have yet to try this, but will make a post if it succeeds.
we can then make things simpler by dividing the weak form to a base integral that contains the whole computation cell and an additional integral on the design domain.
As for turning this into an optimization problem, that seems like a waste of computation power: the additional constraint for current stated above should translate to a single extra equation in the linear system, with a marginal effect on the solving of the system.
Even though the tutorial 18 is labeled “Topology optimization”, I think with the section referred to above, there is no need to actually do any optimization.
If you want to use Ferrite, your code in the jupyter notebook should be quite easily translated to Ferrite with the help of the first tutorial in the docs. But If I understand you correctly, you also want to implement your second equation using periodic boundary constraints, instead of implementing it in a weak form? It is possible with AffineConstraints, but the tricky part might be to find which dofs to couple on the surfaces…
I agree that the first equation (for the magnetic potential) should be pretty easy to implement using Ferrite. But I’m not sure if I understand what you mean by implementing the second one with periodic boundary constraints.
I have no preference for the specific implementation, as long as I can do the following:
I have a physical problem: namely, to constrain the total current through a wire to a certain value (I_{coil}).
Mathematically, this translates to the second equation stated in my original post.
What do you suggest, implementation-wise, to attain this goal?
In my own implementation, I have added an additional equation and degree of freedom (additional row and column in the A matrix) for this “constraint”. This is, as far as I am aware, an accepted way of including electrical circuit equations and coupling to the degrees of freedom on the mesh.
I thought you wanted to use periodic boundary conditions because you linked to the it in the Ferrite docs, but perhaps I misunderstood.
Unfortunately I dont know the best method to implement the constraint in Eq2, but the way you have done it with a sort of Lagrange-multiplier methods is pretty standard I think.
After defining the linear and bilinear terms of the weak form, we can use the low-level API to derive the cell contributions and assemble the linear system. This linear system can then be freely modified (e.g., extended with additional rows and columns). The cell contributions to these extra vectors can also be calculated with the low-level API in quite a nice way.