Télécharger biot_iseult.dgibi

Retour à la liste

optio echo 1;
************************************************************************
*        Test de quelques operatéurs pour le magnétisme (BIOT)         *
*                    et la mécanique des structures                    *
*                                                                      *
* Calcul du champ d'induction magnétique crée par un aimant d'IRM dans *
* l'espace autour de l'aimant, puis des forces de Laplace et des       *
* contraintes mécaniques dans l'aimant                                 *
*                                                                      *
* Application à l'aimant de l'IRM Iseult                               *
*                                                                      *
* Référence :                                                          *
* F. Nunio, G. Aubert, C. Berriaud, T. Schild and P. Vedrine,          *
* "Full Resolution Structural Finite Element Analysis of the CEA/      *
*  Iseult 11.75T MRI Coil"                                             *
* Publié dans IEEE Transactions on Applied Superconductivity,          *
* vol. 18, no. 2, pp. 916-919, Juin 2008                               *
* doi: 10.1109/TASC.2008.922544                                        *
*                                                                      *
************************************************************************
 
 
 
** Désactivation de l'echo des commandes
*OPTI 'ECHO' 0 ;
 
** Indicateur de tracé
itrac    = FAUX ;
 
 
 
*------------------------- P A R A M È T R E S ------------------------*
 
** Dimensions de la bobine principale
rint1    = 0.498 ;
rext1    = 0.941 ;
zinf1    = -1.915 ;
zsup1    =  1.915 ;
 
** Dimensions des bobines de blindage
rint2    = 1.97 ;
rext2    = 2.14 ;
zinf2    = -1.072 ;
zsup2    = -0.564 ;
rint3    = 1.97 ;
rext3    = 2.14 ;
zinf3    = 0.564 ;
zsup3    = 1.072 ;
 
** Perméabilité magnétique du milieu
mu       = 12.566E-7 ;
 
** Paramètres géométriques des inducteurs
*  Les bobines sont modelisées en un ensemble d'inducteurs circulaires à
*  section rectangulaire. Les sections sont décrites par les listes :
*  la1/la2 : liste des rayon internes/externes
*  lb1/lb2 : liste des hauteurs
la1      = PROG rint1 rint2 rint2 ;
la2      = PROG rext1 rext2 rext2 ;
lb1      = PROG zinf1 zinf2 zinf3 ;
lb2      = PROG zsup1 zsup2 zsup3 ;
 
