Télécharger biot_iseult.dgibi

Retour à la liste

************************************************************************
*        Test de quelques operateurs pour le magnetisme (BIOT)         *
*                    et la mecanique des structures                    *
*                                                                      *
* Calcul du champ d'induction magnetique cree par un aimant d'IRM dans *
* l'espace autour de l'aimant, puis des forces de Laplace et des       *
* contraintes mecaniques dans l'aimant                                 *
*                                                                      *
* Application a l'aimant de l'IRM Iseult                               *
*                                                                      *
* Reference :                                                          *
* 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"                                             *
* Publie dans IEEE Transactions on Applied Superconductivity,          *
* vol. 18, no. 2, pp. 916-919, Juin 2008                               *
* doi: 10.1109/TASC.2008.922544                                        *
*                                                                      *
************************************************************************
 
** Indicateurs de trace/calcul
itrac = FAUX ;
i2dax = VRAI ;
i3d   = VRAI ;
 
 
 
*------------------------- P A R A M E T R E S ------------------------*
 
** Dimensions de la bobine principale
rint1 = 0.498 ;
rext1 = 0.941 ;
zinf1 = -1.915 ;
zsup1 =  1.915 ;
 
** Dimensions des bobines d'ecrantage
rint2 = 1.97 ;
rext2 = 2.14 ;
zinf2 = -1.072 ;
zsup2 = -0.564 ;
rint3 = 1.97 ;
rext3 = 2.14 ;
zinf3 = 0.564 ;
zsup3 = 1.072 ;
 
** Permeabilite magnetique du milieu
mu = 12.566E-7 ;
 
** Parametres geometriques des inducteurs
*  Les bobines sont modelisees en un ensemble d'inducteurs circulaires a
*  section rectangulaire. Les sections sont decrites 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 densites 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 a UR^UZ (vers le fond de l'ecran)
*  en 3d     : le sens positif s'enroule autour de l'axe -1 * y selon la regle de la main droite
jprin = -25.13E6 ;
jecra =  32.56E6 ;
lj  = PROG jprin jecra jecra ;
 
** Taille du domaine ou calculer le champ magnetique
r1 =  0. ; r2 = 5. ;
h1 = -5. ; h2 = 5. ;
 
 
*------------------ 2 D   A X I S Y M E T R I Q U E -------------------*
 
SI i2dax ;
 
** Options de calcul
  OPTI 'DIME' 2 'MODE' 'AXIS' ;
 
** On place les points cibles (ou calculer le champ d'induction magnetique)
  OPTI 'ELEM' 'QUA4' 'DENS' 0.1 ;
  cib    = (DROI (r1 h1) (r2 h1)) TRAN (0. (h2 - h1)) ;
  ccib   = CONT cib ;
 
** Calcul du champ d'induction magnetique provoque par chaque inducteur
  ni     = DIME lj ;
  ind    = VIDE 'MAILLAGE' ;
  tind   = TABL ;
  moind  = VIDE 'MMODEL' ;
  chbind = VIDE 'MCHAML' ;
  chjind = VIDE 'MCHAML' ;
  chb    = VIDE 'CHPOINT' ;
  REPE bi ni ;
* densite de courant dans l'inducteur &bi
    j1   = EXTR lj &bi ;
* coordonnees 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 = 0. btot ;
    de1  = VALE 'DENS' ;
    OPTI 'ELEM' 'QUA8' ;
    ind1 = DROI pi1 pi2 TRAN 'DINI' de1 'DFIN' de1 (0.5 * pi3) TRAN 'DINI' de1 'DFIN' de1 (0.5 * pi3) ;
    ind  = ind ET ind1 ;
    tind . &bi = ind1 ;
* modele pour l'inducteur &bi
    moi1 = MODE ind1 'MECANIQUE' 'ELASTIQUE' 'ORTHOTROPE' ;
    moind = moind ET moi1 ;
* champ de densite de courant dans l'inducteur &bi
    chji = MANU 'CHML' moi1 'JT' j1 'STRESSES' ;
* champ d'induction magnetique dans les elements de l'inducteur &bi (pour le calcul mecanique)
    bind = BIOT 'INDU' moi1 'STRESSES' 'CERC' bmoy a1 a2 btot j1 mu ;
* champ d'induction magnetique aux points cibles (contribution de l'inducteur &bi)
    bcib = BIOT 'INDU' cib 'CERC' bmoy a1 a2 btot j1 mu ;
* assemblage des champs
    chb    = chb    + bcib ;
    chbind = chbind + bind ;
    chjind = chjind + chji ;
  FIN bi ;
  cind = CONT ind ;
 
** Trace du champ d'induction magnetique
  SI itrac ;
    vb  = @VECOUL chb (MOTS 'BR' 'BZ') 3.E-2 ;
    TRAC vb (ccib ET cind) 'TITR' 'Champ B dans le domaine [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]' ;
  FINSI ;
 
** Calcul mecanique (deplacements/contraintes dans les inducteurs)
  ind1 = tind . 1 ;
  ind2 = tind . 2 ;
  ind3 = tind . 3 ;
* modele elastique 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 rigidite
  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')) ;
  SI itrac ;
    fln = CHAN 'CHPO' moind fl 'SUPP' ;
    vfl = @VECOUL fln 1.E-9 ;
    TRAC vfl cind 'TITR' 'Forces de Laplace' ;
  FINSI ;
  feq = CNEQ moind fl ;
* champ de deplacements des bobines
  u = RESO (ri ET blz) feq ;
* deformee et contraintes
  def = DEFO u ind 200. ;
  sig = SIGM moind maind u ;
  SI itrac ;
    TRAC sig moind def 'TITR' 'Contraintes' ;
  FINSI ;
 
** Tests
  OPTI 'ECHO' 0 ;
  MESS ;
  MESS '----------------------------------------------' ;
  MESS '|               Calcul 2D axis               |' ;
  MESS '----------------------------------------------' ;
  MESS ;
* test du champ d'induction calcule 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 deplacement au ryaon interne de l'aimant
  ptest2 = ind POIN 'PROC' (rint1 0.) ;
  urt    = EXTR u 'UR' ptest2 ;
  uzt    = EXTR u 'UZ' ptest2 ;
  MESS '     Ur   |' urt '     | 7.78E-4' ;
  MESS '     Uz   |' uzt '     | 0.' ;
  MESS ;
* erreur si trop d'ecart aux valeurs attendures
* on doit trouver Bz = 11,7 Tesla au centre
* et Ur = 7.8E-4 metre 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.8E-4) / 7.8E-4) uzt ;
  SI (err2 > 1.E-3) ;
    ERRE 'Erreur de calcul du champ de deplacements' ;
  FINSI ;
 
