Test of MQS + gauge

1. Introduction

In this page, I present the test case of Maxwell Quasi Static Problem with the A-V Formulation with gauge condition for a geometry of torus surrounded by air.

2. Run the calculation

The command line to run this case is :

    mpirun -np 32 feelpp_toolbox_coefficientformpdes --config-file=mqs.cfg --cfpdes.gmsh.hsize=5e-3

3. Data Files

The case data files are available in Github here :

CFG file - Edit the file
JSON file - Edit the file
GEO file - Edit the file

4. Equation

Thus \(\Omega\) the domain, comprising the conductor (or superconductor) domain \(\Omega_c\) and non conducting materials \(\Omega_n\) (\(\mathbf{J} = 0\)) like air. Let \(\Gamma = \partial \Omega\) the bound of \(\Omega\), \(\Gamma_c = \partial \Omega_c\) the bound of \(\Omega_c\), \(\Gamma_D\), \(\Gamma_{Dc}\) the bound with Dirichlet boundary condition and \(\Gamma_N\), \(\Gamma_{Nc}\) the bound with Neumann boundary condition, such that \(\Gamma = \Gamma_D \cup \Gamma_N\) and \(\Gamma_c = \Gamma_{Dc} \cup \Gamma_{Nc}\).

We introduce :

Magnetic potential field \(\mathbf{A}\) : the magnetic field writes \(\textbf{B} = \nabla \times \textbf{A}\)
Electric potential scalar : \(\nabla V = - \textbf{E} - \frac{d \textbf{A}}{d t}\)

We have the conditions :

\(\nabla \cdot \mathbf{A} = 0\)
\(\mu\) is independant of space

We want to resolve the electromagnetism problem ( with \(\mathbf{A}\) and \(V\) the unknows) :

A-V Formulation with gauge condition

\[\text{(AV+gauge)} \left\{ \begin{matrix} \sigma \frac{d \textbf{A}}{d t} + \frac{1}{\mu} \Delta \textbf{A} &=& - \sigma \nabla V & \text{ on } \Omega & \text{(AV+gauge-1)} \\ \nabla \cdot \left( \sigma \nabla V \right) &=& - \nabla \cdot \left( \sigma \frac{d \mathbf{A}}{d t} \right) & \text{ on } \Omega_c & \text{(AV'-2)} \\ \mathbf{A} \times \mathbf{n} &=& 0 & \text{ on } \Gamma_D & \text{(DA)} \\ \left( \nabla \times \mathbf{A} \right) \times \mathbf{n} &=& 0 & \text{ on } \Gamma_N & \text{(NA)} \\ V &=& V_d & \text{ on } \Gamma_{Dc} & \text{(DV)} \\ \frac{\partial V}{\partial \mathbf{n}} &=& 0 & \text{ on } \Gamma_{Nc} & \text{(NV)} \end{matrix} \right.\]

This equation is obtain from the A-V Formulation with the gauge conditions. The detail calculus are in page MQS + gauge condition.

5. Geometry

The geometry is a quart of torus of conductor surrounded by air.

Geometry

The geometrical domains are :

Conductor (\(\Omega_c\)) : the torus, it is composed of conductive materials
- V0 : enter of electrical potential
- V1 : exit of electrical potential

(V0 \(\cup\) V1 \( = \Gamma_{Dc}\))

air (\(\Omega/\Omega_c\)) : the air surround Conductor
- OXOZ : \(OxOz\) plan
- OYOZ : \(OyOz\) plan
- Infty : the rest of bound of Air

Symbol

Description

value

unit

\(r_{int}\)

interior radius of torus

\(75e-3\)

\(r_{ext}\)

exterior radius of torus

\(100.2e-3\)

\(z_1\)

half-height of torus

\(25e-3\)

\(r_{infty}\)

radius of infty border

\(5*r_{ext}\)

6. Initial/Boundary Conditions

We impose the Dirichlet boundary conditions :

For \(V\) equation :
- On V0 : \(V = 0\)
- On V1 : \(V = \frac{1}{4} \left\{ \begin{matrix} t \text{ on } 0s \leq t < 1s \\ 1 \text{ on } 1s \leq t < 20s \\ t \text{ on } 20s \leq t \end{matrix} \right.\)
For \(\mathbf{A}\) :
- On OXOZ and V0 : \(A_x = A_z = 0\), we want \(\mathbf{A}\) orthogonal to OXOZ and V0
- On OYOZ and V1 : \(A_y = A_z = 0\), we want \(\mathbf{A}\) orthogonal to OYOZ and V1
- Infty : We approximate the problem, Infty is the physical infty thus \(\mathbf{B}=0\) at Infty thus \(\mathbf{A} = 0\).

We initialize at \(t=0s\) :

\(V = 0\)
\(\mathbf{A} = 0\)

On JSON file, the boundary conditions are writed :

Boundary conditions on JSON file

    "BoundaryConditions":
    {
        "magnetic":
        {
            "Dirichlet":
            {
                "Infty":
                {
                    "expr":"{0,0,0}"
                }
            },
            "Dirichlet_x":
            {
                "mydir_y":
                {
                    "markers":["V0","OXOZ"],
                    "expr":0
                }
            },
            "Dirichlet_y":
            {
                "mydir_y":
                {
                    "markers":["V1","OYOZ"],
                    "expr":0
                }
            },
            "Dirichlet_z":
            {
                "mydir_z":
                {
                    "markers":["V0","OXOZ","V1","OYOZ"],
                    "expr":0
                }
            }
        },
        "electric":
        {
            "Dirichlet":
            {
                "V0":
                {
                    "expr":"V0:V0"
                },
                "V1":
                {
                    "expr":"V1:V1"
                }
            }
        }
    }

For initial condition, we write nothing in JSON, by default the application put zeros on initialization.

7. Weak Formulation

The weak formulation of A-V Formulation + gauge is :

\[\text{(AV+gauge Weak)} \left\{ \begin{eqnarray*} \int_{\Omega}{ \sigma \frac{d \mathbf{A}}{d t} \mathbf{ϕ} } + \int_{\Omega}{ \frac{1}{\mu} \nabla \mathbf{A} \cdot \nabla \mathbf{ϕ} } &=& - \int_{\Omega_c}{\sigma \nabla V \, \mathbf{ϕ}} \\ \int_{\Omega_c}{ \sigma \nabla V \cdot \nabla \psi } &=& - \int_{\Omega_c}{\sigma \frac{d \mathbf{A}}{d t} \cdot \nabla \psi} \end{eqnarray*} \right.\]

8. Parameters

The parameters of problem are :

On Conductor :

Symbol

Description

Value

Unit

\(V0\)

scalar electrical potential on V0

\(0\)

\(Volt\)

\(V1\)

scalar electrical potential on V1

\(\frac{1}{4} \left\{ \begin{matrix} t \ \text{on} \ 0s \leq t < 1s \\ 1 \ \text{on} \ 1s \leq t < 20s \\ 0 \ \text{on} \ 20s \leq t \end{matrix} \right.\)

\(Volt\)

\(\sigma\)

electrical conductivity

\(58e6\)

\(S/m\)

\(\mu=\mu_0\)

magnetic permeability of vacuum

\(4\pi.10^{-7}\)

\(kg \, m / A^2 / S^2\)

On Air :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

\(4\pi.10^{-7}\)

\(kg \, m / A^2 / S^2\)

On JSON file, the parmeters are writed :

Parameters on JSON file

    "Parameters":
    {
        "V0":0,
        "V1":"1/4*(t/(0.1*10)*(t<(0.1*10))+((t<(0.5*40))+0*(t>(0.5*40)))*(t>(0.1*10))):t"
    }

9. Coefficient Form PDEs

We use the application Coefficient Form PDEs. The coefficient associate to Weak Formulation are :

In the json file, we note magnetic the equation (AV'-1) and electric the equation (AV'-2).

On Conductor :

Coefficient

Description

Expression

magnetic_d

damping or mass coefficient

\(\sigma\)

magnetic_c

diffusion coefficient

\(\frac{1}{\mu}\)

magnetic_f

source term

\(- \sigma \nabla V \)

electric_c

damping or mass coefficient

\(\sigma\)

electric_gamma

source term

\(\sigma \, \frac{d \mathbf{A}}{d t}\)

On Air (the lectrical potential isn’t computed on Air) :

Coefficient

Description

Expression

magnetic_c

diffusion coefficient

\(\frac{1}{\mu}\)

On JSON file, the coefficients are writed :

CFPDEs coefficients on JSON file

    "Materials":
    {
        "Conductor":
        {
            "sigma":58e6,       // S.m-1
            "mu":"4*pi*1e-7",   // kg.m/A2/S2

            "magnetic_d":"sigma:sigma",
            "magnetic_c":"1/mu:mu",
            "magnetic_f":"{-sigma*electric_grad_V_0,-sigma*electric_grad_V_1,-sigma*electric_grad_V_2}:sigma:electric_grad_V_0:electric_grad_V_1:electric_grad_V_2",

            "electric_c":"sigma:sigma",
            "electric_gamma":"{sigma*magnetic_dA_dt_0,sigma*magnetic_dA_dt_1,sigma*magnetic_dA_dt_2}:sigma:magnetic_dA_dt_0:magnetic_dA_dt_1:magnetic_dA_dt_2"
        },
        "Air":
        {
            "physics":"magnetic",
            "mu":"4*pi*1e-7",   // kg.m/A2/S2
            "magnetic_c":"1/mu:mu"
        }
    }

10. Numeric Parameters

This section show the parameters used to compute the simulation.

Size of mesh :
- On Conductor : \(0.5 \, mm\)
- On Infty : \(100 \, mm\)
Time Parameters :
- Time step : \(0.1 \, s\)
- Initial Time : \(0 \, s\)
- Final Time : \(22 \, s\)
Element type : \(P1\)
Solver : automatic
Number of CPU core : \(32\)

Mesh of Geometry

11. Result

11.1. Electric Potential Scalar \(V\)

TODO : plot curves of \(V\) in function of \(\theta\)

11.2. Potential Magnetic Field \(\mathbf{A}\)

I plot the the components of \(\mathbf{A}\) at \(t=20s\) in stationnary state :

Components of \(\mathbf{A} \, (Tm)\) on Or axis at \(t=20s\)

Magnetic Potential Field \(A_x (T.m)\)

Magnetic Potential Field \(A_y (T.m)\)

As predicted, the z component of potential magentic field \(A_z\) is equal to zeros.

11.3. Magnetic Field \(\mathbf{B}\)

I plot the the components of \(\mathbf{B}\) on a radius of tore at \(t=20s\) in stationnary state :

Components of \(\mathbf{B} \, (T)\) on Or axis at \(t=20s\)

Horizontal Magnetic Field \(B_z(T)\)

12. Validation

12.1. Exact solution

The analytical solve of potential magnetic field and magnetic field (the cyan curve on plots) is computed by python3’s module MagnetTools.MagnetTools. The module based on paper : spire.

The analytical solve of intensity (and magnetic field) is computed on python3’s script Edit the file.

12.2. Convergence Test

I plot the result of z component of magnetic field \(B_z\) on a radius of tore at \(t=20s\) in stationnary state with exact solution computed with python3’s module MagnetTools.MagnetTools to "validate" the solution :

\(B_z (T)\) on Or axis at \(t=20s\)

The two curves overlap.

I plot the result of z component of magnetic field \(B_z\) on origin point \(P_0=(0,0,0)\) across time with numerical solution computed with getdp software to "validate" the solution :

\(B_z (T)\) on \(P_0=(0,0,0)\) across the time

The two curves overlap.

We do the same for \(A_{\theta}\) (the orthoradial component in cylindrical coordinates) :

\(A_{\theta} (Tm)\) on Or axis at \(t=20s\)

It exist a little difference near of Infty. To obtain better precision, we can increase the radius of Infty to have approximatively 0.

With the last comparisons, we have trust of our method.

13. References

Calcul du champ magnétique pour les géométries axisymétriques simples, Christophe Trophime, 2002, unpublished, Download the PDF