Download mag_nlin.procedur

Back to the list

   1 : * MAG_NLIN  PROCEDUR  CHAT      12/08/07    21:15:03     7481           
   2 :  'DEBPROC' MAG_NLIN   ETAB*'TABLE  '    ;
   3 : *********************************************************************   
   4 : *      estimation de mu aux points de gauss                             
   5 : *      multi-materiaux ferro_magnetiques                                        
   6 : *      materiaux  isotropes ou orthotropes
   7 : *********************************************************************   
   8 : *----------------------------------------------------------------------*
   9 : *                                                                      *
  10 : *
  11 : *                     PB  MAGNETOSTATIQUE  2D ou 3D                             
  12 : *                             ---------------                          *
  13 : *                                                                      *
  14 : *     RESOLUTION D'UN PROBLEME DE MAGNETOSTATIQUE  NON-LINEAIRE        *
  15 : *     EN REGIME PERMANENT,A L'AIDE DE LA METHODE DU POINT FIXE         *
  16 : *     ETAB, TABLE CONTENANT EN ENTREE :                                *
  17 : * OBLIGATOIRE                                                          *
  18 : *     INDICE 'SOUSTYPE'     THERMIQUE                                  *
  19 : *     INDICE 'AXI '   LOGIQUE VRAI   EN 2D  SI AXISYM ( PLAN DEFAUT)   *
  20 : *                                                                      *
  21 : *     INDICE 'TABNUSEC'      TABLE DES MATERIAUX FER INDICEE /MODELE   *
  22 : * OPTIONNEL                                                            *
  23 : *     INDICE 'CRITERE'      CRITERE DE CONVERGENCE  (10E-5 PAR DEFAUT) *        
  24 : *     INDICE 'OME'          POUR 2D UNIQUEMENT                         *
  25 : *                           COEFF AMORTISSEMENT OSCI  0< OME < 1.      *
  26 : *                                                                      *
  27 : *     INDICE 'NITER'        REACTUALISATION DE LA CONDUCTIVITE TOUTES  *
  28 : *                           LES NITER ITERATIONS (NITER=1 PAR DEFAUT)  *
  29 : *     INDICE 'NIVEAU'       NIVEAU DE MESSAGES (NIVEAU=0 PAR DEFAUT)   *
  30 : *     INDICE 'ITERMAX'      NOMBRE D'ITERATIONS MAXIMUM                *
  31 : *                           (ITERMAX=10 PAR DEFAUT)                    *
  32 : ************************************************************************
  33 : *  arguments fabriques dans les passages soit ds pot_vect ou pot_scal  *
  34 : *     INDICE 'FLUX'         FLUX EQUIVALENTS                           *
  35 : *     INDICE 'BLOQUE'    MATRICE DE BLOCAGE (CREEE PAR "BLOQUE")    *   
  36 : *     INDICE 'IMPOSE'       VALEURS IMPOSEES (CREE PAR "DEPI")         *
  37 : *     INDICE 'RIGCON '      RAIDEUR CONSTANTE                          *
  38 : *     INDICE 'RIGFER '      RAIDEUR VARIABLE                           *
  39 : *     ETAB CONTENANT EN SORTIE :                                       *
  40 : *                                                                      *
  41 : *     INDICE 'POTENTIEL'    POTENTIEL  RESULTAT                        *
  42 : *----------------------------------------------------------------------*
  43 :  CONVERGE = FAUX ;ETAB.'CONV'= FAUX ;                                   
  44 : 'REPETER' PROC 1 ;                                                      
  45 : 'SI' ( 'NEG' ( ETAB.'SOUSTYPE' ) 'THERMIQUE' ) ;                        
  46 :    'MESS' 'SOUS TYPAGE INCORRECT DE LA TABLE EN ENTREE|' ;              
  47 :    'QUITTER' PROC ;                                                     
  48 : 'FINSI' ;                                                               
  49 : 'SI' ( 'EXISTE' ETAB 'NIVEAU' ) ;                                       
  50 :    NIV_MESS = ETAB.'NIVEAU' ;                                           
  51 : 'SINON' ;                                                               
  52 :    NIV_MESS = 0 ;                                                       
  53 : 'FINSI' ;                                                               
  54 : 'SI' ( NIV_MESS '>EG' 1 ) ;                                             
  55 :    'SAUTER' 1 'LIGNE' ;                                                 
  56 : *                                                                       
  57 : *---  RECUPERATION DE L'INFORMATION CONTENUE DANS "ETAB"                
  58 : *                                                                       
  59 :  IVAL = 'VALE' 'DIME' ;
  60 : *  pour le 3D les blocages sont deja dans la matrice RIGCON
  61 :  'SI' (('EXISTE' ETAB 'BLOQUE' ) ET ('EGA' IVAL 2 )) ;
  62 :  MAT_BLO  = ETAB.'BLOQUE';
  63 :  'FINSI';
  64 :  RIG_CON  = ETAB.'RIGCON';                                              
  65 :  'SI' ( 'EXISTE' ETAB 'IMPOSE' );                                       
  66 :  VAL_IMPO = ETAB.'IMPOSE'   ;                                           
  67 :  'FINSI' ;                                                              
  68 :    'MESS' '***  DEBUT DE LA PROCEDURE "MAG_NONLIN"  ***' ;              
  69 : 'FINSI' ;                                                               
  70 : ***************                                                         
  71 : *---   Lecture de MU=MU(B)   mu secant   suprimee                       
  72 : ***************                                                         
  73 :  'SI' ( 'EXISTE' ETAB 'TABNUSEC' ) ;                                    
  74 :     TA_NUSEC = ETAB.'TABNUSEC' ;                                        
  75 :     'SINON' ;                                                           
  76 :     'MESS' '*  PRECISER la TABLE des NUsecants = 1/MU(B)  *';           
  77 :  'FINSI' ;                                                              
  78 :      AXI = FAUX ;                                                       
  79 :      'SI' ( 'EXISTE' ETAB 'AXI') ;     
  80 :       AXI = ETAB.'AXI' ;         
  81 :      'FINSI' ;
  82 :   si  axi ; mess 'calcul axisymetrique ' ;finsi ;    
  83 :      
  84 :      
  85 :                                                                
  86 : *MAIL_CHP= FER 'ET' ( 'EXTRA' RIG_CON 'MAIL' );                         
  87 :   MAIL_CHP = EXTR (ETAB.'POTENTIEL') MAILLAGE ;                         
  88 :   NBRE_NOE = 'NBNO' MAIL_CHP ;                                          
  89 :  VEC1= 'MANU' 'CHPO' MAIL_CHP 1 'T' ('PROG' NBRE_NOE * 1. ) ;           
  90 :  QTE_FLUX = ETAB.'RHS'  ;                                               
  91 :  'MENAGE' ;                                                             
  92 : 'SI' ( 'EXISTE' ETAB 'CRITERE' ) ;                                      
  93 :    EPSILON  = ETAB.'CRITERE'  ;                                         
  94 : 'SINON' ;                                                               
  95 :    EPSILON = 1.E-5 ;                                                    
  96 : 'FINSI' ;                                                               
  97 : 'SI' ( 'EXISTE' ETAB 'NITER' ) ;                                        
  98 :    NBRE_ITE = ETAB.'NITER' ;                                            
  99 :    'SI' ( NIV_MESS '>EG' 1 ) ;                                          
 100 :       'SAUTER' 1 'LIGNE' ;                                              
 101 :       'MESS' 'REACTUALISATION DE LA MATRICE DE CONDUCTIVITE A L ENTREE '
 102 :              'PUIS TOUTES LES ' NBRE_ITE 'ITERATIONS' ;                 
 103 :    'FINSI' ;                                                            
 104 : 'SINON' ;                                                               
 105 :    NBRE_ITE = 1 ;                                                       
 106 : 'FINSI' ;                                                               
 107 : 'SI' ( 'EXISTE' ETAB 'ITERMAX' ) ;                                      
 108 :    ITER_MAX = ETAB.'ITERMAX' ;                                          
 109 : 'SINON' ;                                                               
 110 :    ITER_MAX = 10 ;                                                      
 111 : 'FINSI' ;                                                               
 112 : **************************                                              
 113 : 'SI' ( 'EXISTE' ETAB 'IMPOSE');                                         
 114 :  FF1 = QTE_FLUX 'ET' VAL_IMPO ;                                         
 115 : 'SINON';                                                                
 116 :  FF1 = QTE_FLUX;                                                        
 117 : 'FINSI';                                                                
 118 :  'SI' ('EXISTE' ETAB 'OME') ; OME= ETAB.'OME' ; 'FINSI' ;
 119 :   'SI'  ('NON' ( 'EXISTE' ETAB 'POTENTIEL')) ;                          
 120 :    'SINON' ;                                                            
 121 :  'MESS'  '  ******************************************************';    
 122 :  'MESS'  '  ******************       REPRISE     *****************';    
 123 :  'MESS'  '  ******************************************************';    
 124 :        U1_T = 'EXCO' ( ETAB.'POTENTIEL')  'T' 'NOID' 'T' ;              
 125 :    'MESS' ' menage ' ; 'MENAGE' ;                                       
 126 :   'FINSI' ;                                                             
 127 :         'SI' ( NIV_MESS '>EG' 2 ) ;                                     
 128 :         'SAUTER' 1 'LIGNE' ;                                            
 129 :         'MESS' 'CHAMP THERMIQUE AVANT ITERATION ' ;                     
 130 :         'LISTE' U1_T ;                                                  
 131 :         'FINSI' ;                                                       
 132 :       DAN= 1.;                                                          
 133 : *                                                                       
 134 :        MOESP='REA' ;                                                    
 135 : *                                                                       
 136 :  'MESS' '            AMAX         AMIN          DU/U ';                 
 137 : *                                                                       
 138 :    MAXA= 'MAXI' U1_T ;MIXA= 'MINI' U1_T ; DAN = 0.;                     
 139 : ***********************************************************             
 140 : *---                          ... ITERATIONS ...                        
 141 : ***********************************************************             
 142 : NUM_ITE = 0 ;                                                           
 143 : IFOIS = 0 ;                                                             
 144 :  'MESS' IFOIS  MAXA MIXA    DAN ;                                       
 145 : 'REPETER' BOUC_1 ;                                                      
 146 :    NUM_ITE = NUM_ITE + 1 ;                                              
 147 :    IFOIS = IFOIS + 1 ;                                                  
 148 : *                                                                       
 149 :  INDE_TAB = 'INDE' TA_NUSEC ;                                           
 150 :  NPERM = 'DIME' INDE_TAB ;                                              
 151 :  i = 0 ;                                                                
 152 : *******************                                                     
 153 :  'REPETER' BOUPERM NPERM ;                                              
 154 :     i = i + 1 ;                                                         
 155 :     MOD_FERi = INDE_TAB.i ;                                             
 156 :     STN = TA_NUSEC.MOD_FERi ;
 157 :     ORTT = FAUX ;                                                       
 158 :          'SI' ('EXISTE' STN 'EV2' ) ;                                   
 159 :     ORTT = VRAI ;                                                       
 160 :     EV2  = STN.'EV2'  ;    DIR2 = STN.'DIR2' ;                          
 161 :     EV1 = STN.'EV1'   ;    DIR1 = STN.'DIR1' ;                          
 162 :          'SINON' ;                                                      
 163 :     EV1 = STN.'EV1'   ;                                                 
 164 :          'FINSI'      ;                                                 
 165 :     MAI_FERi = EXTRAIRE MOD_FERi 'MAIL' ;                               
 166 : *   calcul du champ dans le fer  ----->  modif de  mufer                
 167 : *      mufer (BMOD )   potentiel vecteur  on a donne MU(B)              
 168 : *      mufer (HMOD )   potentiel scalaire on a donne MU(H)              
 169 :       U1_FER='REDU' U1_T  MAI_FERi ;                                    
 170 :       DERIV= ( 'GRAD' MOD_FERi  U1_FER) ;                               
 171 :         DAX= 'EXCO'  'T,X' DERIV 'SCAL';                                
 172 :         DAY= 'EXCO'  'T,Y' DERIV 'SCAL';
 173 :         DAX= DAX 'CHAN' 'TYPE' 'SCALAIRE';
 174 :         DAY= DAY 'CHAN' 'TYPE' 'SCALAIRE';
 175 :                           
 176 :      'SI' AXI ;
 177 :         RFERi = (COOR 1 DERIV ) ** -1. ;
 178 :         DAX   = DAX  * RFERi;
 179 :         DAY   = DAY  * RFERi;
 180 :                                     
 181 :            'FINSI' ; 
 182 :                                                    
 183 :      'SI' ('EXISTE'  DERIV  'T,Z' ) ;                                   
 184 :     DAZ= 'EXCO'  'T,Z' DERIV 'SCAL'; 
 185 :     DAZ= DAZ 'CHAN' 'TYPE' 'SCALAIRE';                                   
 186 :     BB= (( DAY * DAY ) + ( DAX * DAX ) +( DAZ * DAZ )) ** .5 ;          
 187 :      'SINON' ;                                                          
 188 :         BB= ((( DAY * DAY ) + ( DAX * DAX ) ) ** .5 )        ;          
 189 :      'FINSI' ;                                                          
 190 :       BMAX= 'MAXIMUM'  BB ;BMIN= 'MINI' BB ;                                    
 191 :   'MESS'  IFOIS  MOESP ' B ou H max min  ' BMAX BMIN ;
 192 : *==============================================================                 
 193 : *  reactualisation  des permeabilites sur les sous modeles 
 194 : *==============================================================             
 195 :          'SI' ORTT  ;                                                   
 196 :    LIS_TE1  = EXTR EV1  ABSC ; LIS_TE2  = EXTR EV2  ABSC ;              
 197 :    LIS_CO1  = EXTR EV1  ORDO ; LIS_CO2  = EXTR EV2  ORDO ;              
 198 :     K1 = 'CHAN' 'NOEUD'  ('IPOL' BB  LIS_TE1 LIS_CO1) MOD_FERi ;                
 199 :     K2 = 'CHAN' 'NOEUD' ('IPOL' BB  LIS_TE2 LIS_CO2) MOD_FERi  ;                
 200 :                   'SI' ('EXISTE' DERIV 'T,Z' ) ;                  
 201 :  MAT_FERi= 'MATR' MOD_FERi  DIRECTION DIR1 DIR2 PARALLELE
 202 :      'K1'  K1      'K2'  K2    'K3'  K2 ;
 203 :                   'SINON' ;                                             
 204 : *  orthoptrope 2D
 205 :          'SI'   AXI        ;
 206 :  RFERL= 'CHAN' 'NOEUD'  MOD_FERi RFERI ; 
 207 :  MAT_FERi=  'MATR' MOD_FERi  DIRECTION DIR1
 208 :      'K1' ( K1  *  RFERL  )
 209 :      'K2' ( K2  *  RFERL  ) ;
 210 :          'SINON' ;
 211 :  MAT_FERi=  'MATR' MOD_FERi  DIRECTION DIR1
 212 :       'K1'  K1  'K2'  K2   ;     
 213 :          'FINSI' ;
 214 :                    'FINSI' ;                                            
 215 :          'SINON' ;                                                      
 216 : *  isotrope
 217 :       LIS_TE1  = EXTR EV1  ABSC ;                                       
 218 :       LIS_CO1  = EXTR EV1  ORDO ;
 219 :       K1 = 'CHAN' 'NOEUD' ('IPOL' BB  LIS_TE1 LIS_CO1) MOD_FERi ;               
 220 :     'SI' AXI ; 
 221 :   RFERL= 'CHAN' 'NOEUD' MOD_FERi RFERI ;                                        
 222 :   MAT_FERi =  'MATR'   MOD_FERi  'K' (  K1 * RFERL  ) ;
 223 :     'SINON' ;
 224 :   MAT_FERi= 'MATR'   MOD_FERi  'K'   K1 ;                  
 225 :     'FINSI' ;
 226 :             'FINSI' ;                                                   
 227 : ******** 'FINSI' ;
 228 :       'SI' (OME < 1. ) ;                                                
 229 :  KLAST = STN.'KLAST' ;                                                  
 230 :  MAT_FERi = COLI MAT_FERi OME   KLAST (1. - OME ) ;                     
 231 :             'FINSI' ;                                                   
 232 : *   on stoke la derniere carte de mu                                    
 233 :  STN.'KLAST' = MAT_FERi ;                                               
 234 : *    sommation sur les sous objets                                      
 235 :      'SI' ( 'EGA' i  1 ) ;                                              
 236 :  OBJM_FER = MOD_FERi ;                                                  
 237 :  CHAM_CND = MAT_FERi ;                                                  
 238 :     'SINON' ;                                                           
 239 :  CHAM_CND = CHAM_CND ET MAT_FERi ;                                      
 240 :  OBJM_FER =  OBJM_FER ET MOD_FERI ;                                     
 241 :     'FINSI' ;                                                           
 242 :      'FIN' BOUPERM ;                                                    
 243 : *******************                                                     
 244 :                                                                         
 245 :       CND2 = 'CONDUCTIVITE' OBJM_FER CHAM_CND ;                         
 246 :       RR2= CND2 'ET' RIG_CON  ;                                         
 247 :       RESID= ( FF1 - ( RR2 * U1_T )  )  'ENLEVER' 'FLX' ;               
 248 : *                                                                       
 249 : *                           tests  de convergence                       
 250 : *                                                                       
 251 :    RESID= RESID 'NOMC' 'SCAL';                                          
 252 :    NORES= 'XTX' RESID   ** .5 ;                                         
 253 :        ERRMAX= 'MAXI' ( 'ABS' RESID ) ;                                 
 254 :    MAXA= 'MAXI' U1_T ;MIXA= 'MINI' U1_T ;                               
 255 :    'DETR' RESID ;                                                       
 256 :  'MESS' IFOIS  MAXA MIXA    DAN ;                                       
 257 : *                                                                       
 258 : *                                                                       
 259 :    'SI'  (( NUM_ITE 'EGA' NBRE_ITE ) 'OU'                               
 260 :             ( IFOIS 'EGA' 1 ) );                                        
 261 : *****************************************************************       
 262 : * ---  REACTUALISATION DE LA MATRICE DE CONDUCTIVITE                    
 263 : *****************************************************************       
 264 :        MOESP='REA' ;                                                    
 265 : *      CI= ('LUMP' RR2 ) * VEC1 ; C2= ('LUMP' RR2 ( 'MOTS' T )) * VEC1; 
 266 : *      CI = (C2 - CI ) 'NOMC' 'SCAL';                                   
 267 :       'SI' (( 'EXISTE' ETAB 'BLOQUE') ET ( EGA IVAL 2 ));
 268 :      RIG1 = RR2  'ET' MAT_BLO ;
 269 :      'SINON';
 270 :       RIG1 = RR2 ;                                                      
 271 :      'FINSI' ;
 272 : *  mess ' avant resou  1 ' ;
 273 :          'SI' ( 'EXIST' ETAB 'GRAD' ) ;                                 
 274 :       U2 = 'RESOUDRE' 'GRAD' RIG1 FF1 ;                                 
 275 :           'SINON' ;                                                     
 276 :       U2 = 'RESOUDRE' RIG1 FF1 ;                                        
 277 :           'FINSI' ;                                                     
 278 : *  mess 'apres resou   1 ' ;
 279 :       NUM_ITE = 0 ;                                                     
 280 :    'SINON' ;                                                            
 281 : *****************************************************************       
 282 : * ---  RE-EQUILIBRAGE DU SECOND MEMBRE                                  
 283 : *****************************************************************       
 284 :       MOESP='   ' ;                                                     
 285 :       FF2 = ( RIG1  * U1_T ) - ( RR2 * U1_T ) ;                         
 286 :       FF3 = FF1 + FF2 ;                                                 
 287 : *mess ' avant resou reequilibrage ' ;
 288 :           'SI' ( EXISTE ETAB 'GRAD' ) ;                                 
 289 :       U2 = 'RESOUDRE' 'GRAD' RIG1 FF3 ;                                 
 290 :           'SINON' ;                                                     
 291 :       U2 = 'RESOUDRE' RIG1 FF3 ;                                        
 292 :           'FINSI' ;                                                     
 293 : *mess 'apres resou reequilibrage ';
 294 : *     'DETR' CND2 ;'DETR' RR1 ; 'DETR' RR2 ;                            
 295 :    'FINSI' ;                                                            
 296 : 
 297 : *****************************************************************       
 298 : *****************************************************************       
 299 :    U2_T = 'EXCO' U2 'T' 'NOID' 'T' ;                                    
 300 : *                                                                       
 301 :    'SI' ( NIV_MESS '>EG' 2 ) ;                                          
 302 :       'SAUTER' 1 'LIGNE' ;                                              
 303 :       'MESS' 'CHAMP THERMIQUE A L ITERATION :' IFOIS ;                  
 304 :       'LISTE' U2_T ;                                                    
 305 :    'FINSI' ;                                                            
 306 : *                                                                       
 307 :       CDIF= U2_T - U1_T ;                                               
 308 :       DAN= ('XTX'  CDIF)  / ( 'XTX' U1_T) ;                             
 309 :       DAN  = DAN  ** .5 ;                                               
 310 : *                                                                       
 311 : *        'SI' ( ERROR   < EPSILON ) ;                                   
 312 :         'SI' ( DAN    < EPSILON ) ;                                     
 313 :          CONVERGE = VRAI ;                                              
 314 :          'SINON' ;                                                      
 315 :          CONVERGE = FAUX ;                                              
 316 :          'FINSI' ;                                                      
 317 : *                                                                       
 318 : *---  LE CRITERE DE CONVERGENCE EST-IL SATISFAIT ?                      
 319 :      'MENAGE'  ;                                                        
 320 :                                                                         
 321 :    'SI'  CONVERGE  ;                                                    
 322 :       'SI' ( NIV_MESS '>EG' 1 ) ;                                       
 323 : *         'SAUTER' 1 'LIGNE' ;                                          
 324 :          'MESS' 'CONVERGENCE A L ITERATION    :' IFOIS    ;             
 325 :          'MESS' 'CRITERE DE CONVERGENCE       :' EPSILON  ;             
 326 :       'FINSI' ;                                                         
 327 :       'QUITTER' BOUC_1 ;                                                
 328 :    'FINSI' ;                                                            
 329 :    U1_T = U2_T ;                                                        
 330 :     'SI' ( 'EGA' IFOIS  ITER_MAX ) ;                                    
 331 : *      'SAUTER' 1 'LIGNE' ;                                             
 332 :       'MESS' 'PAS DE CONVERGENCE A L ITERATION :' ITER_MAX ;            
 333 :        ETAB.'CONV' = CONVERGE  ;                                        
 334 :       'QUITTER' BOUC_1 ;                                                
 335 :    'FINSI' ;                                                            
 336 : *              ON FAIT LE MENAGE                                        
 337 :  'FIN' BOUC_1 ;                                                         
 338 :         ETAB.'NBITER'= IFOIS;                                           
 339 : *                                                                       
 340 : *---  ARCHIVAGE DES RESULTATS DANS "ETAB"                               
 341 : *                                                                       
 342 : *ETAB.'KLAST'= EXCO 'SCAL' KLAST 'K' ;                                  
 343 :  ETAB.'POTENTIEL' = U2_T ;                                              
 344 : * ETAB.'CI'= CI ;                                                       
 345 : 'SI' ( NIV_MESS '>EG' 1 ) ;                                             
 346 :    'SAUTER' 1 'LIGNE' ;                                                 
 347 :    'MESS' '***  FIN DE LA PROCEDURE "MAG_NLIN"  ***' ;                  
 348 : 'FINSI' ;                                                               
 349 : 'FIN' PROC ;                                                            
 350 : 'FINPROC' ETAB ;                                                        
 351 :  
 352 :  
 353 :  
 354 :  
 355 :  

© Cast3M 2003 - All rights reserved.
Disclaimer