Skip to content
Snippets Groups Projects
Commit bb10fad3 authored by François Henrotte's avatar François Henrotte
Browse files

new models

parent 726809aa
No related branches found
No related tags found
No related merge requests found
Pipeline #5681 canceled
Stage: test
Hide whitespace changes
Inline Side-by-side
Conductors3D/w2d.geo 0 → 100644
+ 125
− 0
View file @ bb10fad3
Include "wires_common.pro";
// global Gmsh options
Mesh.Optimize = 1; // optimize quality of tetrahedra
Mesh.VolumeEdges = 0; // hide volume edges
Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
lc1 = (Flag_Thin==0) ? MeshSizeWire*mm : rs ;
lc2 = box/10; // mesh size at outer surface
llWires[] = {};
centerWires[] = {};
surfWires[] = {};
If( Flag_Thin != 1)
// then circles are in the geometry
// their radius is R= WX~{i} in the real case (Flag_Thin == 0)
// or R=rs in the "regular sleeve case" (Flag_Thin == 2)
For i In {1:NumWires}
R = (Flag_Thin == 0 ) ? rw : rs;
pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
p1 = newp; Point(p1) = {WX~{i}+R, WY~{i} , 0, lc1};
p2 = newp; Point(p2) = {WX~{i} , WY~{i}+R, 0, lc1};
p3 = newp; Point(p3) = {WX~{i}-R, WY~{i} , 0, lc1};
p4 = newp; Point(p4) = {WX~{i} , WY~{i}-R, 0, lc1};
c1 = newl; Circle(c1) = {p1,pC,p2};
c2 = newl; Circle(c2) = {p2,pC,p3};
c3 = newl; Circle(c3) = {p3,pC,p4};
c4 = newl; Circle(c4) = {p4,pC,p1};
l1 = newl; Line(l1) = {pC,p1};
l2 = newl; Line(l2) = {pC,p2};
l3 = newl; Line(l3) = {pC,p3};
l4 = newl; Line(l4) = {pC,p4};
ll1 = newll; Line Loop(ll1) = {l1,c1,-l2};
s1 = news; Plane Surface(s1) = {ll1};
ll2 = newll; Line Loop(ll2) = {l2,c2,-l3};
s2 = news; Plane Surface(s2) = {ll2};
ll3 = newll; Line Loop(ll3) = {l3,c3,-l4};
s3 = news; Plane Surface(s3) = {ll3};
ll4 = newll; Line Loop(ll4) = {l4,c4,-l1};
s4 = news; Plane Surface(s4) = {ll4};
llWires[] += {ll1,ll2,ll3,ll4};
surfWires~{i-1}[] = {s1,s2,s3,s4};
surfWires[] += surfWires~{i-1}[];
centerWires[] += {pC};
EndFor
Else
For i In {1:NumWires}
pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
centerWires[] += {pC};
EndFor
EndIf
// create air box around wires
BoundingBox; // recompute model bounding box
cx = (General.MinX + General.MaxX) / 2;
cy = (General.MinY + General.MaxY) / 2;
cz = (General.MinZ + General.MaxZ) / 2;
/* lx = 2*inf + General.MaxX - General.MinX; */
/* ly = 2*inf + General.MaxY - General.MinZ; */
lx = 2*box;
ly = 2*box;
p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, 0, lc2};
p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, 0, lc2};
p3 = newp; Point (p3) = {cx+lx/2, cy+ly/2, 0, lc2};
p4 = newp; Point (p4) = {cx-lx/2, cy+ly/2, 0, lc2};
l1 = newl; Line(l1) = {p1, p2};
l2 = newl; Line(l2) = {p2, p3};
l3 = newl; Line(l3) = {p3, p4};
l4 = newl; Line(l4) = {p4, p1};
ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
Physical Line("INF", 2) = { l1, l2, l3, l4 };
If( !Flag_3D )
If( Flag_Thin == 0 ) // Real
s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
Physical Surface("AIR", 1) = { s1 };
For i In {1:NumWires}
Physical Surface(Sprintf("VWIRE_%g",i), 10+i) = { surfWires~{i-1}[] };
Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
EndFor
ElseIf( Flag_Thin == 1 ) // 1D raw sleeve
s1 = news; Plane Surface(s1) = {ll1};
Physical Surface("AIR", 1) = { s1 };
For i In {1:NumWires}
Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
Point { centerWires[i-1] } In Surface { s1 };
EndFor
Else // 1D regular sleeve
s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
Physical Surface("AIR", 1) = { s1, surfWires[] };
For i In {1:NumWires}
Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
Point { centerWires[i-1] } In Surface { s1 };
EndFor
EndIf
Else
For i In {1:NumWires}
e[] = Extrude {0,0,Lz} { Surface{ surfWires[i-1] } ; Layers{1}; Recombine; } ;
Printf("%g", #e[]);
Printf("%g %g %g %g %g %g", e[0], e[1], e[2], e[3], e[4], e[5]);
Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[1] };
Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires[i-1] };
Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[0] };
e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; Layers{1}; Recombine; } ;
Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
EndFor
e[] = Extrude {0,0,Lz} { Surface{ s1 } ; Layers{1}; Recombine; } ;
Physical Volume ("AIR", 1) = { e[1] };
Physical Surface("INF", 2) = { s1, e[0], e[2], e[3], e[4], e[5] };
//Physical Surface("INF", 1) = { CombinedBoundary{ Volume{ e[1]}; } };
EndIf
Conductors3D/w2d.pro 0 → 100644
+ 707
− 0
View file @ bb10fad3
Include "wires_common.pro";
Struct S::SOLVER [ Enum, mumps, gmres, fgmres,bcgs ];
DefineConstant[
FileId = ""
DataDir = StrCat["data/" , FileId]
ElekinDir = StrCat["Elekin_Line_Const/" , FileId]
MagDynDir = StrCat["MagDyn_Line_Const/" , FileId]
FullWaveDir = StrCat["FullWave_Line_Const/", FileId]
type_solver = {
S::SOLVER.mumps,
Choices {
S::SOLVER.mumps ="mumps",
S::SOLVER.gmres="gmres",
S::SOLVER.fgmres= "fgmres",
S::SOLVER.bcgs="bcgs"
},
Name "Solver/type_solver"
}
];
DefineConstant[
R_ = {"Dynamic", Name "GetDP/1ResolutionChoices", Visible 0},
C_ = {"-sol -pos", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
ResDir = "Res/"
];
Group{
// Geometrical regions
AIR = Region[ 1 ];
INF = Region[ 2 ];
VWIRES = Region[ {} ];
LWIRES = Region[ {} ];
For i In {1:NumWires}
LWIRE~{i} = Region[ (50+i) ];
LWIRES += Region[ LWIRE~{i} ];
If( Flag_Thin == 0 )
VWIRE~{i} = Region[ (10+i) ];
VWIRES += Region[ VWIRE~{i} ];
SWIRE~{i} = Region[ {} ];
Else
VWIRE~{i} = Region[ {} ];
SWIRE~{i} = ElementsOf[ AIR, OnOneSideOf LWIRE~{i} ];
EndIf
EndFor
// Abstract regions
// Cette formulation peut contenir des conducteurs volumiques Vol_C
// et liniques lVol_C.
If( Flag_Thin == 0 )
Vol_C = Region[ VWIRES ];
Vol_CC = Region[ AIR ];
lVol_C = Region[ {} ];
Else
Vol_C = Region[ {} ];
Vol_CC = Region[ {VWIRES, AIR} ];
lVol_C = Region[ LWIRES ];
EndIf
Vol_nu = Region[ {Vol_C, Vol_CC, lVol_C} ];
Sur_Dirichlet_a = Region[ INF ];
// 3D
If( Flag_3D )
Sur_Dirichlet_a += Region[ { BOT, TOP} ];
EndIf
Sur_Complete_Tree = Region[ { lVol_C, Sur_Dirichlet_a } ];
Electrodes = Region[ {} ];
Dom_Hcurl_a = Region[ Vol_nu ];
Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf lVol_C ]; // Sleeve
Dom_Hgrad_u = Region[ Vol_C ];
Dom_Hregion_i = Region[ lVol_C ];
// For integration on the volume elements (excluding line elements) of the sleeve
Vol_Sleeve = ElementsOf[ {Vol_CC}, OnOneSideOf lVol_C ];
}
Function{
omega = 2*Pi*Freq;
mu0 = Pi*4e-7;
mur = 1.;
sigma = 5.96e7;
i[] = Complex[0,1];
mu[] = mu0;
nu[] = 1/mu[];
sigma[] = sigma;
skin_depth = Sqrt[ 2/(omega*sigma*mu0*mur) ];
tau[] = Complex[-1,1]/skin_depth*rw;
J0[] = JnComplex[0,$1];
J1[] = JnComplex[1,$1];
dJ1[] = JnComplex[0,$1]-JnComplex[1,$1]/$1;
For i In {1:NumWires}
// Current density
J[ Region[ {VWIRE~{i}, LWIRE~{i}} ] ] =
Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
// local radial coordinate with origin on wire i
R~{i}[] = Hypot[ X[]-WX~{i} , Y[]-WY~{i} ];
EndFor
// shape function of the current, wI[r]
wI[] = tau[] * J0[tau[]*$1/rw] / J1[tau[]] /(2*A_c);
//radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
Analytic_A[] = // per Amp
(($1>rw) ? mu0/(2*Pi) * Log[rs/$1] :
mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) );
AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
(($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
// Impedance of thin wire p.u. length
R_DC = 1./(sigma*A_c);
Analytic_R[] = Re[ wI[rw]/sigma ];
Analytic_L[] = mu0/(2*Pi) * Log[rs/rw] + Re[ wI[rw]/(i[]*omega*sigma) ];
If( NumWires == 1 )
Exact_A[] = WI_1*(Analytic_A[R_1[]] + mu0/(2*Pi) * Log[rout/rs]) ;
Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] : 0 ) ;
EndIf
If ( NumWires == 2 )
Exact_A[] = WI_1*Analytic_A[R_1[]] + WI_2*Analytic_A[R_2[]]
+ (WI_1+WI_2)*mu0/(2*Pi) * NumWires * Log[rout/rs] ;
Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] :
((R_2[]<rs) ? WI_2*Analytic_A[R_2[]] : 0 ));
EndIf
If ( NumWires == 3 )
Exact_a[] = WI_1*Analytic_A[R_1[]] + WI_2*Analytic_A[R_2[]]
+ WI_3*Analytic_A[R_3[]] + (WI1_+WI_2+WI_3)*mu0/(2*Pi) * NumWires * Log[rout/rs] ;
Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] :
((R_2[]<rs) ? WI_2*Analytic_A[R_2[]] :
((R_3[]<rs) ? WI_3*Analytic_A[R_3[]] : 0 )));
EndIf
}
Jacobian {
{ Name Vol ;
Case {
{ Region lVol_C ; Jacobian Lin ; Coefficient A_c; }
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 3 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Prism ; NumberOfPoints 21 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Pyramid ; NumberOfPoints 8 ; }
}
}
}
}
}
Constraint{
// 2D
{ Name Hcurl_a_2D_ac; Type Assign;
Case {
{ Region Region[ Sur_Dirichlet_a ]; Value 0; }
}
}
{ Name Impose_I ;
Case {
For i In {1:NumWires}
{ Region VWIRE~{i} ; Value WI~{i} ; }
{ Region LWIRE~{i} ; Value WI~{i} ; }
EndFor
}
}
// 3D
{ Name Hcurl_a_3D_ac; Type Assign;
Case {
{ Region Region[ Sur_Dirichlet_a ]; Value 0.; }
}
}
{ Name GaugeCondition_a; Type Assign;
Case {
{ Region Vol_CC; SubRegion Sur_Complete_Tree; Value 0.; }
}
}
}
FunctionSpace {
// 2D
{ Name Hcurl_a_2D; Type Form1P;
BasisFunction {
{ Name sc; NameOfCoef ac; Function BF_PerpendicularEdge;
Support Dom_Hcurl_a; Entity NodesOf[All]; }
{ Name sw; NameOfCoef aw; Function BF_PerpendicularEdge;
Support Dom_Hthin_a; Entity NodesOf[lVol_C]; }
}
SubSpace {
{ Name ac ; NameOfBasisFunction { sc } ; }
{ Name aw ; NameOfBasisFunction { sw } ; }
}
Constraint {
{ NameOfCoef ac;
EntityType Auto; NameOfConstraint Hcurl_a_2D_ac; }
}
}
{ Name Hgrad_u_2D ; Type Form1P ;
BasisFunction {
{ Name sr ; NameOfCoef ur ; Function BF_RegionZ ;
Support Dom_Hgrad_u ; Entity Vol_C ; }
}
GlobalQuantity {
{ Name U ; Type AliasOf ; NameOfCoef ur ; }
{ Name I ; Type AssociatedWith ; NameOfCoef ur ; }
}
Constraint {
{ NameOfCoef I ; EntityType GroupsOfNodesOf ; NameOfConstraint Impose_I ; }
}
}
{ Name Hregion_i_2D; Type Vector;
BasisFunction {
{ Name sr; NameOfCoef ir; Function BF_RegionZ;
Support Dom_Hregion_i; Entity lVol_C; }
}
GlobalQuantity {
{ Name Is; Type AliasOf; NameOfCoef ir; }
{ Name Us; Type AssociatedWith; NameOfCoef ir; }
}
Constraint {
{ NameOfCoef Is;
EntityType Region; NameOfConstraint Impose_I; }
}
}
// 3D
{ Name Hcurl_a_3D; Type Form1;
BasisFunction {
{ Name sc; NameOfCoef ac; Function BF_Edge;
Support Dom_Hcurl_a; Entity EdgesOf[All]; }
{ Name sw; NameOfCoef aw; Function BF_Edge;
Support Dom_Hthin_a; Entity EdgesOf[lVol_C]; }
}
SubSpace {
{ Name ac ; NameOfBasisFunction { sc } ; }
{ Name aw ; NameOfBasisFunction { sw } ; }
}
Constraint {
{ NameOfCoef ac;
EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
}
}
{ Name Hgrad_u_3D; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Node;
Support Dom_Hgrad_u; Entity NodesOf[ Vol_C, Not Electrodes]; }
{ Name sn2; NameOfCoef un2; Function BF_GroupOfNodes;
Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }
}
GlobalQuantity {
{ Name U; Type AliasOf ; NameOfCoef un2; }
{ Name I; Type AssociatedWith; NameOfCoef un2; }
}
Constraint {
//{ NameOfCoef U; EntityType Auto; NameOfConstraint Impose_U; }
{ NameOfCoef I; EntityType Auto; NameOfConstraint Impose_I; }
}
}
}
Formulation {
{ Name MagnetoDynamics; Type FemEquation;
If( !Flag_3D )
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_2D; }
{ Name ac; Type Local; NameOfSpace Hcurl_a_2D[ac]; }
{ Name aw; Type Local; NameOfSpace Hcurl_a_2D[aw]; }
{ Name dv; Type Local; NameOfSpace Hgrad_u_2D; }
{ Name U; Type Global; NameOfSpace Hgrad_u_2D [U]; }
{ Name I; Type Global; NameOfSpace Hgrad_u_2D [I]; }
{ Name i; Type Local; NameOfSpace Hregion_i_2D; }
{ Name Is; Type Global; NameOfSpace Hregion_i_2D [Is]; }
{ Name Us; Type Global; NameOfSpace Hregion_i_2D [Us]; }
}
Equation {
Integral { [ nu[] * Dof{d ac} , {d ac} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {ac} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ nu[] * Dof{d aw} , {d aw} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [-Dof{i}/A_c , {aw} ];
In lVol_C; Jacobian Vol; Integration I1; }
/*
Integral { [ -J[] , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
*/
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{dv} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {dv} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{dv} , {dv} ];
In Vol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [ Dof{I} , {U} ]; In Vol_C; }
GlobalTerm { [ Dof{Us} , {Is} ];
In lVol_C; }
GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
In lVol_C;}
Integral { DtDof [ Dof{ac}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1; }
If( Flag_Thin != 3)
Integral { DtDof [ -Dof{aw}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1;}
GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
In lVol_C;}
EndIf
}
Else // ############### 3D
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_3D; }
{ Name ac; Type Local; NameOfSpace Hcurl_a_3D[ac]; }
{ Name aw; Type Local; NameOfSpace Hcurl_a_3D[aw]; }
// { Name v; Type Local; NameOfSpace Hgrad_u_3D; }
// { Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
// { Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
}
Equation {
/*
Integral { [ nu[] * Dof{d ac} , {d ac} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ nu[] * Dof{d aw} , {d aw} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Vector[ 0, 0, I0/A_c] , {ac} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -Vector[ 0, 0, I0/A_c] , {aw} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -Vector[ 0, 0, I0/A_c] , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [Dof{U} , {I} ]; In Electrodes; }
*/
}
EndIf
}
}
Resolution {
{ Name Dynamic;
System {
{ Name S; NameOfFormulation MagnetoDynamics;
Type ComplexValue; Frequency Freq; }
}
Operation {
CreateDir[DataDir];
// SOLVER OPTION
If(type_solver == S::SOLVER.gmres)
SetGlobalSolverOptions[
"-pc_type ilu -ksp_type gmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
ElseIf(type_solver == S::SOLVER.fgmres)
SetGlobalSolverOptions[
"-pc_type ilu -ksp_type fgmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
ElseIf(type_solver == S::SOLVER.bcgs)
SetGlobalSolverOptions[
"-pc_type ilu -ksp_type bcgs -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
EndIf
SetGlobalSolverOptions["-petsc_prealloc 900"];
Generate[S]; Solve[S]; SaveSolution[S];
PostOperation[integaz];
If( Flag_Maps == 1 )
PostOperation[mapaz];
PostOperation[cutaz];
DeleteFile["U_full.dat"];
DeleteFile["U_raw.dat"];
DeleteFile["U_reg.dat"];
DeleteFile["U_naive.dat"];
DeleteFile["joule_full.dat"];
DeleteFile["joule_raw.dat"];
DeleteFile["joule_reg.dat"];
DeleteFile["joule_naive.dat"];
EndIf
If( Flag_Thin != 0 )
For i In {1:NumWires}
Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}},
Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
EndFor
EndIf
}
}
}
PostProcessing {
{ Name PostProcessings; NameOfFormulation MagnetoDynamics;
PostQuantity {
{ Name az;
Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name acz;
Value { Local { [ CompZ[ {ac}] ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name awz;
Value { Local { [ CompZ[ {aw}] ];
In Dom_Hthin_a; Jacobian Vol; }}}
{ Name acwz;
Value { Local { [ CompZ[ {ac} - {aw} ] ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name exact;
Value { Local { [ Exact_A[] ] ;
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name corr;
Value { Local { [ Correction_A[] ] ;
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name bcw;
Value { Local { [ {d ac} - {d aw} ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name j;
Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ];
In Vol_C; Jacobian Vol; }}}
{ Name un;
Value{ Local { [ 1 ] ;
In Vol_nu; Integration I1 ; Jacobian Vol; }}}
{ Name surf;
Value{ Integral{ [ 1 ] ;
In Vol_nu; Integration I1 ; Jacobian Vol; }}}
If( Flag_Thin == 0 )
{ Name flux;
Value { Integral { [ CompZ[ {a} ] / A_c ];
In Vol_C; Integration I1 ; Jacobian Vol; }}}
{ Name JouleLosses;
Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
{ Name U;
Value { Term { [ Cart2Pol [ {U} ] ]; In Vol_C; }}}
{ Name L;
Value { Term { [ Im [ {U}/{I}/omega ] ]; In Vol_C; }}}
Else
{ Name flux; // naive flux
Value { Integral { [ CompZ[ {ac} ] ];
In lVol_C; Integration I1 ; Jacobian Vol; }}}
{ Name JouleLosses;
Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
{ Name U;
Value { Term { [ Cart2Pol [ {Us} ] ] ; In lVol_C; }}}
{ Name L; // only valid with one current
Value { Term { [ Im [ {Us}/{Is}/omega ] ] ; In lVol_C; }}}
For i In {1:NumWires}
{ Name area~{i};
Value{ Integral{ [ 1 ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
{ Name intbx~{i}; // integral of b, 1st step for computing b average
Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
{ Name intby~{i}; // integral of b, 1st step for computing b average
Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
EndFor
EndIf
}
}
}
PostOperation mapaz UsingPost PostProcessings {
Print [az, OnElementsOf Vol_nu, File "az.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly];
If(Flag_Thin != 0)
Print [bcw, OnElementsOf Vol_Sleeve, File "b.pos"];
Print [acwz, OnElementsOf Vol_nu, File "acz.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly] ;
EndIf
}
NbPoints = 1000;
Xmax = 0.15*box;
PostOperation cutaz UsingPost PostProcessings {
// cut of az displayed in GUI
If( Flag_Thin == 0 )
Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut az full';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
Else
Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut az not corrected';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut az corrected';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
EndIf
// cuts in txt format for Gnuplot
Print [ exact, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_analytic.txt" ];
If( Flag_Thin == 0 )
Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_fullfem.txt" ];
Else
Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_thinfem.txt" ];
Print [ acz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_ac.txt" ];
Print [ awz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_aw.txt" ];
Print [ acwz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_acw.txt" ];
Print [ corr, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_correction.txt" ];
EndIf
}
PostOperation integaz UsingPost PostProcessings {
For i In {1:NumWires}
If( Flag_Thin == 0 )
Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
Print[ L, OnRegion Vol_C, File > "L_full.dat", Format Table,
SendToServer Sprintf["}Inductance/Wire %g/1Full",i]];
Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
Format Table, File > "joule_full.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
EndIf
If( Flag_Thin == 1 )
Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/2Raw",i]];
Print[ L, OnRegion lVol_C, File > "L_raw.dat", Format Table,
SendToServer Sprintf["}Inductance/Wire %g/2Raw",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_raw.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/2Raw",i]];
EndIf
If( Flag_Thin == 2 )
Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
Print[ L, OnRegion lVol_C, File > "L_reg.dat", Format Table,
SendToServer Sprintf["}Inductance/Wire %g/3Reg",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_reg.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
EndIf
If( Flag_Thin == 3 )
Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
Print[ L, OnRegion lVol_C, File > "L_naive.dat", Format Table,
SendToServer Sprintf["}Inductance/Wire %g/4naive",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_naive.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
EndIf
If( Flag_Thin != 0 )
Print[ area~{i}[SWIRE~{i}], OnGlobal,
Format Table, File >"zzz",
StoreInVariable $sarea~{i}];
Print[ intby~{i}[SWIRE~{i}], OnGlobal,
Format Table, File > "zzz",
StoreInVariable $intby~{i}];
EndIf
EndFor
}
/*{ Name flux_corr; // corrected flux
Value {
Integral { [ CompZ[ {ac} - {aw} ] ];
In lVol_C; Integration I1 ; Jacobian Vol; }
Integral { [ Analytic_L[]*{Is} ];
In lVol_C; Integration I1 ; Jacobian Vol; }
}
}*/
/* Print [awz, OnElementsOf Vol_nu, File "awz.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly] ;
Print [acz, OnElementsOf Vol_nu, File "acz.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly] ;
Print[ flux[ VWIRE~{i} ], OnGlobal,
Format Table, File > "flux_full.dat",
SendToServer Sprintf["}Flux/Wire %g/1Full",i]];
Print[ flux[ LWIRE~{i} ], OnGlobal,
Format Table, File > "flux_raw.dat",
SendToServer Sprintf["}Flux/Wire %g/2Raw",i]];
Print[ flux[ LWIRE~{i} ], OnGlobal,
Format Table, File > "flux_reg.dat",
SendToServer Sprintf["}Flux/Wire %g/4Reg",i]];
*/
/*
X = rw/skin_depth;
pI[] = (X/2)*Re[(1+i[])*(J_0[(1+i[])*X]/J_1[(1+i[])*X])];
Resist[] = Re[pI[]*R_dc];
// radial analytical solution with a zero flux condition imposed at r=rs
Analytic_a~{i}[] = mu0/(2*Pi) * WI~{i} *
( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] :
Log[rs/WR~{i}] + (J0[tau[]*R~{i}[]/rw] - J0[tau[]]) / J1[tau[]] / tau[] );
Analytic_L~{i}[] = mu0 / (2*Pi) *
( 2/tau[]/tau[]- J0[tau[]] / J1[tau[]] / tau[] - Log[rs/WR~{i}] ) ;
// Idem in static case (thus particular case of the above)
AnalyticStatic_a~{i}[] = mu0 * WI~{i} / (2*Pi) *
( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] : Log[rs/WR~{i}] + (1-(R~{i}[]/WR~{i})^2)/2. );
Analytic_e~{i}[] = mu0 * WI~{i} / (2*Pi) * i[] * omega
* J0[tau[]*R~{i}[]/WR~{i}] / J1[tau[]] / tau[] ;
Analytic_flux~{i}[] = mu0 * WI~{i} / (2*Pi) *
( i[] * (skin_depth/WR~{i})^2 + Log[rs/WR~{i}] - J0[tau[]] / J1[tau[]] / tau[] ) ;
*/
Conductors3D/w2s_common.pro 0 → 100644
+ 47
− 0
View file @ bb10fad3
mm = 1e-3;
deg = 180/Pi;
DefineConstant[
NumWires = {1, Name "Parameters/0Number of wires",
Choices {1="1", 2="2", 3="3"}, Visible 1}
Flag_Thin = {1, Name "Parameters/1Wire representation",
Choices {0="Real", 1="1D (raw sleeve)", 2="1D (regular sleeve)", 3="1D (naive)"}}
Flag_3D = {0, Name "Parameters/2Extruded model", Choices {0,1}, Visible 0}
Flag_Maps = {1, Name "Input/2Display maps", Choices {0,1}, Visible 1}
LogFreq = {4, Min 0, Max 8, Step 1, Name "Input/4Log of frequency"}
Freq = 10^LogFreq
WireRadius = {0.564189, Name "Parameters/5Wire radius [mm]"}
MeshSizeWire = {0.2, Name "Parameters/6Mesh size on wire [mm]"}
SleeveRadius = {2, Name "Parameters/7Sleeve radius [mm]"}
box = {100*mm, Name "Parameters/7box half-width [m]"}
rw = WireRadius*mm
rs = SleeveRadius*mm
A_c = Pi*rw^2 // wire cross section
rout = box
Lz = 1*mm
/* 'WireRadius' is the actual radius of the wire.
'MeshSizeWire' is the imposed mesh size at the nodes of the wire,
and thus also the practical approximative value of the radius of the sleeve.
The 'rs' and 'rw' parameters are assumed identical for all wires.
*/
];
Xlocs() = { 0, 5*mm, -5*mm};
Ylocs() = { 0, 0, 0};
For i In {1:NumWires}
DefineConstant[
WX~{i} = { Xlocs(i-1), // place wires symmetrically around x=0
Name Sprintf("Parameters/Wire %g/0X position [m]", i), Closed }
WY~{i} = { Ylocs(i-1), Name Sprintf("Parameters/Wire %g/0Y position [m]", i) }
WI~{i} = { 1, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
WP~{i} = { 0*NumWires*(i-1),
Name Sprintf("Parameters/Wire %g/3Phase [deg]", i) }
];
EndFor
Conductors3D/w3d.geo 0 → 100644
+ 149
− 0
View file @ bb10fad3
Include "w3d_common.pro";
// global Gmsh options
Mesh.Optimize = 1; // optimize quality of tetrahedra
Mesh.VolumeEdges = 0; // hide volume edges
Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
lc1 = (Flag_Thin) ? MeshSizeThinWire*mm : MeshSizeThinWire*mm ;
lc2 = box/10; // mesh size at outer surface
llWires[] = {};
centerWires[] = {};
surfWires[] = {};
If( !Flag_Thin || (Flag_Thin && Flag_Corr==2))
// then circles are in the geometry
// their radius is R=rw if Flag_Thin != 0
// or R=rs in the "regular sleeve case" (Flag_Corr == 2)
R = ( !Flag_Thin ) ? rw : rs;
For i In {1:NumWires}
pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
p1 = newp; Point(p1) = {WX~{i}+R, WY~{i} , 0, lc1};
p2 = newp; Point(p2) = {WX~{i} , WY~{i}+R, 0, lc1};
p3 = newp; Point(p3) = {WX~{i}-R, WY~{i} , 0, lc1};
p4 = newp; Point(p4) = {WX~{i} , WY~{i}-R, 0, lc1};
c1 = newl; Circle(c1) = {p1,pC,p2};
c2 = newl; Circle(c2) = {p2,pC,p3};
c3 = newl; Circle(c3) = {p3,pC,p4};
c4 = newl; Circle(c4) = {p4,pC,p1};
l1 = newl; Line(l1) = {pC,p1};
l2 = newl; Line(l2) = {pC,p2};
l3 = newl; Line(l3) = {pC,p3};
l4 = newl; Line(l4) = {pC,p4};
ll1 = newll; Line Loop(ll1) = {l1,c1,-l2};
s1 = news; Plane Surface(s1) = {ll1};
ll2 = newll; Line Loop(ll2) = {l2,c2,-l3};
s2 = news; Plane Surface(s2) = {ll2};
ll3 = newll; Line Loop(ll3) = {l3,c3,-l4};
s3 = news; Plane Surface(s3) = {ll3};
ll4 = newll; Line Loop(ll4) = {l4,c4,-l1};
s4 = news; Plane Surface(s4) = {ll4};
ll5 = newll;
Line Loop(ll5) = {c1,c2,c3,c4};
llWires[] += { ll5 };
surfWires~{i-1}[] = {s1,s2,s3,s4};
surfWires[] += surfWires~{i-1}[];
centerWires[] += {pC};
EndFor
Else
// then geometry only contains line conductors
For i In {1:NumWires}
pC = newp; Point(pC) = {WX~{i}, WY~{i}, 0, lc1};
centerWires[] += {pC};
EndFor
EndIf
// create air box around wires
BoundingBox; // recompute model bounding box
cx = (General.MinX + General.MaxX) / 2;
cy = (General.MinY + General.MaxY) / 2;
cz = (General.MinZ + General.MaxZ) / 2;
// lx = 2*inf + General.MaxX - General.MinX;
// ly = 2*inf + General.MaxY - General.MinZ;
lx = 2*box;
ly = 2*box;
p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, 0, lc2};
p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, 0, lc2};
p3 = newp; Point (p3) = {cx+lx/2, cy+ly/2, 0, lc2};
p4 = newp; Point (p4) = {cx-lx/2, cy+ly/2, 0, lc2};
l1 = newl; Line(l1) = {p1, p2};
l2 = newl; Line(l2) = {p2, p3};
l3 = newl; Line(l3) = {p3, p4};
l4 = newl; Line(l4) = {p4, p1};
ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
Skin[] = {};
If( !Flag_Thin ) // Wires with finite radius
For i In {1:NumWires}
e[] = Extrude {0,0,Lz}
{ Surface{ surfWires~{i-1}[] } ; /*Layers{1}; Recombine;*/ } ;
Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[ {1,6,11,16} ] };
Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires~{i-1}[] };
Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[ {0,5,10,15} ] };
Skin[] += e[ {3,8,13,18} ];
EndFor
s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
Physical Volume ("AIR", 1) = { e[1] };
ElseIf( Flag_Corr == 2 ) // regular sleeve
Wires[]= {};
For i In {1:NumWires}
e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; /*Layers{1};*/ } ;
Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
e[] = Extrude {0,0,Lz} { Surface{ surfWires~{i-1}[] } ;
/*Layers{1}; Recombine;*/ } ;
Wires() += e[ {1,6,11,16} ];
Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[ {1,6,11,16} ] };
Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires~{i-1}[] };
Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[ {0,5,10,15} ] };
EndFor
s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
Physical Volume ("AIR", 1) = { e[1], Wires[] };
Else
s1 = news; Plane Surface(s1) = {ll1};
e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
s2 = e[0];
v1 = e[1];
For i In {1:NumWires}
e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; /*Layers{1};*/ } ;
//Printf("e=", e[]);
Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
Physical Volume ("AIR", 1) = { v1 };
// embed line conductors in mesh
Point { centerWires[i-1] } In Surface { s1 };
Point { e[0] } In Surface { s2 };
Curve { e[1] } In Volume { v1 };
EndFor
EndIf
Physical Surface("INF", 2) = { CombinedBoundary{ Volume{ : }; } };
Physical Surface("SKIN", 3) = Skin[];
Electrodes[] = {};
If( !Flag_Thin ) // Wires with finite radius
For i In {1:NumWires}
Electrodes[] += {20+i, 30+i};
EndFor
EndIf
Physical Line("LINTREE", 4) = { Boundary { Physical Surface { Electrodes[] }; } };
Conductors3D/w3d.pro 0 → 100644
+ 1015
− 0
View file @ bb10fad3
Include "w3d_common.pro";
Struct S::SOLVER [ Enum, mumps, gmres, fgmres, bcgs ];
DefineConstant[
FileId = ""
DataDir = StrCat["data/" , FileId]
ElekinDir = StrCat["Elekin_Line_Const/" , FileId]
MagDynDir = StrCat["MagDyn_Line_Const/" , FileId]
FullWaveDir = StrCat["FullWave_Line_Const/", FileId]
type_solver = {
S::SOLVER.mumps,
Choices {
S::SOLVER.mumps ="mumps",
S::SOLVER.gmres="gmres",
S::SOLVER.fgmres= "fgmres",
S::SOLVER.bcgs="bcgs"
},
Name "Solver/type_solver"
}
];
DefineConstant[
R_ = {"Dynamic", Name "GetDP/1ResolutionChoices", Visible 0},
C_ = {"-sol -pos", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
ResDir = "Res/"
];
Group{
// Geometrical regions
AIR = Region[ 1 ];
INF = Region[ 2 ];
SKIN = Region[ 3 ];
LINTREE = Region[ 4 ]; // not used
VWIRES = Region[ {} ];
LWIRES = Region[ {} ];
ANODES = Region[ {} ];
CATHODES = Region[ {} ];
For i In {1:NumWires}
LWIRE~{i} = Region[ (50+i) ];
LWIRES += Region[ LWIRE~{i} ];
If( !Flag_Thin )
VWIRE~{i} = Region[ (10+i) ];
VWIRES += Region[ VWIRE~{i} ];
SWIRE~{i} = Region[ {} ];
ANODE~{i} = Region[ (20+i) ];
CATHODE~{i} = Region[ (30+i) ];
Else
VWIRE~{i} = Region[ {} ];
ANODE~{i} = Region[ (60+i) ];
CATHODE~{i} = Region[ (70+i) ];
SWIRE~{i} = ElementsOf[ AIR, OnOneSideOf LWIRE~{i} ];
EndIf
ANODES += Region[ ANODE~{i} ];
CATHODES += Region[ CATHODE~{i} ];
EndFor
// Abstract regions
// Cette formulation peut contenir des conducteurs volumiques Vol_C
// et liniques lVol_C.
If( !Flag_Thin )
Vol_C = Region[ VWIRES ];
Vol_CC = Region[ AIR ];
//lVol_C = Region[ {} ];
Else
Vol_C = Region[ {LWIRES} ];
Vol_CC = Region[ {VWIRES, AIR} ];
//lVol_C = Region[ LWIRES ];
EndIf
Vol_nu = Region[ {Vol_C, Vol_CC} ];
Sur_Dirichlet_a = Region[ INF ];
Electrodes = Region[ {ANODES, CATHODES} ];
Vol_Tree = Region[ { Vol_nu } ];
Sur_Tree = Region[ { Sur_Dirichlet_a /*, SKIN*/ } ];
// If( Flag_Corr!= 0 )
// Lin_Tree = Region[ { LWIRES /*, LINTREE*/ } ];
// Else
Lin_Tree = Region[ {} ];
Dom_Hcurl_a = Region[ Vol_nu ];
Dom_Hgrad_u = Region[ Vol_C ];
//Dom_Hregion_i = Region[ Vol_nu ];
// additional Groups for the local correction method
//Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf LWIRES ];
Dom_Hthin_a = Region[ Vol_nu ];
// For integration on the volume elements (excluding line elements)
// of the sleeve
Vol_Sleeve = ElementsOf[ {Vol_CC}, OnOneSideOf LWIRES ];
}
Function{
omega = 2*Pi*Freq;
mu0 = Pi*4e-7;
mur = 1.;
sigma = 5.96e7;
i[] = Complex[0,1];
mu[] = mu0;
nu[] = 1/mu[];
sigma[] = sigma;
skin_depth = Sqrt[ 2/(omega*sigma*mu0*mur) ];
tau[] = Complex[-1,1]/skin_depth*rw;
J0[] = JnComplex[0,$1];
J1[] = JnComplex[1,$1];
dJ1[] = JnComplex[0,$1]-JnComplex[1,$1]/$1;
For i In {1:NumWires}
// Current density
J[ Region[ {VWIRE~{i}, LWIRE~{i}} ] ] =
Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
// local radial coordinate with origin on wire i
R~{i}[] = Hypot[ X[]-WX~{i} , Y[]-WY~{i} ];
EndFor
// shape function of the current, wI[r]
wI[] = tau[] * J0[tau[]*$1/rw] / J1[tau[]] /(2*A_c);
//radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
Analytic_A[] = // per Amp
(($1>rw) ? mu0/(2*Pi) * Log[rs/$1] :
mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) );
AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
(($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
Analytic_B[] = // per Amp
(($1>rw) ? mu0/(2*Pi*$1) :
mu0*mur/(2*Pi) * J1[tau[]*$1/rw] / J1[tau[]] );
// Impedance of thin wire p.u. length
R_DC = 1./(sigma*A_c);
Analytic_R[] = Re[ wI[rw]/sigma ];
Analytic_L[] = mu0/(2*Pi) * Log[rs/rw] + Re[ wI[rw]/(i[]*omega*sigma) ];
If( NumWires == 1 )
Exact_B[] = Analytic_B[R_1[]] ;
Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] : 0 ) ;
EndIf
If ( NumWires == 2 )
Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] ;
Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] :
((R_2[]<rs) ? Analytic_B[R_2[]] : 0 ));
EndIf
If ( NumWires == 3 )
Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] + Analytic_B[R_3[]] ;
Correction_A[] = ((R_1[]<rs) ? Analytic_B[R_1[]] :
((R_2[]<rs) ? Analytic_B[R_2[]] :
((R_3[]<rs) ? Analytic_B[R_3[]] : 0 )));
EndIf
}
Jacobian {
{ Name Vol ;
Case {
{ Region LWIRES ; Jacobian Lin ; Coefficient A_c; }
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 3 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Prism ; NumberOfPoints 21 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Pyramid ; NumberOfPoints 8 ; }
}
}
}
}
}
Constraint{
{ Name Impose_U ;
Case {
For i In {1:NumWires}
If( Flag_Form == 1) // massive
{ Region ANODE~{i} ; Value 0 ; }
EndIf
If( Flag_U )
If( Flag_Form == 1) // massive
{ Region CATHODE~{i} ; Value -R_DC*Lz*WI~{i} ; }
Else // stranded
{ Region VWIRE~{i} ; Value -R_DC*Lz*WI~{i} ; }
EndIf
EndIf
EndFor
}
}
{ Name Impose_I ;
Case {
For i In {1:NumWires}
If( !Flag_U )
If( Flag_Form == 1) // massive
{ Region ANODE~{i} ; Value WI~{i} ; }
{ Region CATHODE~{i} ; Value WI~{i} ; }
Else // stranded
{ Region VWIRE~{i} ; Value WI~{i} ; }
{ Region LWIRE~{i} ; Value WI~{i} ; }
EndIf
EndIf
EndFor
}
}
{ Name Hcurl_a_3D_ac; Type Assign;
Case {
{ Region Region[ Sur_Dirichlet_a ]; Value 0.; }
}
}
{ Name GaugeCondition_a; Type Assign;
Case {
{ Region Vol_Tree; SubRegion Region[ { Sur_Tree, Lin_Tree } ];
Value 0.; }
}
}
}
FunctionSpace {
// Magnetic vector potential
{ Name Hcurl_a_3D; Type Form1;
BasisFunction {
{ Name sc; NameOfCoef ac; Function BF_Edge;
Support Dom_Hcurl_a; Entity EdgesOf[All]; }
}
Constraint {
{ NameOfCoef ac;
EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
{ NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
NameOfConstraint GaugeCondition_a ; }
}
}
// Electric scalar potential for massive conductors
{ Name Hgrad_u_3D; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Node;
Support Dom_Hgrad_u; Entity NodesOf[ Dom_Hgrad_u, Not Electrodes]; }
{ Name sn2; NameOfCoef un2; Function BF_GroupOfNodes;
Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }
}
GlobalQuantity {
{ Name U; Type AliasOf ; NameOfCoef un2; }
{ Name I; Type AssociatedWith; NameOfCoef un2; }
}
Constraint {
{ NameOfCoef U; EntityType GroupsOfNodesOf; NameOfConstraint Impose_U; }
{ NameOfCoef I; EntityType GroupsOfNodesOf; NameOfConstraint Impose_I; }
}
}
{ Name Hregion_i_3D; Type Vector;
BasisFunction {
{ Name sr; NameOfCoef ir; Function BF_RegionZ;
Support Dom_Hthin_a; Entity Vol_C; }
}
GlobalQuantity {
{ Name Is; Type AliasOf ; NameOfCoef ir; }
{ Name Us; Type AssociatedWith; NameOfCoef ir; }
}
Constraint {
{ NameOfCoef Us; EntityType Region; NameOfConstraint Impose_U; }
{ NameOfCoef Is; EntityType Region; NameOfConstraint Impose_I; }
}
}
// For the local correction method
{ Name Hcurl_acorr_3D; Type Form1;
BasisFunction {
{ Name sw; NameOfCoef aw; Function BF_Edge;
Support Vol_nu; Entity EdgesOf[ Vol_nu, ConnectedTo LWIRES ]; }
}
Constraint {
// { NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
// NameOfConstraint GaugeCondition_a ; }
}
}
}
Formulation {
{ Name MagnetoDynamics; Type FemEquation;
If( Flag_Form == 1) // Massive conductor
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_3D; }
{ Name v; Type Local; NameOfSpace Hgrad_u_3D; }
{ Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
{ Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
}
Equation {
Integral { [ nu[] * Dof{d a} , {d a} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
}
Else // stranded conductor
If( Flag_Corr==0 ) // naive thin wire model
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_3D; }
{ Name i; Type Local; NameOfSpace Hregion_i_3D; }
{ Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
{ Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
}
Equation {
Integral { [ nu[] * Dof{d a} , {d a} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
// Integral { [ -J[] , {a} ];
// In Vol_C; Jacobian Vol; Integration I1; }
// Us is the voltage drop
// grad v = -Dt a - j/sigma
Integral { DtDof [ Dof{a} , {i} ];
In Vol_C; Jacobian Vol; Integration I1;}
Integral { [ Dof{i}/(sigma[]*A_c) , {i} ];
In Vol_C; Jacobian Vol; Integration I1;}
GlobalTerm { [ Dof{Us}*A_c , {Is} ]; In Vol_C; }
}
Else // use local correction method
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_3D; }
{ Name as; Type Local; NameOfSpace Hcurl_acorr_3D; }
{ Name i; Type Local; NameOfSpace Hregion_i_3D; }
{ Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
{ Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
}
Equation {
/*
Integral { [ nu[] * Dof{d a} , {d a} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
*/
Integral { [ nu[] * Dof{d as} , {d as} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {as} ];
In Vol_C; Jacobian Vol; Integration I1; }
/*
GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
In Vol_C;}
Integral { DtDof [ Dof{a}/A_c , {i} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof [ -Dof{as}/A_c , {i} ];
In Vol_C; Jacobian Vol; Integration I1;}
GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
In Vol_C;}
GlobalTerm { [ Dof{Us} , {Is} ];
In Vol_C; }
*/
}
EndIf
EndIf
}
}
Resolution {
{ Name Dynamic;
System {
{ Name S; NameOfFormulation MagnetoDynamics;
Type ComplexValue; Frequency Freq; }
}
Operation {
CreateDir[DataDir];
Generate[S]; Solve[S]; SaveSolution[S];
If( Flag_Maps == 1 )
// DeleteFile["U.dat"];
// DeleteFile["I.dat"];
// DeleteFile["Z.dat"];
PostOperation[map];
EndIf
PostOperation[integaz];
PostOperation[cut];
For i In {1:NumWires}
Print[ {i, rw, rs, A_c, skin_depth, R_DC, WI~{i}/A_c, mu0*WI~{i}/(2*Pi*rw)},
Format "WIRE %2g: rw = %7.3e rs = %7.3e A_c = %7.3e delta = %7.3e R_DC = %7.3e J = %7.3e Bmax = %7.3e"] ;
// Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}},
// Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
EndFor
}
}
}
PostProcessing {
{ Name MagnetoDynamics; NameOfFormulation MagnetoDynamics;
PostQuantity {
{ Name a;
Value { Local { [ {a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name b;
Value { Local { [ {d a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name by;
Value { Local { [ Re[Cart2Pol[CompY [ {d a}]]]];
In Dom_Hcurl_a; Jacobian Vol; }}}
If( Flag_Form == 1 ) // massive
{ Name J;
Value { Local { [ -sigma[] * ( Dt[{a}] + {d v} ) ];
In Vol_C; Jacobian Vol; }}}
{ Name U;
Value { Term { [ {U} ]; In Electrodes; }}}
{ Name I;
Value { Term { [ {I} ]; In Electrodes; }}}
{ Name R;
Value { Term { [ Re[ -{U}/{I}/Lz ] ]; In Electrodes; }}}
{ Name L; // actually total flux/I
Value { Term { [ Im[ -{U}/{I}/omega/Lz ] ]; In Electrodes; }}}
Else // stranded
{ Name J;
Value { Local { [ {i}/A_c ];
In Vol_C; Jacobian Vol; }}}
{ Name U;
Value { Term { [ {Us} ]; In Vol_C; }}}
{ Name I;
Value { Term { [ {Is} ]; In Vol_C; }}}
{ Name R;
Value { Term { [ Re[ -{Us}/{Is}/Lz ] ]; In Vol_C; }}}
{ Name L; // actually total flux/I
Value { Term { [ Im[ -{Us}/{Is}/omega/Lz ] ]; In Vol_C; }}}
EndIf
If( Flag_Corr )
{ Name bs;
Value { Local { [ {d as}]; In Dom_Hthin_a; Jacobian Vol; }}}
{ Name bcorry;
Value { Local { [ Re[Cart2Pol[CompY [{d a}-{d as}]]] + Correction_B[] ];
In Dom_Hthin_a; Jacobian Vol; }}}
EndIf
}
}
}
PostOperation map UsingPost MagnetoDynamics {
Print[ J, OnElementsOf Vol_C, File "j.pos"];
Print[ b, OnElementsOf Vol_nu, File "b.pos"];
If( Flag_Corr )
Print[ bs, OnElementsOf Vol_nu, File "bs.pos"];
PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo LWIRES ],
In Vol_nu, File "group.pos"];
EndIf
}
PostOperation integaz UsingPost MagnetoDynamics {
For i In {1:NumWires}
If( Flag_Form == 1 ) // massive
Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
Print[ I, OnRegion CATHODES, File > "I.dat", Format Table,
SendToServer Sprintf["Results/2Current/Wire %g",i]];
Print[ R, OnRegion CATHODES, File > "Z.dat", Format Table,
SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
Print[ L, OnRegion CATHODES, File > "Z.dat", Format Table,
SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
Else
Print[ U, OnRegion Vol_C, File > "U.dat", Format Table,
SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
Print[ I, OnRegion Vol_C, File > "I.dat", Format Table,
SendToServer Sprintf["Results/2Current/Wire %g",i]];
Print[ R, OnRegion Vol_C, File > "Z.dat", Format Table,
SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
Print[ L, OnRegion Vol_C, File > "Z.dat", Format Table,
SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
EndIf
EndFor
}
NbPoints = 400;
Xmax = 0.15*box;
PostOperation cut UsingPost MagnetoDynamics {
Print [ by, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "by.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut by';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
If( Flag_Corr )
Print [ bcorry, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "by.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut by';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
EndIf
}
/*
PostProcessing {
{ Name MagnetoDynamics_corr; NameOfFormulation MagnetoDynamics;
PostQuantity {
{ Name a;
Value { Local { [ {a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name b;
Value { Local { [ {d a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name by;
Value { Local { [ Cart2Pol[CompY [ {d a}]]]; In Dom_Hcurl_a; Jacobian Vol; }}}
If( Flag_Form <= 1 )
{ Name U;
Value { Term { [ {U} ]; In CATHODES; }}}
{ Name flux;
Value { Term { [ Im[ {U}/omega ] ]; In CATHODES; }}}
Else
{ Name bc;
Value { Local { [ {d ac}]; In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name bw;
Value { Local { [ {d aw}]; In Dom_Hthin_a; Jacobian Vol; }}}
{ Name byw;
Value { Local { [ Cart2Pol[CompY [ {d aw}]]]; In Dom_Hthin_a; Jacobian Vol; }}}
{ Name U;
Value { Term { [ {Us} ]; In CATHODES; }}}
{ Name flux;
Value { Term { [ Im[ {Us}/omega ] ]; In CATHODES; }}}
EndIf
}
}
}
PostOperation map_corr UsingPost MagnetoDynamics {
Print[ b, OnElementsOf Vol_nu, File "b.pos"];
Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g",i]];
Print[ flux, OnRegion CATHODES, File > "L.dat", Format Table,
SendToServer Sprintf["}Inductance/Wire %g",i]];
If( Flag_Form >= 2 )
PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo lVol_C ],
In Vol_nu, File "group.pos"];
Print [bw, OnElementsOf Vol_nu, File "b.pos"];
EndIf
}
FunctionSpace {
{ Name Hcurl_ai_3D; Type Form1;
BasisFunction {
{ Name sc; NameOfCoef ac; Function BF_Edge;
Support Dom_Hcurl_a; Entity EdgesOf[All]; }
{ Name sw; NameOfCoef aw; Function BF_Edge;
Support Dom_Hthin_a; Entity EdgesOf[ Vol_nu, ConnectedTo lVol_C ]; }
}
SubSpace {
{ Name ac ; NameOfBasisFunction { sc } ; }
{ Name aw ; NameOfBasisFunction { sw } ; }
}
Constraint {
{ NameOfCoef ac;
EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
{ NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
NameOfConstraint GaugeCondition_a ; }
{ NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
NameOfConstraint GaugeCondition_a ; }
}
}
}
Formulation {
{ Name Corrected; Type FemEquation;
If( Flag_Form == 1) // A-v formulation
// Conventional A-v formulation
// Works with FE discretized true wires of finite radius (Flag_Thin==0)
// and idealized thin wires (Flag_Thin!=0)
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_3D; }
{ Name v; Type Local; NameOfSpace Hgrad_u_3D; }
{ Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
{ Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
}
Equation {
Integral { [ nu[] * Dof{d a} , {d a} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {d v} ];
In Vol_C; Jacobian Vol; Integration I1; }
//GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {a} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{d v} , {d v} ];
In lVol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
}
Else // A-I formulation
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_ai_3D; }
{ Name ac; Type Local; NameOfSpace Hcurl_ai_3D[ac]; }
{ Name aw; Type Local; NameOfSpace Hcurl_ai_3D[aw]; }
{ Name i; Type Local; NameOfSpace Hregion_i_3D; }
{ Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
{ Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
}
Equation {
Integral { [ nu[] * Dof{d ac} , {d ac} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ nu[] * Dof{d aw} , {d aw} ];
In Vol_nu; Jacobian Vol; Integration I1; }
If( 0 ) // impose current density
Integral { [ -J[], {ac} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -J[] , {aw} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -J[] , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Else
Integral { [ -Dof{i}/A_c , {ac} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {aw} ];
In lVol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [ Dof{Us} , {Is} ]; In lVol_C; }
GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; In lVol_C;}
Integral { DtDof [ Dof{ac}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1; }
If( Flag_Form != 2)
Integral { DtDof [ -Dof{aw}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1;}
GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
In lVol_C;}
EndIf
EndIf
}
EndIf
}
}
*/
// ################################
/*
NbPoints = 1000;
Xmax = 0.15*box;
PostOperation cut UsingPost PostProcessings {
// cut of by displayed in GUI
If( Flag_Thin == 0 )
Print [ by, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut |by|';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
Else
Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut az not corrected';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
Format Gmsh, File "Cut_az.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].Name = 'cut az corrected';",
"View[l].Axes = 3;",
"View[l].LineWidth = 3;",
"View[l].Type = 2;"],
File "tmp.geo", LastTimeStepOnly];
EndIf
// cuts in txt format for Gnuplot
Print [ exact, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_analytic.txt" ];
If( Flag_Thin == 0 )
Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_fullfem.txt" ];
Else
Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_thinfem.txt" ];
Print [ acz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_ac.txt" ];
Print [ awz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_aw.txt" ];
Print [ acwz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_acw.txt" ];
Print [ corr, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
Format SimpleTable, File "Cut_correction.txt" ];
EndIf
}
PostProcessing {
{ Name PostProcessings; NameOfFormulation MagnetoDynamics;
PostQuantity {
If( Flag_3D)
{ Name a;
Value { Local { [ {a}];
In Dom_Hcurl_a; Jacobian Vol; }}}
// { Name ac;
// Value { Local { [ {ac} ];
// In Dom_Hcurl_a; Jacobian Vol; }}}
// { Name aw;
// Value { Local { [ {aw} ];
// In Dom_Hthin_a; Jacobian Vol; }}}
{ Name b;
Value { Local { [ {d a}];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name U;
Value { Term { [ {U} ]; In CATHODES; }}}
{ Name L;
Value { Term { [ Im[ {U}/{I}/omega ] ]; In CATHODES; }}}
Else ## 2D
{ Name az;
Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name acz;
Value { Local { [ CompZ[ {ac}] ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name awz;
Value { Local { [ CompZ[ {aw}] ];
In Dom_Hthin_a; Jacobian Vol; }}}
{ Name acwz;
Value { Local { [ CompZ[ {ac} - {aw} ] ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name exact;
Value { Local { [ Exact_A[] ] ;
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name corr;
Value { Local { [ Correction_A[] ] ;
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name bcw;
Value { Local { [ {d ac} - {d aw} ];
In Dom_Hcurl_a; Jacobian Vol; }}}
{ Name j;
Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ];
In Vol_C; Jacobian Vol; }}}
{ Name un;
Value{ Local { [ 1 ] ;
In Vol_nu; Integration I1 ; Jacobian Vol; }}}
{ Name surf;
Value{ Integral{ [ 1 ] ;
In Vol_nu; Integration I1 ; Jacobian Vol; }}}
If( Flag_Thin == 0 )
{ Name flux;
Value { Integral { [ CompZ[ {a} ] / A_c ];
In Vol_C; Integration I1 ; Jacobian Vol; }}}
{ Name JouleLosses;
Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
// { Name U;
// Value { Term { [ Cart2Pol [ {U} ] ]; In Vol_C; }}}
{ Name U;
Value { Term { [ Im [ {U}/omega ] ]; In Vol_C; }}}
Else
{ Name flux; // naive flux
Value { Integral { [ CompZ[ {ac} ] ];
In lVol_C; Integration I1 ; Jacobian Vol; }}}
{ Name JouleLosses;
Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
// { Name U;
// Value { Term { [ Cart2Pol [ {Us} ] ] ; In lVol_C; }}}
{ Name U;
Value { Term { [ Im [ {Us}/omega ] ] ; In lVol_C; }}}
For i In {1:NumWires}
{ Name area~{i};
Value{ Integral{ [ 1 ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
{ Name intbx~{i}; // integral of b, 1st step for computing b average
Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
{ Name intby~{i}; // integral of b, 1st step for computing b average
Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ;
In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
EndFor
EndIf
EndIf
}
}
}
PostOperation map UsingPost PostProcessings {
Print [az, OnElementsOf Vol_nu, File "az.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly];
If(Flag_Thin != 0)
Print [bcw, OnElementsOf Vol_Sleeve, File "b.pos"];
Print [acwz, OnElementsOf Vol_nu, File "acz.pos"];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
"View[l].RaiseZ = 15000;",
"View[l].NormalRaise = 0;"],
File "tmp.geo", LastTimeStepOnly] ;
EndIf
}
PostOperation integaz UsingPost PostProcessings {
For i In {1:NumWires}
If( Flag_Thin == 0 )
Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
Format Table, File > "joule_full.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
EndIf
If( Flag_Thin == 1 )
Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/2Raw",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_raw.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/2Raw",i]];
EndIf
If( Flag_Thin == 2 )
Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_reg.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
EndIf
If( Flag_Thin == 3 )
Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
Print[ JouleLosses, OnRegion LWIRE~{i},
Format Table, File > "joule_naive.dat",
SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
EndIf
If( Flag_Thin != 0 )
Print[ area~{i}[SWIRE~{i}], OnGlobal,
Format Table, File >"zzz",
StoreInVariable $sarea~{i}];
Print[ intby~{i}[SWIRE~{i}], OnGlobal,
Format Table, File > "zzz",
StoreInVariable $intby~{i}];
EndIf
EndFor
}
*/
/* 2D case
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_2D; }
{ Name ac; Type Local; NameOfSpace Hcurl_a_2D[ac]; }
{ Name aw; Type Local; NameOfSpace Hcurl_a_2D[aw]; }
{ Name dv; Type Local; NameOfSpace Hgrad_u_2D; }
{ Name U; Type Global; NameOfSpace Hgrad_u_2D [U]; }
{ Name I; Type Global; NameOfSpace Hgrad_u_2D [I]; }
{ Name i; Type Local; NameOfSpace Hregion_i_2D; }
{ Name Is; Type Global; NameOfSpace Hregion_i_2D [Is]; }
{ Name Us; Type Global; NameOfSpace Hregion_i_2D [Us]; }
}
Equation {
Integral { [ nu[] * Dof{d ac} , {d ac} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [ -Dof{i}/A_c , {ac} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ nu[] * Dof{d aw} , {d aw} ];
In Vol_nu; Jacobian Vol; Integration I1; }
Integral { [-Dof{i}/A_c , {aw} ];
In lVol_C; Jacobian Vol; Integration I1; }
Integral { [ -J[] , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{dv} , {a} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { DtDof[ sigma[] * Dof{a} , {dv} ];
In Vol_C; Jacobian Vol; Integration I1; }
Integral { [ sigma[] * Dof{dv} , {dv} ];
In Vol_C; Jacobian Vol; Integration I1; }
GlobalTerm { [ Dof{I} , {U} ]; In Vol_C; }
GlobalTerm { [ Dof{Us} , {Is} ];
In lVol_C; }
GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
In lVol_C;}
Integral { DtDof [ Dof{ac}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1; }
If( Flag_Thin != 3)
Integral { DtDof [ -Dof{aw}/A_c , {i} ];
In lVol_C; Jacobian Vol; Integration I1;}
GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
In lVol_C;}
EndIf
}
*/
Conductors3D/w3d_common.pro 0 → 100644
+ 57
− 0
View file @ bb10fad3
mm = 1e-3;
deg = 180/Pi;
/*
Expected inductance of a rectilinear 1mm radius wire: 4 nH/cm
*/
DefineConstant[
NumWires = {1, Name "Parameters/0Number of wires",
Choices {1="1", 2="2", 3="3"}, Visible 1}
Flag_Form = {2, Name "Parameters/1Conductor type", Visible 1,
Choices {1="Massive", 2="Stranded"}}
Flag_Thin = {1, Name "Parameters/2Use thin wires", Choices {0,1}, Visible 1}
Flag_Corr = {1, Name "Parameters/3Correction method",
Visible (Flag_Form==2) && Flag_Thin,
Choices {0="None", 1="Raw sleeve", 2="regular sleeve"}}
Flag_U = {0, Name "Parameters/4Impose U", Choices {0,1}, Visible 1}
Flag_Maps = {1, Name "Input/2Display maps", Choices {0,1}, Visible 1}
LogFreq = {4, Min 1, Max 8, Step 0.25, Name "Input/4Log of frequency"}
Freq = 10^LogFreq
WireRadius = {1, Name "Parameters/5Wire radius [mm]"}
MeshSizeWire = {1, Name "Parameters/6Mesh size in wire [mm]",
Visible !Flag_Thin}
MeshSizeThinWire = {1, Name "Parameters/7Mesh size on thin wire [mm]",
Visible Flag_Thin}
SleeveRadius = {2, Name "Parameters/8Sleeve radius [mm]"}
box = {100*mm, Name "Parameters/9box half-width [m]"}
rw = WireRadius*mm
rs = SleeveRadius*mm
A_c = Pi*rw^2 // wire cross section
rout = box
Lz = 10*mm
/* 'WireRadius' is the actual radius of the wire.
'MeshSizeWire' is the imposed mesh size at the nodes of the wire,
and thus also the practical approximative value of the radius of the sleeve.
The 'rs' and 'rw' parameters are assumed identical for all wires.
*/
];
Xlocs() = { 0, 5*mm, -5*mm};
Ylocs() = { 0, 0, 0};
For i In {1:NumWires}
DefineConstant[
WX~{i} = { Xlocs(i-1), // place wires symmetrically around x=0
Name Sprintf("Parameters/Wire %g/0X position [m]", i), Closed }
WY~{i} = { Ylocs(i-1), Name Sprintf("Parameters/Wire %g/0Y position [m]", i) }
WI~{i} = { 1.23456, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
WP~{i} = { 0*NumWires*(i-1),
Name Sprintf("Parameters/Wire %g/3Phase [deg]", i) }
];
EndFor
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment