Test Case : Comsol - Cylinder

1. Introduction

This is the test case of High-Temperature Superconductors using the H Formulation on a Bulk Cylinder geometry surrounded by air in axisymmetric coordinates.

2. Run the Calculation

To run this case, compute the Study 1 solver on Comsol Multyphysics on the file CylindreBulkFormH_bis.mph.

3. Data Files

The case data file is available on Github here.

4. Equation

The H Formulation in axysimmetric coordinates is :

H Formulation in axisymmetric coordinates
\[\text{(H Axi)} \left\{ \begin{matrix} \nabla\times\left(\rho\nabla\times H\right) + \frac{\partial (\mu H)}{\partial t} &=& 0 & \text{on } \Omega^{axis} & \text{(H Axi)}\\ \mathbf{H} \times \mathbf{n} &=& 0 &\text{ on } \Gamma_D^{axis} & \text{(D Axi)} \\ \left( \nabla \times \mathbf{H} \right) \times \mathbf{n} &=& 0 &\text{ on } \Gamma_N^{axis} & \text{(N Axi)} \end{matrix} \right.\]

With :

  • \(H\) : magnetic field

  • \(\rho\) : electric resistivity \(\Omega \cdot m\) described by the e-j power law : \(\rho=\frac{E_c}{J_c}\left(\frac{\mid\mid J \mid\mid}{J_c}\right)^{(n)}=\frac{E_c}{J_c}\left(\frac{\mid\mid \nabla\times \textbf{H} \mid\mid}{J_c}\right)^{(n)}\)

  • \(\mu\) : electric permeability \(kg/A^2/S^2\)

5. Geometry

The geometry is a set of stacked tapes in axisymmetric coordinates \((r,z)\), surrounded by air. The geometry is created on Comsol, but in order to show the different physical names, this is a recreation on gmsh :

cylinder
Geometry in Axisymmetrical cut

The geometrical domains are :

  • Cylinder : the cylinder, composed by a conductor

  • Air & Spherical Shell : the air surrounding Conductor

    • Symmetry Line : Air 's bound, correspond to \(Oz\) axis (\(\{(z,r), \, z=0 \}\))

    • Exterior Boundary : the rest of the Air 's bound

Symbol

Description

value

unit

\(R\)

radius of cylinder

\(0.00125\)

m

\(H_{cylinder}\)

height of cylinder

\(0.001\)

m

\(R_{inf}\)

radius of infty border

\(0.1\)

m

6. Parameters

  • On Cylinder :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

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

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

\(t_f\)

final time

\(15\)

\(s\)

\(b_{max}\)

maximum applied field

\(1\)

\(T\)

\(rate\)

rate of the applied field raise

\(\frac{3}{tf}b_{max}\)

\(T/s\)

\(hsVal\)

time ramp for the applied field

\(\begin{cases}rate*t &\quad\text{if }t<\frac{t_f}{3}\\b_{max} &\quad\text{if }t<\frac{2t_f}{3}\\b_{max} - (t-\frac{2t_f}{3})*rate &\quad\text{if }t>\frac{2t_f}{3}\end{cases}\)

\(T\)

\(J_c\)

critical current density

\(3.10^8\)

\(A/m^2\)

\(E_c\)

threshold electric field

\(10^{-4}\)

\(V/m\)

\(n\)

material dependent exponant

\(20\)

\(\rho\)

electrical resistivity (described by the \(e-j\) power law)

\(\frac{E_c}{J_c}\left(\frac{\mid\mid J \mid\mid}{J_c}\right)^{(n)}\)

\(\Omega\cdot m\)

  • On Air :

Symbol

Description

Value

Unit

\(\mu=\mu_0\)

magnetic permeability of vacuum

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

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

\(\rho_{air}\)

electrical resistivity of the air

\(100\)

\(\Omega\cdot m\)

On MPH file, the parameters are written in the parameters category, the materials and as ad piecewise function for the applied magnetic field :

Table 1. Parameters on MPH file
comsolaxiparam
parameters
comsolaxiparam2
parameters
comsolaxiparam3
parameters

7. Boundary Conditions

For the Dirichlet boundary conditions, we want to impose the applied magnetic field :

  • On Exterior Boundary : \(H_z = hsVal(t)/\mu\)

On MPH file, the boundary conditions are written :

Table 2. Boundary conditions on MPH file
comsolaxidir
Dirichlet

8. Weak Formulation

Weak formulation of H Formulation
\[\text{(Weak H)} \begin{cases} \int_{\Omega}{ \rho \nabla\times \textbf{H} \cdot \nabla\times\Phi} = \int_{\Omega}{ \partial_t(\mu \textbf{H}) \cdot \Phi}\\ \forall \Phi \in H^{curl}(\Omega)\end{cases}\]

9. Faraday’s law

For the H Formulation, the physic Magnetic Field Formulation (mfh) is used on Comsol. It allows to use the Faraday’s Law, which defined the H formulation :

comsolaxifara
Ampère’s Law

10. Numeric Parameters

  • Time

    • Initial Time : \(0s\)

    • Final Time : \(15s\)

    • Time Step : \(0.1s\)

  • Mesh

pyMPH mesh
Mesh

11. Results

The results that we obtain with this formulation with Feelpp are compared to the results of the article Finite-Element Formulations for Systems With high-temperature Superconductors where the solver getDP is used.

The time evolution of the applied field is :

applied field
Time evolution of the external applied field

With \(t_1=5s\), \(t_2=10s\) and \(t_3=15s\)

11.1. Electric current density

The electric current density \(J\) is defined by :

\[ J=\nabla\times\textbf(H)\]
Electric current density \(j_\theta (A/m^2)\)

We compare the current density profiles with Comsol and getDP on the \(O_r\) axis, at the mid-height of the cylinder, at time \(t_3\) for a maximum applied field of 1 T and \(n=20\), for a mesh of 30199 nodes.

L2 Relative Error Norm : \(10.72 \%\)

11.2. Magnetic flux density

The magnetic flux density \(B\) is defined by:

\[ B=\mu H\]
r_component of the magnetic flux density \(B_r (T)\)
z_component of the magnetic flux density \(B_z (T)\)

We compare the distribution of the z-component of the magnetic flux density 2mm above the cylinder at the instants \(t_1\), \(t_2\) and \(t_3\) with Comsol and getDP.

t1 \(=5s\)

L2 Relative Error Norm : \(1.15 \%\)

t2 \(=10s\)

L2 Relative Error Norm : \(5.59 \%\)

t3 \(=15s\)

L2 Relative Error Norm : \(2.99 \%\)

12. References

  • Finite-Element Formulation for Systems with High-Temperature Superconductors, Julien Dular, Christophe Gauzaine, Benoît Vanderheyden, IEEE Transactions on Applied Superconductivity VOL. 30 NO. 3, April 2020, PDF