** Visualisation des lignes de champ
  MESS ;
  MESS '   Calcul des lignes de champ ...' ;
* calcul de la fonction de courant (potentiel)
  chv = EXCO (MOTS 'BR' 'BZ') chb (MOTS 'UR' 'UZ') ;
  fc  = FCOURANT cib chv ;
* maillage des lignes d'isovaleurs
  maxf = MAXI fc ;
  minf = MINI fc ;
  nlig = 21 ;
  lf   = PROG minf 'PAS' ((maxf - minf) / (nlig + 1)) maxf ;
  lig toto = @ISOSURF cib lf fc ;
* interpolation de l'amplitude du champ B sur les lignes de champ
  chbm = CHAN 'CHAM' chb cib ;
  bm = (((EXCO 'BR' chbm 'SCAL') ** 2) + ((EXCO 'BZ' chbm 'SCAL') ** 2)) ** 0.5 ;
  blig = PROI bm lig ;
  SI itrac ;
    TRAC fc cib 'TITR' 'Fonction de courant' ;
    OPTI 'ISOV' 'LIGN' ;
    TRAC fc cib (ccib ET cind) nlig 'TITR' 'Isovaleurs de la fonction de courant' ;
    OPTI 'ISOV' 'SURF' ;
    TRAC blig lig (ccib ET cind) 'TITR' 'Lignes de champ B' ;
  FINSI ;
 
  SAUT 2 'LIGNE' ;
 
FINSI ;
OPTI 'ECHO' 1 ;
 
 
 
*----------------------- 3 D   V O L U M I Q U E ----------------------*
 
SI i3d ;
 
** Options de calcul
  OPTI 'DIME' 3 ;
 
** On place les points cibles (ou calculer le champ d'induction magnetique)
  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)) ;
 
