[Getdp] OneTooth Example with non linear material
Matthias Lang
lang_matthias at gmx.de
Tue Oct 22 23:55:37 CEST 2013
Hello ererybody,
I took the 3D magnetostatic example OneTooth.pro
and included nonlinear material. (found in Material.pro)
But unfortunatly I get the following Error Message:
.
.
Info : GetDP - E n d P r e - P r o c e s s i n g
Info : GetDP - P r o c e s s i n g . . .
Info : GetDP - IterativeLoop ...
Info : GetDP - Non Linear Iteration Relaxation 1
Info : GetDP - GenerateJac[Sys_Mag]
Error : GetDP - Function _Exp_Fct_17 called with too few arguments
Info : Done running 'GetDP'
Where is the problem? - inside or in front of the computer?
Regards
Matthias Lang
-------------- next part --------------
/* -------------------------------------------------------------------
File "MagSta_a_3D_nl.pro"
Magnetostatics - Magnetic vector potential a formulation (3D)
-------------------------------------------------------------------
I N P U T
---------
GlobalGroup : (Extension Mag is for Magnetic problem)
-----------
Domain_Mag Whole magnetic domain
DomainS_Mag Inductor regions (Source)
Function :
--------
nu[] Magnetic reluctivity
Constraint :
----------
*/
Group {
DefineGroup[ Domain_Mag, DomainS_Mag, Domain_NL];
}
Function {
DefineFunction[ nu, dhdb_NL ];
DefineConstant[
Nb_max_iter = 30,
relaxation_factor = 1,
stop_criterion = 1e-5,
NbT = 1,
Steel_1020_dhdb_NL
];
}
FunctionSpace {
// Magnetic vector potential a (b = curl a)
{ Name Hcurl_a_Mag_3D; Type Form1;
BasisFunction {
// a = a s
// e e
{ Name se; NameOfCoef ae; Function BF_Edge;
Support Domain_Mag; Entity EdgesOf[ All ]; }
}
Constraint {
{ NameOfCoef ae; EntityType EdgesOf;
NameOfConstraint MagneticVectorPotential_3D; }
{ NameOfCoef ae; EntityType EdgesOfTreeIn; EntitySubType StartingOn;
NameOfConstraint Gauge; }
}
}
}
Formulation {
{ Name Magnetostatics_a_3D; Type FemEquation;
Quantity {
{ Name a ; Type Local; NameOfSpace Hcurl_a_Mag_3D; }
}
Equation {
Galerkin { [ nu[] * Dof{Curl a} , {Curl a} ]; In Domain_Mag;
Jacobian Vol; Integration CurlCurl; }
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration CurlCurl ; }
Galerkin { [ - js[] , {a} ]; In DomainS_Mag;
Jacobian Vol; Integration CurlCurl; }
}
}
}
Resolution {
{ Name MagSta_a_3D;
System {
{ Name Sys_Mag; NameOfFormulation Magnetostatics_a_3D; }
}
Operation {
// Generate[Sys_Mag]; Solve[Sys_Mag]; //linear case
IterativeLoop[Nb_max_iter, stop_criterion, relaxation_factor]{
GenerateJac[Sys_Mag] ; SolveJac[Sys_Mag] ; }
SaveSolution[Sys_Mag];
}
}
}
PostProcessing {
{ Name MagSta_a_3D; NameOfFormulation Magnetostatics_a_3D;
Quantity {
{ Name a;
Value {
Local { [ {a} ]; In Domain_Mag; Jacobian Vol; }
}
}
{ Name a_norm ;
Value {
Local { [ Norm[{a}] ] ; In Domain_Mag ; Jacobian Vol ; }
}
}
{ Name b;
Value {
Local { [ {d a} ]; In Domain_Mag; Jacobian Vol; }
}
}
{ Name h;
Value {
Local { [ nu[] * {d a} ]; In Domain_Mag; Jacobian Vol; }
}
}
{ Name js;
Value {
Local { [js[]] ; In DomainS_Mag; Jacobian Vol; }
}
}
}
}
}
-------------- next part --------------
Function{
/*--------------------------------------------------------------------------------------*/
// nu = 100. + 10. * exp ( 1.8 * b * b )
// analytical
nu_1a[] = 100. + 10. * Exp[1.8*SquNorm[$1]] ;
dnudb2_1a[] = 18. * Exp[1.8*SquNorm[$1]] ;
h_1a[] = nu_1a[$1]*$1 ;
dhdb_1a[] = TensorDiag[1,1,1] * nu_1a[$1#1] + 2*dnudb2_1a[#1] * SquDyadicProduct[#1] ;
dhdb_1a_NL[] = 2*dnudb2_1a[$1#1] * SquDyadicProduct[#1] ;
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
// BH data EIcore
MatEIcore_h = {0.,
2.1827e+02,3.8272e+02,5.2036e+02,7.1167e+02,
8.4921e+02,1.0405e+03,1.4741e+03,2.0409e+03,
3.0109e+03,5.0572e+03,7.5335e+03,1.0037e+04,
1.2486e+04,1.5015e+04,1.7464e+04,2.0021e+04,
2.2040e+04,2.5000e+04} ;
MatEIcore_b = {0.,
2.0329e-01,4.0287e-01,6.0986e-01,8.0575e-01,
1.0053e+00,1.1975e+00,1.4008e+00,1.5154e+00,
1.5967e+00,1.6706e+00,1.7076e+00,1.7335e+00,
1.7446e+00,1.7593e+00,1.7667e+00,1.7815e+00,
1.7889e+00,1.7963e+00} ;
MatEIcore_b2 = List[MatEIcore_b]^2 ;
MatEIcore_nu = List[MatEIcore_h]/List[MatEIcore_b] ;
MatEIcore_nu(0) = MatEIcore_nu(1);
MatEIcore_nu_b2 = ListAlt[MatEIcore_b2, MatEIcore_nu] ;
//nu_EIcore[] = InterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ; // used later in *.pro file
MatEIcore_interp_nu[] = InterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ; // used later in *.pro file
//dnudb2_EIcore[] = dInterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ;
MatEIcore_interp_dnudb2[] = dInterpolationLinear[SquNorm[$1]]{List[MatEIcore_nu_b2]} ;
//h_EIcore[] = nu_EIcore[$1] * $1 ;
MatEIcore_interp_h[] = MatEIcore_interp_nu[$1] * $1 ;
//dhdb_EIcore[] = TensorDiag[1,1,1]*nu_EIcore[$1#1] + 2*dnudb2_EIcore[#1] * SquDyadicProduct[#1] ;
MatEIcore_dhdb[] = TensorDiag[1,1,1]*MatEIcore_interp_nu[$1#1] + 2*MatEIcore_interp_dnudb2[#1] * SquDyadicProduct[#1] ;
//dhdb_EIcore_NL[] = 2*dnudb2_EIcore[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
MatEIcore_dhdb_NL[] = 2*MatEIcore_interp_dnudb2[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
// BH data for: 1020 Steel
Steel_1020_h = {0.000000,
79.577472, 100.182101, 126.121793, 158.777930,
199.889571, 251.646061, 316.803620, 398.832128,
502.099901, 632.106325, 795.774715, 1001.821011,
1261.217929, 1587.779301, 1998.895710, 2516.460605,
3168.036204, 3988.321282, 5020.999013, 6321.063250,
7957.747155, 10018.210114, 12612.179293, 15877.793010,
19988.957103, 25164.606052, 31680.362037, 39883.212823,
50209.990127, 63210.632497, 79577.471546, 100182.101136,
126121.792926, 158777.930096, 199889.571030, 251646.060522,
316803.620370};
Steel_1020_b = {0.000000,
0.166252, 0.208725, 0.261635, 0.327146,
0.407495, 0.504604, 0.619366, 0.750502,
0.893195, 1.038258, 1.173276, 1.286731,
1.373533, 1.437658, 1.489189, 1.538344,
1.591321, 1.649723, 1.711927, 1.774712,
1.834976, 1.891442, 1.945262, 1.998705,
2.052852, 2.106184, 2.155149, 2.196375,
2.229008, 2.255450, 2.280053, 2.307136,
2.339975, 2.381044, 2.432708, 2.497748,
2.579627};
Steel_1020_b2 = List[Steel_1020_b]^2 ;
Steel_1020_nu = List[Steel_1020_h]/List[Steel_1020_b] ;
Steel_1020_nu(0) = Steel_1020_nu(1);
Steel_1020_nu_b2 = ListAlt[Steel_1020_b2, Steel_1020_nu] ;
//nu_EIcore[] = InterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ; // used later in *.pro file
Steel_1020_interp_nu[] = InterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ; // used later in *.pro file
//dnudb2_EIcore[] = dInterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ;
Steel_1020_interp_dnudb2[] = dInterpolationLinear[SquNorm[$1]]{List[Steel_1020_nu_b2]} ;
//h_EIcore[] = nu_EIcore[$1] * $1 ;
Steel_1020_interp_h[] = Steel_1020_interp_nu[$1] * $1 ;
//dhdb_EIcore[] = TensorDiag[1,1,1]*nu_EIcore[$1#1] + 2*dnudb2_EIcore[#1] * SquDyadicProduct[#1] ;
Steel_1020_dhdb[] = TensorDiag[1,1,1]*Steel_1020_interp_nu[$1#1] + 2*Steel_1020_interp_dnudb2[#1] * SquDyadicProduct[#1] ;
//dhdb_EIcore_NL[] = 2*dnudb2_EIcore[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
Steel_1020_dhdb_NL[] = 2*Steel_1020_interp_dnudb2[$1] * SquDyadicProduct[$1] ; // used later in *.pro file
/*--------------------------------------------------------------------------------------*/
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: OneTooth.geo
Type: application/vnd.dynageo
Size: 7369 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/getdp/attachments/20131022/186a4a61/attachment.geo>
-------------- next part --------------
// Fichier "OneTooth.pro"
// Magnétostatique
Include "OneTooth-dat.pro";
Include "Materials.pro"; // nonlinear BH caracteristic of magnetic material
Group {
Air = #AIR;
Core = #CORE;
Ind = #COIL;
Plate = #PLATE;
SurfaceGe0 = #CUT_SYMMETRY ;
SurfaceGInf = #AIR_INF;
SkinCore = #SKIN_CORE;
SkinPlate = #SKIN_PLATE;
DomainS_Mag = Region[{Ind}];
Domain_Mag = Region[{Core, Air, Ind, Plate}];
Domain_NL = Region[{Core, Plate}];
}
Function {
mu0 = 4.e-7 * Pi;
murSteel = 100.;
nu [ Region[{Air, Ind}] ] = 1. / mu0;
//linear Case
// nu [ Region[{Core}]] = 1. / (murSteel*mu0);
// nu [ Region[{Plate}]] = 1. / (murSteel*mu0);
//Nonlinear Case
// nu [ Domain_NL ] = MatEIcore_interp_nu[$1] ;
// dhdb_NL [ Domain_NL ] = MatEIcore_dhdb_NL[$1];
nu [ Domain_NL ] = Steel_1020_interp_nu[$1] ;
dhdb_NL [ Domain_NL ] = Steel_1020_dhdb_NL[$1];
Sc = th_coil*ly_coil; //cross section of coil
Ns = 100; //Number of turns of the inductor Ind
Val_Current = 1.0; //current in the inductor
//calculation of current vectors
xc = 0.; // center of the coil
zc = 0.; // center of the coil
js[ Ind ] = Ns * Val_Current * Vector[-Sin[Atan2[Z[]-zc,X[]-xc]],
0.,
Cos[Atan2[Z[]-zc,X[]-xc]]] ;
}
Constraint {
{ Name MagneticVectorPotential_3D; Type Assign;
Case {
{ Region SurfaceGInf ; Value 0.; }
{ Region SurfaceGe0 ; Value 0.; }
}
}
{ Name Gauge; Type Assign;
Case {
{ Region Domain_Mag; SubRegion #{SurfaceGInf, SurfaceGe0,
SkinPlate, SkinCore}; Value 0.; }
}
}
}
Jacobian {
{ Name Vol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name CurlCurl ;
Case { {Type Gauss ;
Case {
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
}
}
}
}
}
Include "Magsta_a_3D_nl.pro"
PostOperation {
{ Name Map_a; NameOfPostProcessing MagSta_a_3D;
Operation {
Print[ a, OnElementsOf Domain_Mag, File "Magsta_a_3D_a.pos" ];
Print[ a_norm, OnElementsOf Domain_Mag, File "Magsta_a_3D_a_norm.pos" ];
Print[ b, OnElementsOf Domain_Mag, File "Magsta_a_3D_b.pos" ];
Print[ js, OnElementsOf DomainS_Mag, File "Magsta_a_3D_js.pos" ];
Print[ js, OnElementsOf DomainS_Mag, Format Table,File "js.txt" ];
}
}
}
// EOF
-------------- next part --------------
// Fichier "OneTooth-dat.pro"
mm = 1.e-3; // mm -> m
/*
x -> longueur -- gauche, droite
y -> hauteur -- haut, bas
z -> ?paisseur -- avant, arri?re
*/
// core
lx_core = 30. * mm;
ly_core = 100. * mm;
lz_core = 30. * mm;
// coil
th_coil = 10. * mm;
ri_coil = (lx_core^2. + lz_core^2.)^0.5 / 2. + 1. * mm;
re_coil = ri_coil + th_coil;
ly_coil = 50. * mm;
// plate
airgap = 2. * mm;
th_plate = 20. * mm;
ri_plate = ((lx_core/2.)^2. + (ly_core/2.)^2.)^0.5 + airgap;
re_plate = ri_plate+th_plate;
lz_plate = lz_core + 10. * mm;
// la bo?te ? air
bordure = 20. * mm;
axl = -re_plate - bordure; //left of airbox
axr = re_plate + bordure; //right of airbox
ayb = -re_plate - bordure; //bottom of airbox
ayt = re_plate + bordure; //top of airbox
azb = 0.; //back of airbox
azf = re_coil + bordure; //front of airbox
// les entit?s physiques
AIR = 9001;
CORE = 9002;
COIL = 9003;
PLATE = 9004;
AIR_INF = 9011;
SKIN_PLATE = 9012;
SKIN_CORE = 9013;
CUT_SYMMETRY = 9014;
// EOF