$$$$ TRANSNON * TRANSNON PROCEDUR PASCAL 14/01/23 21:15:03 7913 'DEBPROC' TRANSNON PRECED*'TABLE'; *----------------------------------------------------------------------* * * * T R A N S N O N * * --------------- * * * * RESOLUTION D'UN PROBLEME DE THERMIQUE TRANSITOIRE NON-LINEAIRE * * METHODE A UN PAS DE TEMPS ( THETA SCHEMA ) * * * * ETAB, TABLE CONTENANT EN ENTREE : * * * * INDICE 'TET1' CHAMP DE TEMPERATURE INITIAL AU PAS 0 * * INDICE 'RAYO' LOGIQUE VALANT VRAI POUR UNE CONDITION * * DE RAYONNEMENT * * INDICE 'EMISSIVITE' MCHAML DECRIVANT LES FACTEURS D'EMISSIVITE* * NOM DE LA COMPOSANTE : EMIS * * INDICE 'CELSIUS' LOGIQUE VALANT VRAI SI L'UNITE EST LE * * DEGRE CELSIUS (CAPITAL SI RAYONNEMENT) * * INDICE 'MOD_THE' OBJET MODELE THERMIQUE * * INDICE 'MOD_CON' OBJET MODELE CONVECTION * * INDICE 'BLOCAGES_THERMIQUES' MATRICE DE BLOCAGE * * INDICE 'MAT_THE' OBJET MATERIAU THERMIQUE. * * CE CHAMP PEUT AVOIR DES COMPOSANTES DE * * TYPE FLOTTANT OU EVOLUTION (ABS-ORD). * * ABS : VARIABLE QUELCONQUE * * ORD : VALEUR DE LA COMPOSANTE CONCERNE POUR* * LA VALEUR ABS. * * INDICE 'MAT_CON' OBJET MATERIAU CONVECTION * * INDICE 'CHARGEMENT' CHARGEMENT DECRIVANT LES * * VALEURS DES VARIABLES EXTERNES (EX: TE, * * FLUX,TEMPERATURES IMPOSEES ,...) * * INDICE 'PHASE' LOGIQUE VALANT VRAI SI CHGT DE PHASE * * INDICE 'TEMPS0' TEMPS INITIAL (CORRESPOND A TET1) * * INDICE 'TEM_CALC' LISTREEL : TEMPS DES RESULTATS A CALCULER * * INDICE 'RELAXATION_THETA' VALEUR DU COEFFICIENT DE RELAXATION * * (VALEUR PAR DEFAUT 0.5) * * INDICE 'SOUS_RELAXATION' VALEUR DU COEFF. DE SOUS-RELAXATION * * (VALEUR PAR DEFAUT 0.5) * * INDICE 'PRECISION' VALEUR DU CRITERE DE CONVERGENCE * * (VALEUR PAR DEFAUT 1.E-4) * * INDICE 'MAXITERATION' NOMBRE MAXIMUM D'ITERATIONS AUTORISEES * * INDICE PROJECTION LOGIQUE VALANT VRAI SI COUPLAGE ET SI LE * * MAILLAGE DE LA MECANIQUE ET DE LA THERMIQUE* * EST DIFFERENT * * * * ETAB CONTIENT EN SORTIE * * * * INDICE INITIAL(2) DERNIER CHAMP DE TEMPERATURE CALCULE * * * * REMARQUE : CERTAINES VARIABLES ONT DES NOMS IMPOSES. * * LE TEMPS : TEMPS * * LA TEMPERATURE : T * * LA TEMPERATURE EXTERIEURE : TE * * LES TEMPERATURES IMPOSEES : TIMP * * LES FLUX : Q * * * * LE NOM DES AUTRES VARIABLES EST LAISSE AU CHOIX DE * * L'UTILISATEUR * * * ******************************************************************** * Recuperation de l'information contenue dans ETAB * ******************************************************************** ERR_VAL = FAUX;STAB = TABLE ; ETAB=PRECED.'WTABLE'; * * LINCO1 = WTAB . 'INCO_TOT' ; LDUAL1 = WTAB . 'DUAL_TOT' ; LINCO2 = ('MOTS' 'TEMP') 'ET' LINCO1 ; NBDUA1 = 'DIME' LDUAL1 ; * * MAIL_TOT = 'EXTR' ETAB.'MOD_TOT' 'MAIL'; IRIB = FAUX ; 'SI' ('NEG' ETAB.'RIBLO_T' 'INCONNU'); IRIB = VRAI ; RIBL_VAL = ETAB.'RIBLO_T' ; LISE_VAL = ETAB.'LISEA_T' ; 'FINSI' ; * Parallelisation automatique des calculs ? IPARAL = 'EGA' ETAB.'PROCESSEURS' 'AUTOMATIQUE' ; * ------------- pour le changement de phase ----------- IPHASE=ETAB.'PHASE'; 'SI' IPHASE; BLOPHA=ETAB.'BLO_PHASE'; PRO_VAL=ETAB.'PROPPHAS' ; MATPHA= ETAB.'MAT_PHA' ;MODPHA= ETAB.MOD_PHA; 'FINSI'; * -------------- presence de l'advection conduction convection IADVEC=ETAB.'ADVECTION'; 'SI' IADVEC; MATADV= ETAB.'MAT_ADV';MODADV=ETAB.'MOD_ADV'; 'FINSI'; ICONDU= ETAB.'CONDUCTION'; 'SI' ICONDU; MATCOND=ETAB.'MAT_COND';MODCOND=ETAB.'MOD_COND';'FINSI'; ICONVE= ETAB.'CONVECTION' ; 'SI' ICONVE; MATCON=ETAB.'MAT_CON'; MODCON=ETAB.'MOD_CON'; 'FINSI'; IRADI = ETAB.'RAYO'; 'SI' IRADI ; MATRAY=ETAB.'MAT_RAY'; MODRTHM=ETAB.'MOD_RAY'; 'FINSI'; ITHM1=WTAB.'THM1'; 'SI' ITHM1;MATTHM1=ETAB.'MAT_THM'; MODTHM= ETAB.'MOD_THM';'FINSI'; *--------------- Traitement special pour le rayonnement ----------- 'SI'('EXISTE' ETAB 'CTE_STEFAN_BOLTZMANN') ; CTE_SB = ETAB.'CTE_STEFAN_BOLTZMANN' ; 'SINON' ; CTE_SB = 5.673E-8 ; ETAB.'CTE_STEFAN_BOLTZMANN' =CTE_SB; 'MESS' 'ATTENTION ! Constante de Stefan-Boltzmann en unités SI' ; 'MESS' ' soit ' CTE_SB ; 'FINSI' ; *------- Initialisation de la liste des temps à calculer --------- TE_CALC = ETAB.'TEM_CALC' ; DCALC = 'DIME' TE_CALC; *-------------- Calcul du premier pas de calcul ------------------- TEM0 = ETAB.'TEMPS0' ; TEM1 = 'EXTR' TE_CALC 1; T_PAS = TEM1 - TEM0; *-------------- modif pour récuperer le temps du dernier pas de calcul TEM_VAL = TEM1 ; *-------------- fin de modif ----------------------------------------- *---- Initialisation de la valeur du critère de fin d'itération --- ZPREC1 = ETAB . 'PRECISION' ; *--------- Initialisation du nombre d'iterations autorisees ------- 'SI' ('EXISTE' ETAB 'MAXITERATION'); MAX1 = ETAB.'MAXITERATION'; 'SINON'; MAX1 = 150; 'FINSI'; *-------- Recherche des variables selon lesquelles varie ---------- *---------------------- le champ de matériau ---------------------- LVAREX = 'EXTR' ETAB.'MAT_TOT' 'DEVA'; ILVA = 'DIME' LVAREX; 'SI' ('NEG' ILVA 0); NONLIN1 = VRAI; 'SI' (ETAB.'COUPLAGE'); CONTI=PRECED.'CONTINUATION'; *--- Creation du MCHAML contenant tous les parametres mecaniques -- *-------- Ces parametres sont determines au temps TEM0 ------------ *----- Ces parametres sont projetes sur le maillage thermique ----- DEP1 = 'CHAN' 'CHAM' CONTI.'DEPLACEMENTS' ETAB.'MO_TOT' 'NOEUD'; CON1 = 'CHAN' 'NOEUD' CONTI.'CONTRAINTES' ETAB.'MO_TOT' ; CHAR1_ME = DEP1 'ET' CON1; * 'SI' ('EXISTE' CONTI 'VITESSES'); * VIT1 = 'CHAN' 'CHAM' CONTI.'VITESSES' * ETAB.'MO_TOT' 'NOEUD'; * CHAR1_ME = CHAR1_ME 'ET' VIT1; * 'FINSI'; 'SI' ('EXISTE' CONTI 'VARIABLES_INTERNES'); VAR1 = 'CHAN' 'NOEUD' CONTI.'VARIABLES_INTERNES' ETAB.'MO_TOT' ; * DEF1 ='CHAN' 'NOEUD' (ETAB.'DEI_EST') * ETAB.'MO_TOT' ; CHAR1_ME = CHAR1_ME 'ET' VAR1 ; * CHAR1_ME = CHAR1_ME 'ET' VAR1 'ET' DEF1; 'FINSI'; 'SI' ETAB.'PROJECTION'; CHAR_ME = 'CHAN' 'CHAM' ('PROI' MAIL_TOT CHAR1_ME) (ETAB.'MOD_TOT') 'RIGIDITE'; 'SINON'; CHAR_ME = 'CHAN' 'CHAM' ('CHAN' 'CHPO' CHAR1_ME ETAB.'MO_TOT') ETAB.'MO_TOT' 'RIGIDITE'; 'FINSI'; 'FINSI'; 'SINON'; NONLIN1 = FAUX; 'FINSI'; *-----------recherche des composantes de materiau qui varient et initia * --------- des logiques de recalcul des matrices.------------------- IKDEP=FAUX;ICDEP=FAUX; IECDEP=FAUX;IKVITE=FAUX; IKDDEP=FAUX;ICDDEP=FAUX; LCOREX= 'EXTR' ETAB.'MAT_TOT' 'COVA'; ILCOVA= 'DIME' LCOREX; 'SI' ('NEG' ILCOVA 0); 'REPETER' BOUVAR ILCOVA; ILO='EXTR' LCOREX &BOUVAR; 'SI'('EGA' ILO 'VITE') ; IKVITE=VRAI; 'FINSI'; 'SI'('EGA' ILO 'K ') ; IKDEP=VRAI ; 'FINSI' ; 'SI'('EGA' ILO 'C ') ; ICDEP=VRAI ; 'FINSI' ; 'SI'('EGA' ILO 'H ') ; IECDEP=VRAI ; 'FINSI' ; 'SI'('EGA' ILO 'RHO ') ; ICDEP=VRAI ; 'FINSI'; 'SI'('EGA' ILO 'KD ') ; IKDDEP=VRAI ; 'FINSI' ; 'SI'('EGA' ILO 'CDIF') ; ICDDEP=VRAI ; 'FINSI' ; 'FIN' BOUVAR; 'FINS'; 'SI' ETAB.'THM1' ; IKDEP=VRAI; ICDEP=VRAI; 'FINS' ; *------- Cas particulier ou il n'y a pas de CAPACITE ------------------* *- Ceci correspond au cas ou, au moins, une des 2 composantes 'C ' -* *- et/ou 'RHO ' n'est pas definie dans le materiau (facultatives). -* LCOMPM = 'EXTRAIRE' ETAB.'MAT_TOT' 'COMP' ; ICAPATH = ('EXISTE' LCOMPM 'C ') 'ET' ('EXISTE' LCOMPM 'RHO ') 'ET' ETAB.'FOR_THER' ; ICAPADI = ('EXISTE' LCOMPM 'CDIF') 'ET' ETAB . 'FOR_DIFF' ; ICAPA = ICAPATH 'OU' ICAPADI ; 'SI' ('NON' ICAPATH) ; ICDEP = FAUX ; 'FINSI' ; 'SI' ('NON' ICAPADI) ; ICDDEP = FAUX ; 'FINSI' ; *---------- Initialisation de l'unité de température -------------- 'SI' ('EXISTE' ETAB 'CELSIUS'); ICELS = ETAB.'CELSIUS' 'ET' ETAB.'FOR_THER' ; 'SINON'; ICELS = FAUX; 'FINSI'; *---------- Initialisation de la temperature de reference---------- 'SI' ICELS; TREF = 273.15 ; 'SINON'; TREF = 0. ; 'FINSI'; ETAB.'TREF'=TREF; *----------- Initialisation du champ de temperatures ------------ 'SI' ETAB.'FOR_THER' ; U0 = 'REDU' ('COPIER' ETAB.'TET1' 'GEOM') MAIL_TOT ; 'SI' ICELS ; CHTREF = (0. * U0) + TREF ; U0 = U0 + CHTREF ; 'FINSI' ; ETAB.'THER_COURANT' = ETAB.'TET1' ; 'SI' ETAB . 'FOR_DIFF' ; U0 = U0 'ET' ETAB.'CO1' ; ETAB.'THER_COURANT' = ETAB.'THER_COURANT' 'ET' ETAB.'CO1' ; 'FINS' ; 'SINO' ; U0 = 'REDU' ('COPI' ETAB.'CO1' 'GEOM') MAIL_TOT ; ETAB.'THER_COURANT' = ETAB.'CO1' ; 'FINS' ; *---------- y-a-t-il definition des coefficients lambda ---------- 'SI' ( 'EXISTE' ETAB 'RELAXATION_THETA' ) ; LAMBDA1 = ETAB.'RELAXATION_THETA' ; 'SINON' ; LAMBDA1=1.; ETAB.'RELAXATION_THETA' =LAMBDA1; 'FINSI' ; * En regime stationnaire (sans CAPACITE), on force LAMBDA1 a 1. * Sinon on ne resout pas le bon probleme ! 'SI' ('NON' ICAPA) ; LAMBDA1 = 1. ; 'FINSI' ; 'SI' IADVEC; LAMBDA1 = 1. ; 'FINSI' ; I_LAMBD1 = 0 ; 'SI' (' Calcul du premier pas par la theta-methode ' ; 'FINSI'; 'MESS' ' *** Données utilisees du probleme :'; 'MESS' ' > Champ de temperatures a l instant : ' TEM0; 'SI' NONLIN1; 'MESS' ' > K,C ou H variables.'; 'SINON'; 'MESS' ' > K,C et H constants.'; 'FINSI'; 'SI' ICELS; 'MESS' ' > Unité de temperature utilisee : Celsius.'; 'SINON'; 'MESS' ' > Unité de temperature utilisee : Kelvin.'; 'FINSI'; 'SI' ('EGA' (ETAB.'PROCEDURE_THERMIQUE') 'NONLINEAIRE'); 'MESS' ' > Coefficient de relaxation : ' LAMBDA1; 'SINON'; 'MESS' ' > Coefficient de relaxation pour la theta-méthode : ' LAMBDA1 ; 'MESS' ' > Coefficient de relaxation pour la méthode Dupont2 : ' ETAB.'RELAXATION_DUPONT' ; 'MESS' ' > Coefficient de sous-relaxation : ' LAMBDA2; 'FINSI'; 'SI' IPHASE; 'MESS' ' > Changement de phase prevu.'; 'SINON'; 'MESS' ' > Pas de changement de phase prevu.'; 'FINSI'; 'SI' ('NEG' ETAB.'BLOCAGES_THERMIQUES' 'INCONNU'); 'MESS' ' > Presence de temperatures imposees.'; 'SINON'; 'MESS' ' > Pas de temperatures imposees.'; 'FINSI'; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'Q '); 'MESS' ' > Présence de flux.'; 'SINON'; 'MESS' ' > Pas de presence de flux.'; 'FINSI'; 'SI' ICONDU; 'MESS' ' > Presence de conduction'; 'SINON'; 'MESS' ' > Pas de presence de conduction '; 'FINSI'; 'SI' ICONVE; 'MESS' ' > Presence de convection.'; 'SINON'; 'MESS' ' > Pas de terme de convection.'; 'FINSI'; 'SI' IRADI; 'MESS' ' > Presence de rayonnement.'; 'SINON'; 'MESS' ' > Pas de presence de rayonnement.'; 'FINSI'; 'SI' IADVEC; 'MESS' ' > Presence d advection.'; 'SINON'; 'MESS' ' > Pas de presence d advection'; 'FINSI'; 'MESS' ' '; 'FINSI'; ******************************************************************** * Boucle sur les pas de temps * ******************************************************************** 'REPETER' BOUC_SO ; ITE3 = &BOUC_SO; * 'SI' ('EGA' I_LAMBD1 0) ; TEMM = (LAMBDA1 * TEM1) + ((1 - LAMBDA1) * TEM0); * 'SINON' ; * 'SI' ('EGA' I_LAMBD1 1) ; TEMM = TEM0 ; 'FINSI' ; * 'SI' ('EGA' I_LAMBD1 2) ; TEMM = 0.5 * (TEM0+TEM1) ; 'FINSI' ; * 'SI' ('EGA' I_LAMBD1 3) ; TEMM = TEM1 ; 'FINSI' ; * 'FINSI' ; U_BOU = U0 ; TH_COUR = U_BOU ; 'SI' ICELS ; TH_COUR = U_BOU - CHTREF ; 'FINSI' ; ETAB.'THER_COURANT' = TH_COUR ; 'MESS' ' '; 'MESS' '---- Calcul du champ de température au temps :'TEM1; 'MESS' ' '; *---- initialisation du materiau 'SI' NONLIN1; *-------------- Cas de la temperature et du temps --------------- CHVAR = 'CHAN' 'CHAM' TH_COUR ETAB.'MOD_TOT' 'RIGIDITE' ; CHVAR2 = 'MANU' 'CHML' ETAB.'MOD_TOT' 'TEMP' TEMM 'RIGIDITE' ; CHVAR = CHVAR 'ET' CHVAR2 ; *-------------- Cas des variables mecaniques -------------------- *---- si il y a couplage on les "rentre" toutes dans CHVAR ------ 'SI' (ETAB.'COUPLAGE'); CHVAR = CHVAR 'ET' CHAR_ME; 'FINSI'; *-------------- Mise a jour des variables externes -------------- 'SI' (ETAB.'PROCEDURE_PARA_THER') ; PARATHER PRECED TEMM ; 'FINSI' ; *-------------- Cas des autres variables -------------------- 'REPETER' BOU_NU2 ILVA; ITE2 = &BOU_NU2; MOTI = 'EXTR' LVAREX ITE2; 'SI' ('EXISTE' LINCO2 MOTI) ; 'ITER' BOU_NU2 ; 'FINS' ; 'SI' (ETAB.'COUPLAGE'); 'SI' ('EXISTE' CHAR_ME MOTI ); 'SI' ('NEG' LAMBDA1 0.); 'MESS' '*** Le materiau thermique depend de parametres mecaniques'; 'ERRE' '*** Il faut initialiser SOUS-RELAXATION à 0.'; 'FINSI'; 'ITERER' BOU_NU2; 'FINSI'; 'FINSI'; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') MOTI); CH_1 = 'TIRE' (ETAB.'CHARGEMENT') MOTI TEMM; TYP_1 = 'TYPE' CH_1; * on suppose que les chargements qui ne sont pas d'origine mecanique * sont appuyes sur le maillage thermique . 'SI' ('EGA' TYP_1 'CHPOINT '); CHVAR2 = 'CHAN' 'CHAM' CH_1 ETAB.'MOD_TOT' 'RIGIDITE'; 'SINON'; CHVAR2 = 'CHAN' 'RIGIDITE' CH_1 ETAB.'MOD_TOT'; 'FINSI'; CHVAR = CHVAR + CHVAR2; 'SINON'; 'MESS' '*** Il manque une donnee du chargement ***'; 'MESS' '*** Variable concernee : 'MOTI; 'ERREUR' ' '; 'FINSI'; 'FIN' BOU_NU2; 'FINSI'; *----------------- Calcul du champ de matériau a T= U_BOU--------------- 'SI' NONLIN1; CHVAR1 = 'REDU' CHVAR (ETAB.'MOD_TOT'); CHMAT = 'VARI' 'NUAG' (ETAB.'MOD_TOT') (ETAB.'MAT_TOT') CHVAR1 'RIGIDITE'; ETAB.'CHMAT' = CHMAT; 'SI' (ETAB . 'FOR_DIFF' 'OU' ETAB.'THE1') ; CHMAT1 = 'REDU' CHMAT ETAB.'MOD_THE' ; 'SI' ETAB.'THM1' ; ma ff0 = @MATETHM ETAB.'MOD_THM' TH_COUR ; CHMAT1 = CHMAT1 'ET' ma ; 'FINS' ; 'SINO' ; 'SI' ETAB.'THM1' ; CHMAT1 ff0 = @MATETHM ETAB.'MOD_THM' TH_COUR ; 'FINS' ; 'FINS' ; 'SINON'; CHMAT = ETAB.'MAT_TOT'; ETAB.'CHMAT' = CHMAT; 'SI' (ETAB . 'FOR_DIFF' 'OU' ETAB.'THE1') ; CHMAT1 = 'REDU' CHMAT ETAB.'MOD_THE' ; 'SI' ETAB.'THM1' ; ma ff0 = @MATETHM ETAB.'MOD_THM' TH_COUR ; CHMAT1 = CHMAT1 'ET' ma ; 'FINS' ; 'SINO' ; 'SI' ETAB.'THM1' ; CHMAT1 ff0 = @MATETHM ETAB.'MOD_THM' TH_COUR ; 'FINS' ; 'FINS' ; 'FINSI'; 'SI' ICONDU; MATCOND='REDU' CHMAT1 ETAB.'MOD_COND';'FINSI'; 'SI' ICONVE; MATCON='REDU' CHMAT1 ETAB.'MOD_CON' ;'FINSI'; 'SI' IRADI ; MATRAD='REDU' CHMAT1 ETAB.'MOD_RAY' ;'FINSI'; 'SI' IADVEC; MATADV='REDU' CHMAT1 ETAB.'MOD_ADV' ;'FINSI'; 'SI' ITHM1 ; MATTHM='REDU' CHMAT1 ETAB.'MOD_THM' ;'FINSI'; 'SI' IPHASE; MATPHA='REDU' CHMAT1 ETAB.'MOD_PHA' ;'FINSI'; *--- Partitionnement des modeles pour la parallelisation automatique --- 'OPTION' 'PARA' IPARAL ; 'SI' IPARAL ; 'SI' ICONDU; MODRCOND= 'PART' MODCOND ETAB.'NBPART' ; MATRCOND= 'REDU' CHMAT1 MODRCOND; 'FINSI'; 'SI' ICONVE; MODRCON = 'PART' MODCON ETAB.'NBPART' ; MATRCON = 'REDU' CHMAT1 MODRCON; 'FINSI'; 'SI' IADVEC; MODRADV = 'PART' MODADV ETAB.'NBPART' ; MATRADV = 'REDU'CHMAT1 MODRADV; 'FINSI' ; * 'SI' IPHASE; MODRPHA = 'PART' MODPHA ETAB.'NBPART' ; * MATRPHA = 'REDU'CHMAT1 MODRPHA; * 'FINSI'; 'SI' ETAB . 'FOR_DIFF' ; MODDRED = 'PART' ETAB.'MOD_DIF' ETAB.'NBPART' ; CHMDRED = 'REDU' CHMAT1 MODDRED ; 'FINSI' ; 'SI' ITHM1; MODRTHM= 'PART' MODTHM ETAB.'NBPART' ; MATRTHM='REDU' CHMAT1 MODRTHM ; 'FINSI' ; 'SINON' ; 'SI' ICONDU; MODRCOND=MODCOND; MATRCOND=MATCOND;'FINSI'; 'SI' ICONVE; MODRCON=MODCON; MATRCON=MATCON; 'FINSI'; 'SI' IADVEC; MODRADV=MODADV; MATRADV=MATADV; 'FINSI'; 'SI' ITHM1; MODRTHM=MODTHM; MATRTHM=MATTHM; 'FINSI'; * 'SI' IPHASE; MODRPHA=MODPHA; MATRPHA=MATPHA; 'FINSI'; 'SI' ETAB . 'FOR_DIFF' ; MODDRED = ETAB.'MOD_DIF' ; CHMDRED = CHMAT1 ; 'FINSI' ; 'FINSI' ; *--- Calcul des matrices d'iteration On les recalcule a cause de perso IKDETO=FAUX;IBCONS=FAUX; 'SI' ICONDU ; KDETO= 'COND' MODRCOND MATRCOND; IKDETO=VRAI; 'FINSI'; 'SI' ITHM1; KDCC= 'COND' MODRTHM MATRTHM; 'SI' IKDETO; KDETO=KDETO 'ET' KDCC; 'SINON'; KDETO=KDCC;IKDETO=VRAI; 'FINSI'; 'FINSI'; 'SI' (('NON' IKDEP) 'ET' IKDETO); BCONS= KDETO;IBCONS=VRAI; 'FINSI'; 'SI' ICONVE; KDCCG= 'COND' MODRCON MATRCON; 'SI' IKDETO; KDETO=KDETO 'ET' KDCCG; 'SINON'; KDETO=KDCCG;IKDETO=VRAI; 'FINSI'; 'SI' ('NON' IECDEP); 'SI' IBCONS ; BCONS=BCONS 'ET' KDCCG; 'SINON'; BCONS=KDCCG;IBCONS=VRAI; 'FINSI'; 'FINSI'; 'FINSI'; 'SI' IADVEC; * la conduction de du model d'advection est deja calculer grace a icondu 'SI' ('NON' IKVITE); KDCC2= 'ADVE' MODRADV MATRADV; 'SI' IKDETO; KDETO=KDETO 'ET' KDCC2; 'SINON'; KDETO= KDCC2 ;IKDETO=VRAI; 'FINSI'; 'SI' IBCONS; BCONS = BCONS 'ET' KDCC2; 'SINON'; BCONS= KDCC2;IBCONS=VRAI; 'FINSI'; 'FINSI'; 'FINSI'; 'SI' ETAB . 'FOR_DIFF' ; KDEDI = 'RIGI' MODDRED CHMDRED ; 'SI' IKDETO; KDETO=KDETO 'ET' KDEDI; 'SINON'; KDETO=KDEDI;IKDETO=VRAI; 'FINSI'; 'SI' ('NON' IKDDEP); 'SI' IBCONS; BCONS=BCONS 'ET' KDEDI; 'SINON'; BCONS= KDEDI;IBCONS=VRAI; 'FINSI'; 'FINSI'; 'FINSI' ; * AA = KDETO * LAMBDA1 ; CCONS = AA ; BB = KDETO ; 'SI' ('NON' IBCONS) ; BCONS = EXTR KDETO 'RIGI' 'MULT'; IBCONS=VRAI; 'FINSI'; * 'SI' ICAPA ; 'SI' ICAPATH ; 'SI' ICONDU; CDET= ( 'CAPA' MODRCOND MATRCOND)/T_PAS; AA= CDET 'ET' AA ; 'SI' ICDEP ; CCONS = CDET 'ET' CCONS ; 'FINSI' ; 'FINSI'; 'SI' ITHM1; CDET= ( 'CAPA' MODRTHM MATRTHM)/T_PAS; AA= CDET 'ET' AA ; 'SI' ICDEP;CCONS = CDET 'ET' CCONS ; 'FINSI' ; 'FINSI'; 'FINSI'; 'SI' ICAPADI ; CDET = ('CAPA' MODDRED CHMDRED) / T_PAS ; AA = CDET 'ET' AA ; 'SI' ICDDEP ; CCONS = CDET 'ET' CCONS ; 'FINSI' ; 'FINSI' ; 'FINSI' ; * * 'SI' IPHASE; T0SSLX = TH_COUR 'ENLEVER' 'LX ' ; DEPPHA ='PHAJ' ETAB.'MOD_PHA' ETAB.'MAT_PHA' T0SSLX BLOPHA; MODPHA=ETAB.'MOD_PHA'; MATPHA= 'REDU' CHMAT MODPHA; AAQP ='EXCO' MATPHA 'QLAT' 'SCAL'; AARH='EXCO' MATPHA 'RHO' 'SCAL'; AASO=AAQP * AARH; QLAPHA=SOURCE ETAB.'MOD_PHA' AASO 'ELEM'; QLAPHA= QLAPHA / T_PAS ; 'FINSI'; * 'OPTION' 'PARA' FAUX ; *----------------- y-a-t-il des blocages --------------------------- * on le laisse ici car perso peut l'avoir modifie de memepour les autres * matrices IKBLOC=FAUX; 'SI' ETAB.'FOR_THER' ; IKBLOC=VRAI; MAT_BLOC = ETAB.'BLOCAGES_THERMIQUES' ; 'SI' ('EGA' MAT_BLOC 'INCONNU') ; MAT_BLOC = 'EXTR' KDETO 'RIGI' 'MULT' ; 'FINSI' ; 'FINSI' ; 'SI' ETAB . 'FOR_DIFF' ; MAT_BLDI = ETAB.'BLOCAGES_DIFFUSIONS' ; 'SI' ('EGA' MAT_BLDI 'INCONNU') ; MAT_BLDI = 'EXTR' KDEDI 'RIGI' 'MULT' ; 'FINSI' ; 'SI' IKBLOC; MAT_BLOC=MAT_BLOC 'ET' MAT_BLDI; 'SINON'; MAT_BLOC = MAT_BLDI; 'FINSI'; 'FINSI' ; ch_z = U0 'ENLEVER' 'LX' ; SECON = MAT_BLOC * ch_z ; AA = AA 'ET' MAT_BLOC; MAT_CHPO= MAT_BLOC * ch_z ; *------------------- y-a-t-il des temperatures imposees-------- ICHUIMP=FAUX; 'SI' ETAB.'FOR_THER' ; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'TIMP') ; CH_UIMP = 'TIRE' (ETAB.'CHARGEMENT') 'TIMP' TEM1 ; 'SI' ICELS ; CH_UIMP = CH_UIMP + TREF ; 'FINSI' ; ICHUIMP=VRAI; 'FINSI'; 'FINSI'; 'SI' ETAB . 'FOR_DIFF' ; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'CIMP') ; CH_CI = 'TIRE' (ETAB.'CHARGEMENT') 'CIMP' TEM1 ; 'SI' ICHUIMP; CH_UIMP= CH_UIMP 'ET' CH_CI; 'SINON'; CH_UIMP=CH_CI; 'FINSI'; ICHUIMP=VRAI; 'FINSI'; 'FINSI'; 'SI' ('NON' ICHUIMP); CH_UIMP='DEPI' MAT_BLOC 0.; 'FINSI'; MAT_CHPO = CH_UIMP - MAT_CHPO ; SECON = CH_UIMP - SECON ; *---------- Si convection on cherche TE ------------ 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'TECO' ) ; CH_TEC = 'TIRE' (ETAB.'CHARGEMENT') 'TECO' TEMM ; 'SI' ICELS ; CH_TEC = CH_TEC + TREF ; 'FINSI' ; 'FINSI'; *----------calcul convection second membre -------------- 'SI' ICONVE; CHMAT2 = 'REDU' CHMAT ETAB.'MOD_CON'; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'TECO' ) ; CHAL_TEI = 'CONV' ETAB.'MOD_CON' CHMAT2 CH_TEC ; MAT_CHPO = MAT_CHPO + CHAL_TEI; 'FINSI'; 'SI' ('NON' IECDEP) ; 'SI' ('EXISTE' (ETAB.'CHARGEMENT') 'TECO' ) ; SECON = SECON + CHAL_TEI; 'FINSI' ; 'FINSI'; 'FINSI'; *--------- SI terme de flux ------------------------------- 'SI' ETAB.'THM1' ; ma FTGC = @MATETHM ETAB.'MOD_THM' TH_COUR ; MAT_CHPO = MAT_CHPO-FTGC; SECON = SECON -FTGC; 'FINSI'; 'REPE' Boucle NBDUA1 ; MOTCOMP = 'EXTR' LDUAL1 &Boucle ; 'SI' ('EXIS' (ETAB.'CHARGEMENT') MOTCOMP) ; CH_DUA = 'TIRE' (ETAB.'CHARGEMENT') MOTCOMP TEMM ; MAT_CHPO = MAT_CHPO + CH_DUA ; SECON = SECON + CH_DUA ; 'FINS' ; 'FIN' Boucle ; *------------- calcul des termes rayonnement- et des autres----------- 'SI' (etab.'PROCEDURE_CHAR_THER'); IAPPEL=3; TAAT=charther preced temm IAPPEL; 'SI'('EXIS' TAAT 'ADDI_MATRICE'); mat_z = lambda1 * taat.'ADDI_MATRICE' ; ccons = mat_z 'ET' ccons ; aa = mat_z 'ET' aa ; bb=taat.'ADDI_MATRICE' 'ET' bb; 'FINSI'; 'SI' ('EXIS' TAAT 'ADDI_SECOND' ); MAT_CHPO= MAT_CHPO et taat.'ADDI_SECOND'; 'FINSI'; 'SI' ('EXIS' TAAT 'RAYO_VALEUR'); STAB.'RAYONNEMENT'=TAAT.'RAYO_VALEUR'; 'FINSI'; 'FINSI'; *------------- si phase insertion des blocages et des jeux --------- 'SI' IPHASE; AA = AA 'ET' BLOPHA; SECON = SECON 'ET' DEPPHA; MAT_CHPO = MAT_CHPO 'ET' DEPPHA; 'FINSI'; *--------------- Calcul du premier résidu -------------------------- B3 = BB * U0 ; MAT_CHPO = MAT_CHPO - B3 ; ******************************************************************** * Boucle sur les iterations * ******************************************************************** * on peut mettre n'importe quoi c'est pour * ne pas faire de tests dans la boucle ZCLIM = MAT_BLOC; ZNACCE =2;ZITAC= 0 ;IT= 0 ;ITDEP=3;'SI' IPHASE; ITDEP=5;'FINSI'; ITACC = 0; ACFP1 = 'COPIER' MAT_CHPO * 0.; ACFP2 = ACFP1 ; MAT_CHPP=ACFP1 ; ACFP3 = ACFP1 ; ACFEP1 = ACFP1; ACFEP2 = ACFP1 ; CORREC=0.; DU0 = 0. * U0; U00 = 'ENLEVER' U0 'LX'; 'SI' IPHASE ; U00MTREF = U00 ; 'SI' ICELS ; U00MTREF = U00 - CHTREF ; 'FINSI' ; 'FINSI' ; U1=U00; 'REPETER' BOU_IT1; *IT est le compteur de ETIQ, ITACC doit etre =< 0 pour qu'on accelere IT= IT + 1 ; ITACC = ITACC - 1; ZITAC = ZITAC + 1 ; ITE_IT = &BOU_IT1; AUTORISA=VRAI; BB= BCONS;CC=CCONS; * * selection des fluxmax pour le changement de phase * 'SI' IPHASE; 'SI' (ITE_IT > 1) ; U1TRE = U1 'ENLEVER' 'LX' ; 'SI' ICELS ; U1TRE = U1TRE - CHTREF ; 'FINSI' ; U1TRE = U1TRE + (U1 'EXCO' 'LX' 'NOID' 'LX') ; 'SINON' ; U1TRE = U1 ; 'SI' ICELS ; U1TRE = U1 - CHTREF ; 'FINSI' ; AUTORISA = FAUX ; 'FINSI'; CHPPHA = 'EXCP' MODPHA MATPHA PRO_VAL QLAPHA U00MTREF U1TRE BLOPHA; 'FINSI'; * nombre max d'iteration autorisé est-il dépassé? 'SI' (ITE_IT '>' MAX1); 'MESS' '*** LE CALCUL N A PAS CONVERGE ***'; ETAB.'ERREUR' = VRAI; 'QUITTER' TRANSNON; 'FINSI'; * *--------------------------------------------------------------------- * La force motrice de l'iteration est fixee: RESIDU * on va calculer un nouveau champ de deplacement * * calcul de l'increment de l'increment de deplacement zdep1 * par resolution lineaire * * * acceleration de convergence * CORRECP = CORREC; CORREC = 0; ACFP0 = (MAT_CHPO - MAT_CHPP) 'ENLE' FLX ; ACFEP0 = ACFP0; ACFEP0 = ACFEP0 - CORRECP ; 'SI' ('MULT' IT ZNACCE) ; 'SI' (IT > ITDEP); CORREC = 'ACT3' ACFEP2 ACFEP1 ACFEP0 ACFP3 ACFP2 ACFP1 ACFP0 ; MAT_CHPO = MAT_CHPO - CORREC; 'FINSI'; 'FINSI'; SI (IT > 3); 'DETR' ACFP3; 'DETR' ACFEP2; 'FINSI'; ACFP3 = ACFP2 ; ACFP2 = ACFP1 ; ACFP1 = ACFP0 ; ACFEP2 = ACFEP1 ; ACFEP1 = ACFEP0 ; MENAGE; 'SI' (IRIB); 'SI' IPHASE; DU1 BID BID BID BID = 'RESOU' AA MAT_CHPO 'INIB' RIBL_VAL LISE_VAL CHPPHA; 'SINON'; DU1 BID BID BID = 'RESOU' AA MAT_CHPO 'INIB' RIBL_VAL LISE_VAL ; 'FINSI'; 'SINON'; 'SI' IPHASE; DU1 = 'RESOUDRE' AA MAT_CHPO CHPPHA; 'SINON'; DU1 = 'RESOUDRE' AA MAT_CHPO ; * mess ' ITeration ' it; * aaa= extr aa rigi mult ; mess ' aaa ' ; list aaa; * mess ' mat_chpo' ; list mat_chpo; * mess ' du1 ' ; list du1; 'FINSI'; 'FINSI'; ZRAID_T = 'EXTR' AA 'CONT'; 'SI' ('EXIS' ZRAID_T 'NITER'); RIBL_VAL = ZRAID_T.7; LISE_VAL = ZRAID_T.6; IRIB=VRAI; 'FINSI'; U1 = U00 + DU1; XXX1= ZCLIM * U1; MAT_CHPP= MAT_CHPO 'ENLEVER' 'FLX'; 'SI' ('EGA' I_LAMBD1 0) ; U_BOU1 = (LAMBDA1 * U1) + ((1 - LAMBDA1) * U0); 'SINON' ; 'SI' ('EGA' I_LAMBD1 1) ; U_BOU1 = U0 ; 'FINSI' ; 'SI' ('EGA' I_LAMBD1 2) ; U_BOU1 = 0.5 * (U0 + U1) ; 'FINSI' ; 'SI' ('EGA' I_LAMBD1 3) ; U_BOU1 = U1 ; 'FINSI' ; 'FINSI' ; TH_COUR = U_BOU1 ; 'SI' ICELS; TH_COUR = U_BOU1 - CHTREF; 'FINSI'; ETAB.'THER_COURANT' = TH_COUR ; *-------- Calcul du critère de fin d'itération ------------------- DUABS2 = 'ABS' (DU1 - DU0) ; ZPREC2 = 0. ; 'REPE' BINCO1 NBDUA1 ; LINCOI1 = 'MOTS' ('EXTR' LINCO1 &BINCO1) ; MDUABSI1 = 'MAXI' DUABS2 'AVEC' LINCOI1 ; NORMI1 = 'MAXI' 'ABS' U1 'AVEC' LINCOI1 ; 'SI' (NORMI1 '' ZPREC2) ; ZPREC2 = ZPRECI2 ; 'FINS' ; 'FIN' BINCO1 ; DU0 = DU1; 'MESS' ' Itération numéro : 'ITE_IT ' Valeur du critère : 'ZPREC2 ; 'SI' ((ZPREC2 '