** Calcul du champ d'induction magnetique provoque par chaque inducteur
  ni     = DIME lj ;
  ind    = VIDE 'MAILLAGE' ;
  tind   = TABL ;
  moind  = VIDE 'MMODEL' ;
  chbind = VIDE 'MCHAML' ;
  chjind = VIDE 'MCHAML' ;
  chb    = VIDE 'CHPOINT' ;
  REPE bi ni ;
* densite de courant dans l'inducteur &bi
    j1 = EXTR lj &bi ;
* coordonnees 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 = 0. btot 0. ;
    de1 = VALE 'DENS' ;
    OPTI 'ELEM' 'CU20' ;
    ind1 = DROI pi1 pi2 TRAN 'DINI' de1 'DFIN' de1 (0.5 * pi3) TRAN 'DINI' de1 'DFIN' de1 (0.5 * pi3) ;
    ind1 = ind1 VOLU 60 'ROTA' 360. (0. 0. 0.) (0. 1. 0.) ;
    ELIM ind1 1.E-10 ;
    ind = ind ET ind1 ;
    tind . &bi = ind1 ;
* modele pour l'inducteur &bi
    moi1 = MODE ind1 'MECANIQUE' 'ELASTIQUE' 'ORTHOTROPE' ;
    moind = moind ET moi1 ;
* champ d'induction magnetique dans l'inducteur &bi (pour le calcul mecanique)
    bind = BIOT 'INDU' moi1 'STRESSES' 'CERC' (0. bmoy 0.) (1. bmoy 0.) (0. bmoy 1.) a1 a2 btot j1 mu ;
* champ d'induction magnetique aux points cibles mp1 (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 densite 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 ;
 
** Trace du champ d'induction magnetique
  SI itrac ;
    vb = @VECOUL chb (MOTS 'BX' 'BY' 'BZ') 7.E-2 ;
    TRAC vb aind 'TITR' 'Champ B dans le domaine [3D]' ;
    bin = CHAN 'CHPO' moind chbind 'SUPP' ;
    vb = @VECOUL bin (MOTS 'BX' 'BY' 'BZ') 7.E-2 ;
    TRAC vb aind 'TITR' 'Champ B dans les inducteurs [3D]' ;
  FINSI ;
 
** Calcul mecanique (deplacements/contraintes dans les inducteurs)
  ind1 = tind . 1 ;
  ind2 = tind . 2 ;
  ind3 = tind . 3 ;
* modele elastique 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 rigidite
  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'))) ;
  SI itrac ;
    fln = CHAN 'CHPO' moind fl 'SUPP' ;
    vfl = @VECOUL fln 1.E-9 ;
    TRAC vfl aind 'TITR' 'Forces de Laplace' ;
  FINSI ;
  feq = CNEQ moind fl ;
* champ de deplacements des bobines
  u = RESO (ri ET blx ET bly ET blz) feq ;
* deformee et contraintes
  def = DEFO u ind 200. ;
  sig = SIGM moind maind u ;
  SI itrac ;
    TRAC sig moind def 'TITR' 'Contraintes' ;
    TRAC sig moind def 'COUP' (0. 0. 0.) (1. 0. 0.) (0. 1. 0.) 'TITR' 'Contraintes' ;
  FINSI ;
 
** Tests
  OPTI 'ECHO' 0 ;
  MESS ;
  MESS '----------------------------------------------' ;
  MESS '|                  Calcul 3D                 |' ;
  MESS '----------------------------------------------' ;
  MESS ;
* test du champ d'induction calcule 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 deplacement au ryaon 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.8E-4' ;
  MESS '     Uy   |' uyt '     | 0.' ;
  MESS '     Uz   |' uzt '     | 0.' ;
  MESS ;
* erreur si trop d'ecart aux valeurs attendures
* on doit trouver By = 11,7 Tesla au centre
* et Ux = 7.8E-4 metre 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.8E-4) / 7.8E-4) uyt uzt ;
  SI (err2 > 1.5E-2) ;
    ERRE 'Erreur de calcul du champ de deplacements' ;
  FINSI ;
 
FINSI ;
OPTI 'ECHO' 1 ;
 
 
FIN ;
 
 
 
 
 
 

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