** Liste des densités de courant (A.m-2) pour chaque bobine
*  ici le courant est uniquement dans la direction circonferentielle theta
*  en 2d axi : le sens positif correspond à UR^UZ (vers le fond de l'écran)
*  en 3d     : le sens positif s'enroule autour de l'axe -1 * y selon la
*              règle de la main droite
jprin    = -25.13E6 ;
jecra    =  32.56E6 ;
lj       = PROG jprin jecra jecra ;
 
** Taille du domaine où calculer le champ magnétique
r1       =  0. ;
r2       =  5. ;
h1       = -5. ;
h2       =  5. ;
 
 
 
*------------------ 2 D   A X I S Y M É T R I Q U E -------------------*
 
** Options de calcul
OPTI 'DIME' 2 'MODE' 'AXIS' ;
 
** Calcul magnétique (champ d'induction magnétique)
* maillage des points cibles (où calculer le champ d'induction)
OPTI 'ELEM' 'QUA4' 'DENS' 0.1 ;
cib      = (DROI (r1 h1) (r2 h1)) TRAN (0. (h2 - h1)) ;
ccib     = CONT cib ;
* définition/maillage des inducteurs
OPTI 'ELEM' 'QUA8' 'DENS' 0.1 ;
ni       = DIME lj ;
ind      = VIDE 'MAILLAGE' ;
tind     = TABL ;
moind    = VIDE 'MMODEL' ;
chbind   = VIDE 'MCHAML' ;
chjind   = VIDE 'MCHAML' ;
chb      = VIDE 'CHPOINT' ;
REPE bi ni ;
* densité de courant dans l'inducteur &bi
  j1       = EXTR lj &bi ;
* coordonnées de l'inducteur &bi
  a1       = EXTR la1 &bi ;
  a2       = EXTR la2 &bi ;
  b1       = EXTR lb1 &bi ;
  b2       = EXTR lb2 &bi ;
  bmoy     = 0.5 * (b1 + b2) ;
  btot     = b2 - b1 ;
* maillage de l'inducteur &bi
  pi1      = a1 b1 ;
  pi2      = a2 b1 ;
  pi3      = a1 b2 ;
  pi4      = a2 b2 ;
  ind1     = REGL (DROI pi1 pi2) (DROI pi3 pi4) ;
  ind      = ind ET ind1 ;
  tind . &bi = ind1 ;
* modèle pour l'inducteur &bi
  moi1 = MODE ind1 'MECANIQUE' 'ELASTIQUE' 'ORTHOTROPE' ;
  moind = moind ET moi1 ;
* calcul du champ d'induction magnétique dans le maillage de l'inducteur &bi
* (utile pour le calcul mécanique qui suit)
  bind     = BIOT 'INDU' moi1 'STRESSES' 'CERC' bmoy a1 a2 btot j1 mu ;
* calcul du champ d'induction magnétique aux points cibles
* (contribution de l'inducteur &bi)
  bcib     = BIOT 'INDU' cib 'CERC' bmoy a1 a2 btot j1 mu ;
* champ de densité de courant dans l'inducteur &bi
  chji     = MANU 'CHML' moi1 'JT' j1 'STRESSES' ;
* assemblage des champs
  chb      = chb    + bcib ;
  chbind   = chbind + bind ;
  chjind   = chjind + chji ;
FIN bi ;
cind     = CONT ind ;
* lignes de champ magnétique
chv      = EXCO (MOTS 'BR' 'BZ') chb (MOTS 'UR' 'UZ') ;
fc       = FCOURANT cib chv ;
maxf     = MAXI fc ;
minf     = MINI fc ;
nlig     = 21 ;
lf       = PROG minf 'PAS' ((maxf - minf) / (nlig + 1)) maxf ;
lig toto = @ISOSURF cib lf fc ;
chbm     = CHAN 'CHAM' chb cib ;
bm       = (((EXCO 'BR' chbm 'SCAL') ** 2) + ((EXCO 'BZ' chbm 'SCAL') ** 2)) ** 0.5 ;
blig     = PROI bm lig ;
 
** Calcul mécanique (déplacements/contraintes dans les inducteurs)
ind1     = tind . 1 ;
ind2     = tind . 2 ;
ind3     = tind . 3 ;
* modèle élastique orthotrope (1,2,3) <--> (r,z,t)
maind    = MATE moind 'YG1 ' 60.781E9 'YG2 ' 45.762E9 'YG3 ' 112.08E9
                      'NU12' 0.23988  'NU13' 0.19574  'NU23' 0.14228
                      'G12 ' 43.857E9 'DIRECTION' (1. 0.) ;
* matrice de rigidité
ri       = RIGI moind maind ;
* blocages
pz       = (ind1 POIN 'DROI' (0. 0.)    (1. 0.)) ET
           (ind2 POIN 'DROI' (0. zsup2) (1. zsup2)) ET 
           (ind3 POIN 'DROI' (0. zinf3) (1. zinf3)) ;
blz      = BLOQ 'UZ' pz ;
* forces de Laplace (= j ^ B)
fl       = (-1. * (EXCO 'JT' chjind 'FR') * (EXCO 'BZ' chbind 'FR')) ET
           (      (EXCO 'JT' chjind 'FZ') * (EXCO 'BR' chbind 'FZ')) ;
feq      = CNEQ moind fl ;
* champ de déplacements des bobines
u        = RESO (ri ET blz) feq ;
* déformée et contraintes
def      = DEFO u ind 200. ;
sig      = SIGM moind maind u ;
 
** Vérification des résultats
MESS ;
MESS '----------------------------------------------' ;
MESS '|               Calcul 2D axis               |' ;
MESS '----------------------------------------------' ;
MESS ;
* test du champ d'induction calculé au centre de l'aimant
ptest1   = cib POIN 'PROC' (0. 0.) ;
MESS '   Comp.  | Valeur calculee | Valeur attendue' ;
MESS '   -------+-----------------+----------------' ;
brt      = EXTR chb 'BR' ptest1 ;
bzt      = EXTR chb 'BZ' ptest1 ;
MESS '     Br   |' brt '     | 0.' ;
MESS '     Bz   |' bzt '     | 11.7' ;
* test du champ de déplacement au rayon interne de l'aimant
ptest2   = ind POIN 'PROC' (rint1 0.) ;
urt      = EXTR u 'UR' ptest2 ;
uzt      = EXTR u 'UZ' ptest2 ;
MESS '     Ur   |' urt '     | 7.75E-4' ;
MESS '     Uz   |' uzt '     | 0.' ;
MESS ;
* erreur si trop d'écart aux valeurs attendues
* on doit trouver Bz = 11,7 T    au centre
*              et Ur = 7.75E-4 m au rayon interne
err1     = MAXI 'ABS' brt ((bzt - 11.7) / 11.7) ;
SI (err1 > 1.E-3) ;
  ERRE 'Erreur de calcul du champ d''induction B' ;
FINSI ;
err2     = MAXI 'ABS' ((urt - 7.75E-4) / 7.75E-4) uzt ;
SI (err2 > 1.E-2) ;
  ERRE 'Erreur de calcul du champ de deplacements' ;
FINSI ;
 
** Quelques tracés
SI itrac ;
  OPTI 'PALE' (PALE 'PLAS') ;
  vb       = @VECOUL chb (MOTS 'BR' 'BZ') 3.E-2 ;
  TRAC vb (ccib ET cind) 'TITR' 'Champ B dans le domaine [2D axis]' ;
  TRAC blig lig (ccib ET cind) 'TITR' 'Lignes de champ B [2D axis]' ;
  bin      = CHAN 'CHPO' moind chbind 'SUPP' ;
  vb       = @VECOUL bin (MOTS 'BR' 'BZ') 1.E-2 ;
  TRAC vb ind 'TITR' 'Champ B dans les inducteurs [2D axis]' ;
  OPTI 'PALE' (PALE 'JET') ;
  fln      = CHAN 'CHPO' moind fl 'SUPP' ;
  vfl      = @VECOUL fln 3.E-10 ;
  TRAC vfl cind 'TITR' 'Forces de Laplace dans les inducteurs [2D axis]' ;
  TRAC sig moind def 'TITR' 'Contraintes (sur geometrie deformee) [2D axis]' ;
FINSI ;
 
 
SAUT 2 'LIGNE' ;
 
 
 
*----------------------- 3 D   V O L U M I Q U E ----------------------*
 
** Options de calcul
OPTI 'DIME' 3 'MODE' 'TRID' ;
 
** Calcul magnétique (champ d'induction magnétique)
* maillage des points cibles (où calculer le champ d'induction)
OPTI 'ELEM' 'CUB8' 'DENS' 0.25 ;
cib      = ((DROI ((-1. * r2) h1 h1) (r2 h1 h1)) TRAN (0. (h2 - h1) 0.))
             VOLU 'TRAN' (0. 0. (h2 - h1)) ;
* définition/maillage des inducteurs
OPTI 'ELEM' 'CU20' 'DENS' 0.15;
ni       = DIME lj ;
ind      = VIDE 'MAILLAGE' ;
tind     = TABL ;
moind    = VIDE 'MMODEL' ;
chbind   = VIDE 'MCHAML' ;
chjind   = VIDE 'MCHAML' ;
chb      = VIDE 'CHPOINT' ;
REPE bi ni ;
* densité de courant dans l'inducteur &bi
  j1       = EXTR lj &bi ;
* coordonnées de l'inducteur &bi
  a1       = EXTR la1 &bi ;
  a2       = EXTR la2 &bi ;
  b1       = EXTR lb1 &bi ;
  b2       = EXTR lb2 &bi ;
  bmoy     = 0.5 * (b1 + b2) ;
  btot     = b2 - b1 ;
* maillage de l'inducteur &bi
  pi1      = a1 b1 0. ;
  pi2      = a2 b1 0. ;
  pi3      = a1 b2 0. ;
  pi4      = a2 b2 0. ;
  ind1     = REGL (DROI pi1 pi2) (DROI pi3 pi4) ;
  ind1     = ind1 VOLU 60 'ROTA' 360. (0. 0. 0.) (0. 1. 0.) ;
  ELIM ind1 1.E-10 ;
  ind      = ind ET ind1 ;
  tind . &bi = ind1 ;
* modèle pour l'inducteur &bi
  moi1     = MODE ind1 'MECANIQUE' 'ELASTIQUE' 'ORTHOTROPE' ;
  moind    = moind ET moi1 ;
* calcul du champ d'induction magnétique dans le maillage de l'inducteur &bi
* (utile pour le calcul mécanique qui suit)
  bind     = BIOT 'INDU' moi1 'STRESSES' 'CERC' (0. bmoy 0.) (1. bmoy 0.) (0. bmoy 1.) a1 a2 btot j1 mu ;
* calcul du champ d'induction magnétique aux points cibles
* (contribution de l'inducteur &bi)
  bcib     = BIOT 'INDU' cib 'CERC' (0. bmoy 0.) (1. bmoy 0.) (0. bmoy 1.) a1 a2 btot j1 mu ;
* champ de densité de courant dans l'inducteur &bi
  x1 y1 z1 = COOR bind ;
  th1      = ATG z1 x1 ;
  chji     = j1 * ((NOMC 'JX' (-1. * (SIN th1))) ET 
                   (NOMC 'JY'        (0. * th1)) ET
                   (NOMC 'JZ'        (COS th1))) ;
* assemblage des champs
  chb      = chb    + bcib ;
  chbind   = chbind + bind ;
  chjind   = chjind + chji ;
FIN bi ;
aind     = ARET ind ;
 
** Calcul mécanique (déplacements/contraintes dans les inducteurs)
ind1     = tind . 1 ;
ind2     = tind . 2 ;
ind3     = tind . 3 ;
* modèle élastique orthotrope (1,2,3) <--> (z,r,t)
maind    = MATE moind 'YG1 ' 45.762E9 'YG2 ' 60.781E9 'YG3 ' 112.08E9
                       'NU12' 0.23988  'NU13' 0.14228  'NU23' 0.19574
                       'G12 ' 43.857E9 'G13 ' 38.833E9 'G23 ' 36.493E9
                       'RADIAL' (0. 0. 0.) (0. 1. 0.) ;
* matrice de rigidité
ri       = RIGI moind maind ;
* blocages
px       = ind POIN 'PLAN' (0. 0. 0.) (0. 1. 0.) (0. 0. 1.) ;
blx      = BLOQ 'UX' px ;
py       = (ind1 POIN 'PLAN' (0. 0. 0.)    (1. 0. 0.)    (0. 0. 1.)) ET
           (ind2 POIN 'PLAN' (0. zsup2 0.) (1. zsup2 0.) (0. zsup2 1.)) ET 
           (ind3 POIN 'PLAN' (0. zinf3 0.) (1. zinf3 0.) (0. zinf3 1.)) ;
bly      = BLOQ 'UY' py ;
pz       = ind POIN 'PLAN' (0. 0. 0.) (1. 0. 0.) (0. 1. 0.) ;
blz      = BLOQ 'UZ' pz ;
* forces de Laplace (= j ^ B)
fl       = (((EXCO 'JY' chjind 'FX') * (EXCO 'BZ' chbind 'FX')) - ((EXCO 'JZ' chjind 'FX') * (EXCO 'BY' chbind 'FX'))) ET
           (((EXCO 'JZ' chjind 'FY') * (EXCO 'BX' chbind 'FY')) - ((EXCO 'JX' chjind 'FY') * (EXCO 'BZ' chbind 'FY'))) ET
           (((EXCO 'JX' chjind 'FZ') * (EXCO 'BY' chbind 'FZ')) - ((EXCO 'JY' chjind 'FZ') * (EXCO 'BX' chbind 'FZ'))) ;
feq      = CNEQ moind fl ;
* champ de déplacements des bobines
u        = RESO (ri ET blx ET bly ET blz) feq ;
* déformée et contraintes
def      = DEFO u ind 200. ;
sig      = SIGM moind maind u ;
 
** Vérification des résultats
MESS ;
MESS '----------------------------------------------' ;
MESS '|                  Calcul 3D                 |' ;
MESS '----------------------------------------------' ;
MESS ;
* test du champ d'induction calculé au centre de l'aimant
ptest1   = cib POIN 'PROC' (0. 0. 0.) ;
MESS '   Comp.  | Valeur calculee | Valeur attendue' ;
MESS '   -------+-----------------+----------------' ;
bxt      = EXTR chb 'BX' ptest1 ;
byt      = EXTR chb 'BY' ptest1 ;
bzt      = EXTR chb 'BZ' ptest1 ;
MESS '     Bx   |' bxt '     | 0.' ;
MESS '     By   |' byt '     | 11.7' ;
MESS '     Bz   |' bzt '     | 0.' ;
* test du champ de déplacement au rayon interne de l'aimant
ptest2   = ind POIN 'PROC' (rint1 0. 0.) ;
uxt      = EXTR u 'UX' ptest2 ;
uyt      = EXTR u 'UY' ptest2 ;
uzt      = EXTR u 'UZ' ptest2 ;
MESS '     Ux   |' uxt '     | 7.75E-4' ;
MESS '     Uy   |' uyt '     | 0.' ;
MESS '     Uz   |' uzt '     | 0.' ;
MESS ;
* erreur si trop d'écart aux valeurs attendues
* on doit trouver By = 11,7 T    au centre
*              et Ux = 7.75E-4 m au rayon interne
err1     = MAXI 'ABS' bxt ((byt - 11.7) / 11.7) bzt ;
SI (err1 > 1.E-3) ;
  ERRE 'Erreur de calcul du champ d''induction B' ;
FINSI ;
err2     = MAXI 'ABS' ((uxt - 7.75E-4) / 7.75E-4) uyt uzt ;
SI (err2 > 1.5E-2) ;
  ERRE 'Erreur de calcul du champ de deplacements' ;
FINSI ;
 
** Quelques tracés
SI itrac ;
  pt1      = (0. zinf1 0.) ET (0. zsup1 0.) ;
  OPTI 'PALE' (PALE 'PLAS') ;
  TRAC blig (lig ET pt1) 'TITR' 'Lignes de champ B [3D]' (PROG 1. 'PAS' 1. 12.) ;
  OPTI 'PALE' (PALE 'JET') ;
  TRAC (sig / 1.E6) (PROG 30. 'PAS' 10. 180.) 'COUP' (0. 0. 0.) (0. 0. 1.) (0. 1. 0.) 
  moind def 'TITR' 'Contraintes (sur geometrie deformee) [3D]' ;
FINSI ;
 
FIN ;
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales