Skip to content
Snippets Groups Projects
Commit 3e5463fd authored by Albert Piwonski's avatar Albert Piwonski
Browse files

getDP solver file for helicoidal symmetric power cable model.

parent ac300fe1
No related branches found
No related tags found
1 merge request!5Albertpiwonski master patch 12427
Hide whitespace changes
Inline Side-by-side
HelicoidalSymmetricCable/helicoidal_power_cable.pro 0 → 100644
+ 707
− 0
View file @ 3e5463fd
Printf("GetDP solver code for helicoidally symmetric power cables");
/*
Documentation of the model:
original IEEE version: https://ieeexplore.ieee.org/document/10006756
arXiv version: https://arxiv.org/abs/2301.03370
The 10 basic GetDP steps are:
01. Group → Generate interface to mesh-file (.msh) via number of physical groups,
s.t. e.g. GetDP routines like FunctionSpace or Formulation can use those regions.
02. Function → Definition of "magnetoquasistatic" material mappings σ (conductivity) and
μ (permeablility) for different regions in the model.
03. Jacobian → Jacobian methods can be referred to in Formulation and PostProcessing objects
to be used in the computation of integral terms and for changes of coordinates.
04. Integration → Integration methods can be referred to in Formulation and PostProcessing
objects to be used in the computation of integral terms ... with options.
05. Constraint → Constraints can be referred to in FunctionSpace objects to be used for
boundary conditions ... (HERE): Fixing the cohomology conditions of the BVP!
06. FunctionSpace → is characterized by the type of its interpolated fields,
one or several basis functions and optional constraints (in space and time).
→ (HERE): V_uvw(Ω) (where magnetic field H_uvw is searched, i.e., H_uvw ∈ V_uvw(Ω))
= V_uv(Ω) ⊕ V_w(Ω)
= V_uv(Ω_i) ⊕ V_uv(Ω_c) ⊕ V_w(Ω_i) ⊕ V_w(Ω_c)
= [grad(nodal fcns.)(Ω_i)_uv & ℋ¹(Ω_i)_uv & edge_fcns.(Ω_c)_uv]
⊕ [const.(Ω_i)_w & perpendicular_edge_fcns.(Ω_c)_w]
07. Formulation → The Formulation tool permits to deal with volume, surface and line
integrals with many kinds of densities to integrate, written in a form that
is similar to their symbolic expressions ...
08. Resolution → "∂ₜ → jω" & Solving systems of equations.
09. PostProcessing → The PostProcessing object is based on the quantities defined in a
Formulation and permits the construction (thanks to the expression syntax)
of any useful piecewise defined quantity of interest.
10. PostOperation → The PostOperation is the bridge between results obtained with GetDP
and the external world.
*/
///////////////////////////////////////////////
//////////// step 01. Group ///////////////////
///////////////////////////////////////////////
Group {
Omega_c = Region[{3}]; // (whole) conducting domain, i.e. contains its boundary
Omega_i = Region[{5}]; // (whole) insulating domain, i.e. contains its boundary
Omega = Region[{Omega_c, Omega_i}]; // (whole) computational domain, i.e. contains its boundary
BdOmega_c = Region[{4}]; // boundary of conducting domain (consists of several bean-shaped outlines of the conductors)
BdOmega_i = Region[{6}]; // boundary of insulating domain (consists of several bean-shaped outlines of the conductors & an outer circle)
BdOmega = Region[{2}]; // boundary of computational domain (is an outer circle)
// NOTE: here comes now N thick cuts for N cohomology basis functions for N conductors:
Cut1TO = Region[{33}]; // thick cut for the 1. cohomology basis function
Cut2TO = Region[{34}]; // thick cut for the 2. cohomology basis function
Cut3TO = Region[{35}]; // thick cut for the 3. cohomology basis function
Cut4TO = Region[{36}]; // thick cut for the 4. cohomology basis function
Cut5TO = Region[{37}]; // thick cut for the 5. cohomology basis function
Cut6TO = Region[{38}]; // thick cut for the 6. cohomology basis function
Cut7TO = Region[{39}]; // thick cut for the 7. cohomology basis function
Cut8TO = Region[{40}]; // thick cut for the 8. cohomology basis function
Cut9TO = Region[{41}]; // thick cut for the 9. cohomology basis function
Cut10TO = Region[{42}]; // thick cut for the 10. cohomology basis function
Cut11TO = Region[{43}]; // thick cut for the 11. cohomology basis function
Cut12TO = Region[{44}]; // thick cut for the 12. cohomology basis function
Cut13TO = Region[{45}]; // thick cut for the 13. cohomology basis function
} // end of Group block
///////////////////////////////////////////////
//////////// step 02. Function ////////////////
///////////////////////////////////////////////
Function {
// Definition of "magnetoquasistatic" material mappings σ (conductivity) and μ (permeablility) in the cartesian chart:
// values:
mu_0 = 4*Pi*1e-7; // vacuum permeablility μ₀ [VsA⁻¹m⁻¹]
sigma_Cu = 58.1e6; // electric conductivity of copper σ_Cu [AV⁻¹m⁻¹]
// permeability μ [VsA⁻¹m⁻¹]:
mu[Omega_i] = mu_0; // ∀ x ∈ Ω_i : μ(x) = μ₀
mu[Omega_c] = mu_0; // ∀ x ∈ Ω_c : μ(x) = μ₀
mu[BdOmega_c] = mu_0; // ∀ x ∈ ∂Ω_c: μ(x) = μ₀
// conductivity σ [AV⁻¹m⁻¹] and resistivity ρ [VmA¯¹]:
sigma[Omega_i] = 0.; // ∀ x ∈ Ω_i: σ(x) = 0
sigma[Omega_c] = sigma_Cu; // ∀ x ∈ Ω_c: σ(x) = σ_Cu
rho[Omega_c] = 1./sigma_Cu; // ∀ x ∈ Ω_c: ρ(x) = 1/σ_Cu
// excitation parameters:
Itot[] = Sqrt[2]/13; // amplitude of the applied AC current through each of the 13 conductors
Freq = 50; // frequency of the applied AC current
// parameter used in the helicoidal transformation (coordinate transformation between (x, y, z) and (u, v, w) coords):
alpha = 1.0;
beta = 0.03183098861837907;
xyz_to_uvw[] = Vector[(+X[]*Cos[alpha/beta*Z[]]+Y[]*Sin[alpha/beta*Z[]]), (-X[]*Sin[alpha/beta*Z[]]+Y[]*Cos[alpha/beta*Z[]]), Z[]]; // coordinate transformation φ: p_xyz → p_uvw
uvw_to_xyz[] = Vector[(+X[]*Cos[alpha/beta*Z[]]-Y[]*Sin[alpha/beta*Z[]]), (+X[]*Sin[alpha/beta*Z[]]+Y[]*Cos[alpha/beta*Z[]]), Z[]]; // coordinate transformation φ⁻¹: p_uvw → p_xyz
// transformation matrices:
// full J:
// cos(w*α/β) -sin(w*α/β) -α*(u*sin(w*α/β) + v*cos(w*α/β))/β
// sin(w*α/β) cos(w*α/β) α*(u*cos(w*α/β) - v*sin(w*α/β))/β
// 0 0 1
J[] = Tensor[1, 0, -alpha/beta*Y[], 0, 1, alpha/beta*X[], 0, 0, 1]; // Jacobian of φ¯¹ at w = 0
// full Jtrans:
// +cos(w*α/β) sin(w*α/β) 0
// -sin(w*α/β) cos(w*α/β) 0
// -α*(u*sin(w*α/β) + v*cos(w*α/β))/β α*(u*cos(w*α/β) - v*sin(w*α/β))/β 1
Jtrans[] = Tensor[1, 0, 0, 0, 1, 0, -alpha/beta*Y[], alpha/beta*X[], 1]; // Transposed Jacobian of φ¯¹ at w = 0
// full Jinv:
// cos(w*α/β) sin(w*α/β) v*α/β
// -sin(w*α/β) cos(w*α/β) -u*α/β
// 0 0 1
Jinv[] = Tensor[1, 0, alpha/beta*Y[], 0, 1, -alpha/beta*X[], 0, 0, 1]; // inverse Jacobian of φ¯¹ at w = 0
// full Jinvtrans:
// cos(w*α/β) -sin(w*α/β) 0
// sin(w*α/β) cos(w*α/β) 0
// v*α/β -u*α/β 1
Jinvtrans[] = Tensor[1, 0, 0, 0, 1, 0, alpha/beta*Y[], -alpha/beta*X[], 1]; // Transposed inverse Jacobian of φ¯¹ at w = 0
// full T: T = Jtrans * J
// 1 0 -v*α/β
// 0 1 u*α/β
// -v*α/β u*α/β (u^2*α^2 + v^2*α^2 + β^2)/β^2
T[] = Tensor[1, 0, -alpha/beta*Y[], 0, 1, alpha/beta*X[], -alpha/beta*Y[], alpha/beta*X[], (X[]^2*alpha^2 + Y[]^2*alpha^2 + beta^2)/beta^2];
// Tinv at w = 0: (Jinv * Jinvtrans)(w=0)
// v^2*α^2/β^2 + 1 -u*v*α^2/β^2 v*α/β
// -u*v*α^2/β^2 u^2*α^2/β^2 + 1 -u*α/β
// v*α/β -u*α/β 1
// Tinv[] = Tensor[1, 0, -alpha/beta*Y[], 0, 1, alpha/beta*X[], -alpha/beta*Y[], alpha/beta*X[], (X[]^2*alpha^2 + Y[]^2*alpha^2 + beta^2)/beta^2];
Tinv[] = Tensor[(alpha/beta)^2 * Y[]^2 + 1, -(alpha/beta)^2 * X[] * Y[], alpha/beta*Y[], -(alpha/beta)^2 * X[] * Y[], (alpha/beta)^2 * X[]^2 + 1, -alpha/beta*X[], alpha/beta*Y[], -alpha/beta*X[], 1];
// material transformation:
mu_uvw[] = mu[] * Tinv[]; // transformation μ_xyz → μ_uvw, see Section III.A in the paper
rho_uvw[] = rho[] * T[]; // transformation ρ_xyz → ρ_uvw, see Section III.A in the paper
} // end of Function block
///////////////////////////////////////////////
//////////// step 03. Jacobian ////////////////
///////////////////////////////////////////////
Jacobian {
{ Name Vol ;
Case { { Region All ; Jacobian Vol ; } }
}
{ Name Sur ;
Case { { Region All ; Jacobian Sur ; } }
}
{ Name Lin ;
Case { { Region All ; Jacobian Lin ; } }
}
} // end of Jacobian block
///////////////////////////////////////////////
//////////// step 04. Integration /////////////
///////////////////////////////////////////////
Integration {
{ Name Int ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 15 ; } // valid choices are 1, 4, 5, 15, 16, ...
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 6 ; }
} // end of inner Case
} // end of "Type Gauss"
} // end of outer Case
} // end of "Int"
} // end of Integration block
///////////////////////////////////////////////
//////////// step 05. Constraint //////////////
///////////////////////////////////////////////
Constraint {
// constraint for magnetic scalar potential:
{ Name phi_constraint ;
// not needed here
} // end Constraint phi_constraint
// constraint for w-component of magnetic field:
{ Name H_w_constraint ;
Case {
{Region Omega_i ; Value Sqrt[2] * 1.0/(0.03183098861837907*2*Pi) -2.335961833704558 ;}
{Region BdOmega_c ; Value Sqrt[2] * 1.0/(0.03183098861837907*2*Pi) -2.335961833704558 ;}
} // end of Case
} // end Constraint H_w_constraint
{ Name VoltageTO ;
Case {
}
} // end of Constraint VoltageTO
// constraint for cohomomology basis functions, which are part of the function space for h_uv:
{ Name CurrentTO ;
Case {
{ Region Cut1TO; Value Itot[] ; }
{ Region Cut2TO; Value Itot[] ; }
{ Region Cut3TO; Value Itot[] ; }
{ Region Cut4TO; Value Itot[] ; }
{ Region Cut5TO; Value Itot[] ; }
{ Region Cut6TO; Value Itot[] ; }
{ Region Cut7TO; Value Itot[] ; }
{ Region Cut8TO; Value Itot[] ; }
{ Region Cut9TO; Value Itot[] ; }
{ Region Cut10TO; Value Itot[] ; }
{ Region Cut11TO; Value Itot[] ; }
{ Region Cut12TO; Value Itot[] ; }
{ Region Cut13TO; Value Itot[] ; }
} // end Case
} // end of Constraint CurrentTO
} // end of Constraint block
///////////////////////////////////////////////
////////// step 06. FunctionSpace /////////////
///////////////////////////////////////////////
FunctionSpace {
// function space for h_w:
{ Name V_w; Type Form1P;
BasisFunction {
{ Name sn; NameOfCoef hn; Function BF_PerpendicularEdge;
Support Omega_c; Entity NodesOf[All, Not BdOmega_c]; }
{ Name rn; NameOfCoef occi; Function BF_RegionZ;
Support Region[Omega_i]; Entity Region[Omega_i]; }
{ Name rn; NameOfCoef occi2; Function BF_GroupOfPerpendicularEdges;
Support Omega_c; Entity GroupsOfNodesOf[BdOmega_c]; }
}
Constraint {
{ NameOfCoef occi ; EntityType Region ; NameOfConstraint H_w_constraint ; }
{ NameOfCoef occi2 ; EntityType GroupsOfNodesOf ; NameOfConstraint H_w_constraint ; }
}
} // end of FunctionSpace V_w
// function space for h_uv:
{ Name V_uv; Type Form1;
BasisFunction {
{ Name sn; NameOfCoef phin; Function BF_GradNode;
Support Omega_i; Entity NodesOf[Omega_i]; }
{ Name sn; NameOfCoef phin2; Function BF_GroupOfEdges;
Support Omega_c; Entity GroupsOfEdgesOnNodesOf[BdOmega_c]; }
{ Name se; NameOfCoef t; Function BF_Edge;
Support Omega_c; Entity EdgesOf[All, Not BdOmega_c]; }
{ Name sc1; NameOfCoef I1; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut1TO]; }
{ Name sc2; NameOfCoef I2; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut2TO]; }
{ Name sc3; NameOfCoef I3; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut3TO]; }
{ Name sc4; NameOfCoef I4; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut4TO]; }
{ Name sc5; NameOfCoef I5; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut5TO]; }
{ Name sc6; NameOfCoef I6; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut6TO]; }
{ Name sc7; NameOfCoef I7; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut7TO]; }
{ Name sc8; NameOfCoef I8; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut8TO]; }
{ Name sc9; NameOfCoef I9; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut9TO]; }
{ Name sc10; NameOfCoef I10; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut10TO]; }
{ Name sc11; NameOfCoef I11; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut11TO]; }
{ Name sc12; NameOfCoef I12; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut12TO]; }
{ Name sc13; NameOfCoef I13; Function BF_GroupOfEdges;
Support Omega; Entity GroupsOfEdgesOf[Cut13TO]; }
} // end of BasisFunction
GlobalQuantity {
{ Name Current1 ; Type AliasOf ; NameOfCoef I1 ; }
{ Name Voltage1 ; Type AssociatedWith ; NameOfCoef I1 ; }
{ Name Current2 ; Type AliasOf ; NameOfCoef I2 ; }
{ Name Voltage2 ; Type AssociatedWith ; NameOfCoef I2 ; }
{ Name Current3 ; Type AliasOf ; NameOfCoef I3 ; }
{ Name Voltage3 ; Type AssociatedWith ; NameOfCoef I3 ; }
{ Name Current4 ; Type AliasOf ; NameOfCoef I4 ; }
{ Name Voltage4 ; Type AssociatedWith ; NameOfCoef I4 ; }
{ Name Current5 ; Type AliasOf ; NameOfCoef I5 ; }
{ Name Voltage5 ; Type AssociatedWith ; NameOfCoef I5 ; }
{ Name Current6 ; Type AliasOf ; NameOfCoef I6 ; }
{ Name Voltage6 ; Type AssociatedWith ; NameOfCoef I6 ; }
{ Name Current7 ; Type AliasOf ; NameOfCoef I7 ; }
{ Name Voltage7 ; Type AssociatedWith ; NameOfCoef I7 ; }
{ Name Current8 ; Type AliasOf ; NameOfCoef I8 ; }
{ Name Voltage8 ; Type AssociatedWith ; NameOfCoef I8 ; }
{ Name Current9 ; Type AliasOf ; NameOfCoef I9 ; }
{ Name Voltage9 ; Type AssociatedWith ; NameOfCoef I9 ; }
{ Name Current10 ; Type AliasOf ; NameOfCoef I10 ; }
{ Name Voltage10 ; Type AssociatedWith ; NameOfCoef I10 ; }
{ Name Current11 ; Type AliasOf ; NameOfCoef I11 ; }
{ Name Voltage11 ; Type AssociatedWith ; NameOfCoef I11 ; }
{ Name Current12 ; Type AliasOf ; NameOfCoef I12 ; }
{ Name Voltage12 ; Type AssociatedWith ; NameOfCoef I12 ; }
{ Name Current13 ; Type AliasOf ; NameOfCoef I13 ; }
{ Name Voltage13 ; Type AssociatedWith ; NameOfCoef I13 ; }
} // end of GlobalQuantity
Constraint {
{ NameOfCoef Current1 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage1 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current2 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage2 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current3 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage3 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current4 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage4 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current5 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage5 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current6 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage6 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current7 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage7 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current8 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage8 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current9 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage9 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current10 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage10 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current11 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage11 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current12 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage12 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef Current13 ; EntityType GroupsOfEdgesOf ; NameOfConstraint CurrentTO ; }
{ NameOfCoef Voltage13 ; EntityType GroupsOfEdgesOf ; NameOfConstraint VoltageTO ; }
{ NameOfCoef phin ; EntityType NodesOf; NameOfConstraint phi_constraint; }
{ NameOfCoef phin2; EntityType NodesOf; NameOfConstraint phi_constraint; }
} // end of Constraint block (in FunctionSpace block)
} // end of FunctionSpace V_uv
} // end of FunctionSpace
///////////////////////////////////////////////
////////// step 07. Formulation ///////////////
///////////////////////////////////////////////
Formulation {
{ Name MagDyn_htot_full; Type FemEquation;
Quantity {
{ Name h_uv; Type Local; NameOfSpace V_uv; }
{ Name h_w; Type Local; NameOfSpace V_w; }
{ Name I1; Type Global; NameOfSpace V_uv[Current1]; }
{ Name V1; Type Global; NameOfSpace V_uv[Voltage1]; }
{ Name I2; Type Global; NameOfSpace V_uv[Current2]; }
{ Name V2; Type Global; NameOfSpace V_uv[Voltage2]; }
{ Name I3; Type Global; NameOfSpace V_uv[Current3]; }
{ Name V3; Type Global; NameOfSpace V_uv[Voltage3]; }
{ Name I4; Type Global; NameOfSpace V_uv[Current4]; }
{ Name V4; Type Global; NameOfSpace V_uv[Voltage4]; }
{ Name I5; Type Global; NameOfSpace V_uv[Current5]; }
{ Name V5; Type Global; NameOfSpace V_uv[Voltage5]; }
{ Name I6; Type Global; NameOfSpace V_uv[Current6]; }
{ Name V6; Type Global; NameOfSpace V_uv[Voltage6]; }
{ Name I7; Type Global; NameOfSpace V_uv[Current7]; }
{ Name V7; Type Global; NameOfSpace V_uv[Voltage7]; }
{ Name I8; Type Global; NameOfSpace V_uv[Current8]; }
{ Name V8; Type Global; NameOfSpace V_uv[Voltage8]; }
{ Name I9; Type Global; NameOfSpace V_uv[Current9]; }
{ Name V9; Type Global; NameOfSpace V_uv[Voltage9]; }
{ Name I10; Type Global; NameOfSpace V_uv[Current10]; }
{ Name V10; Type Global; NameOfSpace V_uv[Voltage10]; }
{ Name I11; Type Global; NameOfSpace V_uv[Current11]; }
{ Name V11; Type Global; NameOfSpace V_uv[Voltage11]; }
{ Name I12; Type Global; NameOfSpace V_uv[Current12]; }
{ Name V12; Type Global; NameOfSpace V_uv[Voltage12]; }
{ Name I13; Type Global; NameOfSpace V_uv[Current13]; }
{ Name V13; Type Global; NameOfSpace V_uv[Voltage13]; }
} // end of Quantity
Equation {
// Weak formumation reads:
// Seek H_uvw ∈ V_uvw(Ω), s.t. ∀ H_uvw' ∈ V₀_uvw(Ω):
// ∫_Ω_uvw ∂ₜ μ_uvw H_uvw ⋅ H_uvw' dV + ∫_Ω_c_uvw ρ_uvw curl_xyz H_uvw ⋅ curl_xyz H_uvw' dV = 0
// ∫_Ω_uvw ∂ₜ μ_uvw H_uvw ⋅ H_uvw' dV:
Galerkin { DtDof [ mu_uvw[] * Dof{h_uv} , {h_uv} ];
In Omega; Integration Int; Jacobian Vol; } // combination 1: uv-components tested by uv-components of test functions
Galerkin { DtDof [ mu_uvw[] * Dof{h_w} , {h_w} ] ;
In Omega; Integration Int; Jacobian Vol; } // combination 2: w-component tested by w-component of test functions
Galerkin { DtDof [ mu_uvw[] * Dof{h_w} , {h_uv} ] ;
In Omega; Integration Int; Jacobian Vol; } // combination 3: w-components tested by uv-components of test functions
Galerkin { DtDof [ mu_uvw[] * Dof{h_uv} , {h_w} ] ;
In Omega; Integration Int; Jacobian Vol; } // combination 4: uv-components tested by w-component of test functions
// ∫_Ω_c_uvw ρ_uvw curl_xyz H_uvw ⋅ curl_xyz H_uvw' dV:
Galerkin { [ rho_uvw[] * Dof{d h_uv} , {d h_uv} ];
In Omega_c; Integration Int; Jacobian Vol; } // combination 1: uv-components tested by uv-components of test functions
Galerkin { [ rho_uvw[] * Dof{d h_w} , {d h_w} ];
In Omega_c; Integration Int; Jacobian Vol; } // combination 2: w-component tested by w-component of test functions
Galerkin { [ rho_uvw[] * Dof{d h_w} , {d h_uv} ];
In Omega_c; Integration Int; Jacobian Vol; } // combination 3: w-components tested by uv-components of test functions
Galerkin { [ rho_uvw[] * Dof{d h_uv} , {d h_w} ];
In Omega_c; Integration Int; Jacobian Vol; } // combination 4: uv-components tested by w-component of test functions
GlobalTerm { [ Dof{V1} , {I1} ] ; In Cut1TO ; }
GlobalTerm { [ Dof{V2} , {I2} ] ; In Cut2TO ; }
GlobalTerm { [ Dof{V3} , {I3} ] ; In Cut3TO ; }
GlobalTerm { [ Dof{V4} , {I4} ] ; In Cut4TO ; }
GlobalTerm { [ Dof{V5} , {I5} ] ; In Cut5TO ; }
GlobalTerm { [ Dof{V6} , {I6} ] ; In Cut6TO ; }
GlobalTerm { [ Dof{V7} , {I7} ] ; In Cut7TO ; }
GlobalTerm { [ Dof{V8} , {I8} ] ; In Cut8TO ; }
GlobalTerm { [ Dof{V9} , {I9} ] ; In Cut9TO ; }
GlobalTerm { [ Dof{V10} , {I10} ] ; In Cut10TO ; }
GlobalTerm { [ Dof{V11} , {I11} ] ; In Cut11TO ; }
GlobalTerm { [ Dof{V12} , {I12} ] ; In Cut12TO ; }
GlobalTerm { [ Dof{V13} , {I13} ] ; In Cut13TO ; }
} // end of Equation
} // end of MagDyn_htot_full
} // end of Formulation
///////////////////////////////////////////////
////////// step 08. Resolution ////////////////
///////////////////////////////////////////////
Resolution {
{ Name MagDynTOComplex;
System {
{ Name A; NameOfFormulation MagDyn_htot_full;
Type ComplexValue; Frequency Freq;} // "∂ₜ → jω"
}
Operation {
Generate[A]; Solve[A]; SaveSolution[A];
}
} // end of MagDynTOComplex block
} // end of Resolution block
///////////////////////////////////////////////
////////// step 09. PostProcessing ////////////
///////////////////////////////////////////////
PostProcessing {
// h-formulation, total field
{ Name MagDyn_htot_full; NameOfFormulation MagDyn_htot_full;
Quantity {
// magnetic scalar potential:
{ Name phi; Value{ Local{ [ {dInv h_uv} ] ; In Omega_i; Jacobian Vol; } } }
// magnetic field in helicoidal and cartesian coords:
{ Name h_uvw; Value{ Local{ [ ({h_uv} + {h_w}) ] ; In Omega; Jacobian Vol; } } }
{ Name h_xyz; Value{ Local{ [ Jinvtrans[]*({h_uv} + {h_w}) ] ; In Omega; Jacobian Vol; } } }
{ Name h_xyz_abs; Value{ Local{ [Norm[ Jinvtrans[]*({h_uv} + {h_w}) ]] ; In Omega; Jacobian Vol; } } }
// current density in helicoidal and cartesian coords:
{ Name j_uvw; Value{ Local{ [ ({d h_uv} + {d h_w}) ] ; In Omega_c; Jacobian Vol; } } }
{ Name j_xyz; Value{ Local{ [J[] * ({d h_uv} + {d h_w}) ] ; In Omega_c; Jacobian Vol; } } }
{ Name j_xyz_abs; Value{ Local{ [Norm[ J[] * ({d h_uv} + {d h_w}) ]] ; In Omega_c; Jacobian Vol; } } }
// loss density:
{ Name power_losses_xyz; Value{ Local{ [ rho[] * SquNorm [ J[] * ({d h_uv} + {d h_w}) ] ]; In Omega_c; Jacobian Vol; } } }
// global quantities:
{ Name dissPower;
Value{
Integral{
[rho[]/2 * SquNorm[ J[] * ({d h_uv} + {d h_w}) ]]; In Omega_c; Integration Int; Jacobian Vol;
} // end of Integral
} // end of Value
}
{ Name I1 ; Value { Term { [ {I1} ] ; In Cut1TO ; } } }
{ Name V1 ; Value { Term { [ {V1} ] ; In Cut1TO ; } } }
{ Name I2 ; Value { Term { [ {I2} ] ; In Cut2TO ; } } }
{ Name V2 ; Value { Term { [ {V2} ] ; In Cut2TO ; } } }
{ Name I3 ; Value { Term { [ {I3} ] ; In Cut3TO ; } } }
{ Name V3 ; Value { Term { [ {V3} ] ; In Cut3TO ; } } }
{ Name I4 ; Value { Term { [ {I4} ] ; In Cut4TO ; } } }
{ Name V4 ; Value { Term { [ {V4} ] ; In Cut4TO ; } } }
{ Name I5 ; Value { Term { [ {I5} ] ; In Cut5TO ; } } }
{ Name V5 ; Value { Term { [ {V5} ] ; In Cut5TO ; } } }
{ Name I6 ; Value { Term { [ {I6} ] ; In Cut6TO ; } } }
{ Name V6 ; Value { Term { [ {V6} ] ; In Cut6TO ; } } }
{ Name I7 ; Value { Term { [ {I7} ] ; In Cut7TO ; } } }
{ Name V7 ; Value { Term { [ {V7} ] ; In Cut7TO ; } } }
{ Name I8 ; Value { Term { [ {I8} ] ; In Cut8TO ; } } }
{ Name V8 ; Value { Term { [ {V8} ] ; In Cut8TO ; } } }
{ Name I9 ; Value { Term { [ {I9} ] ; In Cut9TO ; } } }
{ Name V9 ; Value { Term { [ {V9} ] ; In Cut9TO ; } } }
{ Name I10 ; Value { Term { [ {I10} ] ; In Cut10TO ; } } }
{ Name V10 ; Value { Term { [ {V10} ] ; In Cut10TO ; } } }
{ Name I11 ; Value { Term { [ {I11} ] ; In Cut11TO ; } } }
{ Name V11 ; Value { Term { [ {V11} ] ; In Cut11TO ; } } }
{ Name I12 ; Value { Term { [ {I12} ] ; In Cut12TO ; } } }
{ Name V12 ; Value { Term { [ {V12} ] ; In Cut12TO ; } } }
{ Name I13 ; Value { Term { [ {I13} ] ; In Cut13TO ; } } }
{ Name V13 ; Value { Term { [ {V13} ] ; In Cut13TO ; } } }
} // end of Quantity
} // end of MagDyn_htot_full PostProcessing
} // end of PostProcessing
///////////////////////////////////////////////
////////// step 10. PostOperation /////////////
///////////////////////////////////////////////
PostOperation {
{ Name MagDyn; NameOfPostProcessing MagDyn_htot_full ;
Operation{
Print[ phi, OnElementsOf Omega , File "phi.pos", Name "φ_magnetic [A]" ];
Print[ h_xyz, OnElementsOf Omega , File "h_xyz.pos", Name "H(xyz) [A/m]" ];
Print[ h_uvw, OnElementsOf Omega , File "h_uvw.pos", Name "H(uvw) [A/m]" ];
Print[ h_xyz_abs, OnElementsOf Omega, File "h_xyz_abs.pos", Name "|H|(xyz) [A/m]" ];
Print[ j_xyz, OnElementsOf Omega_c , File "j_xyz.pos", Name "J(xyz) [A/m²]" ];
Print[ j_uvw, OnElementsOf Omega_c , File "j_uvw.pos", Name "J(uvw) [A/m²]" ];
Print[ j_xyz_abs, OnElementsOf Omega_c , File "j_xyz_abs.pos", Name "|J|(xyz) [A/m²]" ];
Print[ power_losses_xyz, OnElementsOf Omega_c , File "p_loss_xyz.pos", Name "p_loss(xyz) [W/m]" ];
// Print[ phi, OnLine {{-0.045, 0, 0}{0.045, 0, 0}} {100}, File "φ_table.txt"];
// Print[ h_xyz, OnLine {{-0.045, 0, 0}{0.045, 0, 0}} {1100}, File "h_xyz_table.txt"];
// Print[ h_uvw, OnLine {{-0.045, 0, 0}{0.045, 0, 0}} {1100}, File "h_uvw_table.txt"];
// Print[ j_xyz, OnLine {{-0.045, 0, 0}{0.045, 0, 0}} {1100}, File "j_xyz_table.txt"];
// Print[ j_uvw, OnLine {{-0.045, 0, 0}{0.045, 0, 0}} {1100}, File "j_uvw_table.txt"];
// print "global" results:
Print[ dissPower[Omega_c], OnGlobal, Format Table, File "jouleLosses.txt" ];
Print[I1, OnRegion Cut1TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V1, OnRegion Cut1TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I2, OnRegion Cut2TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V2, OnRegion Cut2TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I3, OnRegion Cut3TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V3, OnRegion Cut3TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I4, OnRegion Cut4TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V4, OnRegion Cut4TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I5, OnRegion Cut5TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V5, OnRegion Cut5TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I6, OnRegion Cut6TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V6, OnRegion Cut6TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I7, OnRegion Cut7TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V7, OnRegion Cut7TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I8, OnRegion Cut8TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V8, OnRegion Cut8TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I9, OnRegion Cut9TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V9, OnRegion Cut9TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I10, OnRegion Cut10TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V10, OnRegion Cut10TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I11, OnRegion Cut11TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V11, OnRegion Cut11TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I12, OnRegion Cut12TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V12, OnRegion Cut12TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
Print[I13, OnRegion Cut13TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each current
Print[V13, OnRegion Cut13TO, Format FrequencyTable, File >> "currentsAndVoltages.txt"]; // real and imaginary part of each voltage (drop)
} // end of Operation
} // end of PostProcessing "MagDyn"
} // end of PostOperation
// some GetDP constants:
DefineConstant[
R_ = {"MagDynTOComplex" , Name "GetDP/1ResolutionChoices" , Visible 1},
C_ = {"-solve -v2 -pos" , Name "GetDP/9ComputeCommand" , Visible 1},
P_ = {"MagDyn" , Name "GetDP/2PostOperationChoices", Visible 1}
];
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment