* UNPAS PROCEDUR FD218221 21/01/15 21:15:01 10848 * *----------------------------------------------------------------------* * PROCEDURE UNPAS * * * * calcul d'un increment de solution en grand deplacement plastique * * par la methode des residus * *----------------------------------------------------------------------* * En entree * PRECED la table passee à PASAPAS * * kich nota champ de materiau : etablit au debut de procedure UNPAS * estimation des caracteristiques fin de pas (TI ): ZMAT22 -->> ZMATXX * caracteristiques debut de pas (TEMPS0) : ZMAT11 * puis dans la boucle de non-convergence * ZMAT1 initialise avec ZMAT11 * materiau fin de pas de temps : ZMAT2 sorti de COMP, * range dans WTAB.'MAT1' en sortie unpas ************************************************************************ * sortie STAB12 indice : * * * * DEPT increment de deplacement sur le pas * * SIGF contraintes a la fin du pas * * VARF variables internes a la fin du pas * * DFPF deformation inelastique a la fin du pas * * CONV logique valant vrai si pas de probleme de convergence * * DEFF deformations a la fin du pas si grandes deformations * *----------------------------------------------------------------------* * * Appel a PAS_MODL : * Mise a jour des indices de PRECED.WTABLE relatifs aux modeles * si PRECED.WTABLE.MODELE a ete modifie PAS_MODL PRECED ; *CB215821 : Recuperation de XPETIT (07/12/2016) * declanchement du recalcul de la matrice ITRCLC = 1 * WTAB.'DELTAITER'; 'SI' (ITRCLC < 20) ; ITRCLC = 20; 'FINSI'; WTAB = PRECED.'WTABLE' ; LAG_TOT=VRAI; conti= PRECED.'CONTINUATION'; estim= PRECED.'ESTIMATION' ; * * initialisation et reprise des valeurs des tables **************** * MXMYMZ = 'MOTS' 'MX' 'MY' 'MZ' 'MT' 'FP' 'FPQ' 'FTP' 'IMX' 'IMY' 'IMZ' 'IMT' 'IFP' 'IFPQ' 'IFTP' ; MXMFLX = 'MOTS' 'MX' 'MY' 'MZ' 'MT' 'FLX' 'FP' 'FPQ' 'FTP' 'IMX' 'IMY' 'IMZ' 'IMT' 'FLX' 'IFP' 'IFPQ' 'IFTP' ; MNPRIM = 'MOTS' 'UX' 'UY' 'UZ' 'UR' 'UT' 'RX' 'RY' 'RZ' 'RT' 'P' 'PQ' 'TP' 'ALFA' 'BETA' 'IUX' 'IUY' 'IUZ' 'IUR' 'IUT' 'IRX' 'IRY' 'IRZ' 'IRT' 'IP' 'IPQ' 'ITP' 'IALF' 'IBET' ; MNDUAL = 'MOTS' 'FX' 'FY' 'FZ' 'FR' 'FT' 'MX' 'MY' 'MZ' 'MT' 'FP' 'FPQ' 'FTP' 'FALF' 'FBET' 'IFX' 'IFY' 'IFZ' 'IFR' 'IFT' 'IMX' 'IMY' 'IMZ' 'IMT' 'IFP' 'IFPQ' 'IFTP' 'IFAL' 'IFBE' ; MLDEPL = 'MOTS' 'UX' 'UY' 'UZ' 'UR' 'UT' 'IUX' 'IUY' 'IUZ' 'IUR' 'IUT' 'ALFA' 'BETA' 'IALF' 'IBET' ; MLROTA = 'MOTS' 'RX' 'RY' 'RZ' 'RT' 'P' 'PQ' 'TP' 'IRX' 'IRY' 'IRZ' 'IRT' 'IP' 'IPQ' 'ITP'; MLDEFOR = 'MOTS' 'EPXX' 'EPYY' 'EPZZ' 'EPSS' 'EPTT' 'EPRR' 'GAXY' 'GAXZ' 'GAYZ' 'GAST' 'GASN' 'GATN' 'GARZ' 'GART' 'GAZT' 'GXY ' * 'CX ' 'CY ' 'CZ ' 'EPSE' ; MVPRIM = 'MOTS' 'VTX' 'VTY' 'VTZ' 'VTR' 'VTT' 'VWX' 'VWY' 'VWZ' 'VWT' 'VVP' 'VVPQ' 'VVTP' 'VALF' 'VBET' 'IVTX' 'IVTY' 'IVTZ' 'IVTR' 'IVTT' 'IVWX' 'IVWY' 'IVWZ' 'IVWT' 'IVVP' 'IVPQ' 'IVTP' 'IVAL' 'IVBE' 'VLX'; 'V1X ' 'V1Y ' 'V1Z ' 'V2X ' 'V2Y ' 'V2Z ' ; * definition de variables locales PASUNIL = FAUX; *'SI' ( 'EXIS' PRECED 'ECRIT' ) ; list wtab; 'FINS'; ICERAM = WTAB.'CERAMIQUE' ; IDYN = WTAB.'DYNAMIQUE'; IELANL = WTAB.'NON_LINEAIRE'; IENDOM = WTAB.'ENDOMMAGEMENT'; IFEFP = WTAB.'FEFP_FORMULATION' ; IFEFPUL= WTAB.'UPDATE_LAGRANGIAN'; IFTOL = 'NEG' WTAB.'FTOL' 'INCONNU' ; IGRD = WTAB.'GRANDS_DEPLACEMENTS'; IKSIA = WTAB.'K_SIGMA'; IKTAN = WTAB.'K_TANGENT' ; IMPLP = WTAB.'LIAISON_PERSISTANTE'; IMTOL = 'NEG' WTAB.'MTOL' 'INCONNU'; IPILOT = WTAB.'AUTOMATIQUE'; IPLAST = WTAB.'PLASTIQUE'; IPLAVI = WTAB.'IPLAVI'; IPREDIC = WTAB.'PREDICTEUR'; IRCON = WTAB.'RAIDCONST'; IRAUGLU = WTAB.'RAIDAUGM'; IRAUG = IRAUGLU ; AUTAUG = WTAB.'AUTOAUGM'; ISOL = WTAB.'CONSOLIDATION'; ISSTE = WTAB.'SUBSTEPPING'; ITHER = WTAB.'CHAR_THE' 'OU' WTAB.'FOR_THER'; IVIEXT = WTAB.'VISCO_EXTERNE'; IVIDOM = WTAB.'VISCODOMMAGE'; IVISCO = WTAB.'VISCOPLASTIQUE'; LOGDEF = WTAB.'CHAR_DEFI'; LOGPRE = WTAB.'CHAR_PRES' ; LOGPIL = WTAB.'CHAR_PILO' ; NITMA = WTAB.'NITERINTER_MAX'; NSSTE = WTAB.'NMAXSUBSTEPS'; POR1 = WTAB.'POR1' ; TI = WTAB.'T_FINAL'; EKREAC = WTAB.'REAC_GRANDS'; ekreac = 1e50; ZMAXIT = WTAB.'MAXITERATION' ; ZNACCE = VRAI; XCONV = 0.; DEPSTDM = 0.; ZNCONS = WTAB.'NITER_KTANGENT' ; ZPREC = WTAB.'PRECISION' ; ZPREK = WTAB.'PRECISINTER' ; ZPRECD = WTAB.'PRECDECHARGE' ; ZPRECM = WTAB.'PRECFLEX' ; ZCLIM0 = WTAB.'BLOCAGES_MECANIQUES'; *-- Autres initialisations en non-local helmhoktz -- 'SI' (LNLOC 'ET' ('EGA' WTAB.'NON_LOCAL' 'HELM')); TAHELM = WTAB.'HELMHOLTZ' ; NHELM = TAHELM . 'N_VARI_NL' ; 'FINSI' ; *-- On initialise ZCLIM qui va contenir l'ensemble des C.L. ZCLIM = ZCLIM0 ; DT = WTAB.'DT' ; TEMPS0=WTAB.'TEMPS0'; *--- doit on reactualiser la geometrie temp0 ne 0 et grand depl ? --- 'SI' WTAB.'RECALCUL' ; WTAB . 'RECARI' = VRAI ; WTAB . 'RECADET' = VRAI ; WTAB . 'REA_GEOM' = VRAI ; * on suppose que l'on est sur la bonne configuration * 'FORM' GEOM1; 'SINON'; GEOM1 = WTAB.'FOR0' ; WTAB . 'RECARI' = FAUX ; WTAB . 'RECADET' = FAUX ; 'FINS'; 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); WTAB.'RECARI'= VRAI; LAG_TOT=0; 'FINS'; 'SI' IRAUG; RIG_AUG = WTAB.'RIGIDITE_AUGMENTEE'; 'FINS'; 'SI' IRCON; RIG_CONS = WTAB.'RIGIDITE_CONSTANTE'; 'FINS'; *---------- chamelem etat estime pour la fin du pas de temps ------ ZMATFI = ZMAT22 ; *----------------------- nouveau chargement ---------------------- 'SI' ('NEG' TYP_2 'CHPOINT '); 'FINS'; 'SINON'; 'FINS'; 'SI' ('NEG' TYP_2 'CHPOINT '); 'FINS'; ZFEXT2 = ZFEXT2 '+' F2_FOR ; 'FINS'; 'SI' (LOGPRE 'ET' ('NON' IGRD)) ; MOP = WTAB.'MOD_PRE' ; 'SINON' ; 'FINS' ; ZFEXT2 = ZFEXT2 '+' ZFPEXT ; 'FINS' ; *--------------- Si il existe des deplacements imposes -------------- ZFEXT2 = ZFEXT2 + F2_DEP; 'FINS'; *--------------- Si il un increment de deplacements imposes ----------------- * list f2_inc; zclim_mail = zclim extrai 'MAILLAGE'; * list f2_base; ZFEXT2 = ZFEXT2 + F2_INC + F2_base; ; 'FINS'; *---------------- si chargement deformation actualisation DEFOR ------ 'SI' WTAB.'CHAR_DEFI' ; 'FINS'; *---------- dynamique : preparation du second membre -------------- 'SI' ( WTAB.'DYNAMIQUE' ) ; 'SI' ('EGA' WTAB.'FREA1' 'INCONNU'); 'FINS'; F1 = F1 + F1F ; 'SINON'; F1 = F1F ; 'FINS' ; 'FINS'; 'SI' WTAB.'PROCEDURE_CHARMECA'; TFF1 = CHARMECA PRECED TEMPS0 ; FF1=TFF1 .'ADDI_SECOND'; 'FINS'; F1=F1 + FF1; 'SINON'; F1 = FF1; 'FINS'; 'FINS'; 'SI' LOGPRE ; MOP = WTAB.'MOD_PRE' ; 'SINON' ; 'FINS' ; F1=F1 + FF1; 'SINON'; F1 = FF1; 'FINS'; 'FINS'; * 'SI' IRCON ; LAF0 = LAF0 'ET' 'FINS'; wtab.'PROCEDURE_CHARMECA'); LAF1 = F1 - LAF0 ; 'SINON'; LAF1 = -1 * LAF0; 'FINS'; LAF2 = 'ENLEVER' (ZCLIM0 '*' conti.'DEPLACEMENTS') 'FLX'; * forces exterieures + reactions - forces interieures au debut du calcul * c.a.d. (masse*acceleration initiale)+(amortissement*vitesse initiale) WTAB.'FREA1' = LAF1 - LAF2 ; 'FINS'; 'SI' wtab.'LIAISON_PERSISTANTE' ; * forces d'acceleration au debut du pas 'SI' ('NEG' WTAB.'AMORTISSEMENT' 'INCONNU'); FF4 = WTAB.'AMORTISSEMENT'* conti.'VITESSES'; WTAB.'FMAN'= WTAB.'FREA1' - FF4 ; 'SINON' ; WTAB.'FMAN'= WTAB.'FREA1' ; 'FINS'; 'FINS'; *------------ il faut calculer la matrice de masse tout de suite ------------- WTAB.'GRANDS_DEPLACEMENTS'); ZMAT22 WTAB.'MO_TOT' ) ; 'SI' WTAB.'MASSCONST'; WTAB.'MASSE'=WTAB.'MASSE' 'ET' WTAB.'MASSE_CONSTANTE'; 'FINS'; 'FINS'; FF = WTAB.'MASSE' *conti.'VITESSES'; * partie du second membre qui ne depend que des informations du pas prec WTAB.'FREA1' = FF4 + WTAB.'FREA1'; 'FINS'; *--------- consolidation : preparation du second membre --------- 'SI' WTAB.'CONSOLIDATION' ; FF4 = 'EXCO' WTAB.'MOT_POR' FF WTAB.'MOT_POR' 'NOID' 'NATURE' 'DISCRET' ; 'SI' WTAB.'DYNAMIQUE'; WTAB.'FREA1' = WTAB.'FREA1' + FF4 ; 'SINON'; WTAB.'FREA1' = FF4 ; 'FINS'; * ---- traitement des flux si besoin ---- FACFLU = -1. * WTAB.'DT'; FLUXT = ( FACFLU * (1 - WTAB.'TETA') * FLUXT0 ) + ( FACFLU * WTAB.'TETA' *FLUXTI ) ; ZFEXT2 = ZFEXT2 '+' FLUXT ; 'FINS' ; 'FINS'; *----------------- non-local type HELM : preparation -------------- 'SI' (LNLOC 'ET' ('EGA' WTAB.'NON_LOCAL' 'HELM')); PAS_HELM PRECED ; 'FINS' ; *----------------- calcul de la masse si frequentiel ------------- ZMAT22 WTAB.'MO_TOT' ) ; 'FINS'; *----------------- pilotage indirect : preparation -------------- 'SI' (LOGPIL); ZFPILIN = WTAB.'FORCES_PILOTEES' 'ET' WTAB.'DEPLACEMENTS_PILOTES' ; ETA0 = 'EXTR' PRECED.'COEFFICIENT_DE_PILOTAGE' D_ETA=0.; 'FINS' ; *- Second membre ZFCONSTA = ZFEXT2 ; ************************************************************************ * Parallélisation du GIBIANE via les ASSISTANTS * ************************************************************************ ZMODLI = WTAB.'MO_TOT' ; ZMODLP = WTAB.'MO_TOT_PREC'; NBPART = WTAB.'NBPART' ; PARALLEL = FAUX ; PARTLOCA = FAUX ; ZMODL = ZMODLI ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'COMPORTEMENT') ; PARALLEL = VRAI ; PARTLOCA = VRAI ; 'FINS' ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'AUTOMATIQUE') ; PARALLEL = VRAI ; PARTLOCA = FAUX ; 'FINS' ; * Test sur un MODELE qui aurait changé 'SI' ('NEG' ZMODLI ZMODLP) ; 'SI' IPLAVI ; 'FINS'; 'SI' ITHER ; 'SI' POR1; 'FINS'; 'FINS'; WTAB.'DEFOR1' = 'FINS'; WTAB.'DEFOR2' = 'FINS'; 'FINS'; *-----------materiau au debut du pas ------------------------------ ZMAT11 = WTAB.'MAT1' ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'AUTOMATIQUE') ; 'SINO'; ZMAT1 = ZMAT11; 'FINS'; ************************************************************************ * Quelques initialisations * ************************************************************************ STAB12.'ZU1' = CONTI.'DEPLACEMENTS' ; STAB12.'SIGF' = CONTI.'CONTRAINTES' ; STAB12.'DEFF' = CONTI.'DEFORMATIONS' ; STAB12.'FNONL' = WTAB.'FNONL' ; STAB12.'RESIDU' = WTAB.'RESIDU' ; STAB12.'XDENO' = WTAB.'XDENO' ; STAB12.'XDENOM' = WTAB.'XDENOM' ; STAB12.'ETAT1' = WTAB.'ETAT1'; 'FINS'; STAB12.'ETAT2'= ZETAT2; STAB12.'DEFOR1' = WTAB.'DEFOR1' ; STAB12.'DEFOR2' = WTAB.'DEFOR2' ; 'FINS'; STAB12.'FNONL' = WTAB.'FNONL' ; 'FINS'; STAB12.'TET1' = WTAB.'TET1' ; STAB12.'TET2' = WTAB.'TET2' ; 'FINS'; * STAB12.'SUCCES' = VRAI ; 'SI' ('NEG' WTAB.'AUTOCOEF' 'INCONNU') ; STAB12.'AUTOCOEF' = WTAB.'AUTOCOEF' ; 'FINS' ; 'SI' ('NEG' WTAB.'AUTOREDU' 'INCONNU') ; STAB12.'AUTOREDU' = WTAB.'AUTOREDU' ; 'FINS' ; 'SI' ('NEG' WTAB.'SECOND_MEMBRE' 'INCONNU') ; STAB12.'SECOND_MEMBRE' = WTAB.'SECOND_MEMBRE' ; 'FINS' ; 'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ; STAB12.'LASTKTAN' = WTAB.'LASTKTAN' ; 'FINS' ; 'SI' ('NEG' WTAB.'AUTORED1' 'INCONNU') ; STAB12.'AUTORED1' = WTAB.'AUTORED1' ; 'FINS' ; 'SI' ('NEG' WTAB.'LISEA_M' 'INCONNU') ; STAB12.'LISEA_M' = WTAB.'LISEA_M' ; STAB12.'RIBLO_M' = WTAB.'RIBLO_M' ; 'FINS' ; 'SI' ('NEG' WTAB.'INCREMENT' 'INCONNU' ); STAB12.'INCREMENT' = WTAB.'INCREMENT'; INCRPREC = STAB12.'INCREMENT' ; 'FINS'; STAB12.'FFROT' = WTAB.'FFROT' ; STAB12.'INITEMPS' = WTAB.'INITEMPS' ; STAB12.'DT' = WTAB.'DT' ; STAB12.'DTPREC' = WTAB.'DTPREC' ; 'SI' ITHER ; STAB12.'TETA1' = WTAB.'TET1'; STAB12.'TETA2' = WTAB.'TET2'; 'SI' POR1; 'FINS'; 'SINON' ; 'FINS' ; * SP : initialisation DFGRAD en presence d'un modele MECANIQUE * (DFGRAD mis a INCONNU par defaut dans PAS_INIT) 'SI' ('NEG' WTAB.'MOD_MEC' 'INCONNU') ; 'SI' ('NEG' WTAB.'DFGRAD' 'INCONNU') ; **'SINO'; ** STAB12.'DFGRAD' ='REDU' WTAB.'DFGRAD' ZMODL; 'FINS' ; 'SINON' ; STAB12.'DFGRAD' = WTAB.'DFGRAD' ; 'FINS' ; HPP_EPS = FAUX ; * Option a n'utiliser que par les utilisateurs avertis III = PRECED.'ACCELERATION' ; 'FINS' ; * Matrice tangente : non utilisee si IPLAVI a FAUX IKTAN = IKTAN 'ET' IPLAVI ; 'SI' (WTAB.'K_TANGENT' 'ET' ('NON' IPLAVI)) ; ' on utilise la rigidite elastique' ; 'FINS' ; * Matrice tangente par perturbation : * Option non disponible si non local ou si IPLAVI a FAUX IPERT = WTAB.'K_TANGENT_PERT' 'ET' ('NON' LNLOC) 'ET' IPLAVI ; ZPERC1 = WTAB.'K_TANG_PERT_C1' ; ZPERC2 = WTAB.'K_TANG_PERT_C2' ; * Matrice tangente : partie symetrique utilisee 'SI' WTAB.'K_TANGENT_SYME' ; 'SINON' ; ZKTASYM = 'TEXTE' ' ' ; 'FINS'; * Matrice tangente : pas d'acceleration en cas de modele FEFP ou SSTE 'SI' IKTAN ; 'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS'; 'FINS' ; * 'SI' IFTOL ; ZFTOL = 'ABS' WTAB.'FTOL' ; 'FINS'; 'SI' IMTOL ; ZMTOL = 'ABS' WTAB.'MTOL' ; 'FINS'; * * CB215821 : Devrait t-on mutualiser ITCAR et WTAB.ITCAR? ==> * Ce ne sont pas tout a fait les memes tests aujourd'hui 'SI' (ITCAR 'EGA' FAUX); ZMAT2 = ZMATI ; ZMAT2R= ZMAT ; 'FINS'; 'SI' WTAB.'ITCAR'; 'FINS'; CARA1 = ZMAT1; 'FINS'; 'SI' ('OU' ('OU' IVISCO IVIDOM) IVIEXT); ZPREK = 5.E-7 ; 'FINS'; 'SI' IENDOM; ZPREK = ZPREC ; 'FINS'; * on fait ici la séparation poreux .. pour l'avoir sur *les modèles partitionnes 'SI' POR1; ** kich ma_por0 intervient si ISOL ** initialement MA_POR0= 'REDU' MO_POR (STAB12.'MAT1'); 'FINS'; * recuperation de certains champs, si nbpart>1 zmodl est partitionné * sinon c'est le modele initial * * 'SI' IPLAVI ; lnom = com_var ; 'SI' ISOL ; 'FINS' ; 'FINS' ; * * teste t'on les moments ? * teste t'on les POREUX ? 'SI' POR1 ; TSTMOM=VRAI ; 'FINS'; * IKLFFF=VRAI; 'SI'TSTMOM; 'SI' IFTOL; 'SI' IMTOL; IKLFFF=FAUX; 'FINS'; 'FINS'; 'FINS'; * 'SI'('NON' TSTMOM); 'SI' IFTOL; IKLFFF=FAUX; 'FINS'; 'FINS'; * GEOREF0 = WTAB.'FOR0' ; ISOUSPPP=0; WTAB . 'ISOUSPAS' = 0; NSOUSPAS = WTAB . 'MAXSOUSPAS'; ZCCONV = VRAI ; KNOCONV = 0 ; * Initialisation CHArgement SANS T : RED_URG = 0; 'SI' autaug; iraug = faux; 'FINSI'; ************************************************************************ ****** boucle de non convergence ************************************************************************ resmul = 1; augmult = 0.60000000; 'REPETER' BONOCONV NSOUSPAS ; 'SI' ('EGA' RED_URG 0); 'SI' ((XCONV < ZPREC) 'ET' (DEPSTDM < ZPREC) 'ET' (resmul > 0.99) 'ET' (AUGMULT < 100.)); augmult = 0.60000000; *pv resmul = 1; augauto = augmult; augk = augmult; 'SI' autaug ; iraug = faux; 'FINSI'; znacce=VRAI; 'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS'; 'SINON'; augmult = augmult / 1.5; ** si (augmult < 0.6); augmult = 0.6; finsi; 'FINSI'; 'SINON'; augmult = augmult * 1.5; 'SI' autaug; IRAUG = VRAI; 'FINSI'; 'FINSI'; *** wtab.'RECARI' = vrai; KNOCONV = KNOCONV+1 ; DT_INIT = STAB12.'DT' ; DTINI = STAB12.'DT' ; ZSIG0 = ZSIGF ; zsig0_1 = zsig0; 'SI' IPLAVI ; ZEPS0 = ZDEIF ; ZVAR0 = ZVARF ; 'FINS'; ZU1 = STAB12.'ZU1' ; GR_U_DEB = STAB12.'DFGRAD'; 'FINSI'; 'SI' ITHER ; TETA1 = STAB12.'TETA1' ; TETA2 = STAB12.'TETA2' ; DTETD = TETA2 '-' TETA1; 'FINS' ; * materiau au debut du pas en cas de non convergence etat1=etat2 'SI' (knoconv > 1) ; 'SINON'; ZMAT1 = ZMAT1I ; 'FINS' ; *------ caracteristiques initiales en cas de grands deplacements ------ 'SI' WTAB.'ITCAR'; ZMAT1 = ZMAT1 '-' MECAR1 '+' MECAR2 ; 'FINS'; *----------- Calcul du champ de materiau a la fin du pas ------------ 'SI' (knoconv '>' 1) ; 'SINON' ; MMMM = ZMAT22 ; 'FINS' ; *----- Caracteristiques initiales en cas de grands deplacements ----- 'SI' WTAB.'ITCAR'; MMMM = MMMM '-' MECAR1 '+' MECAR2 ; * (fdp) on reporte les variations du materiau sur le pas * dans les autres instances du champ materiau * (ZMAT et ZMATI semblent suffirent) ZMAT = ZMAT '-' MECAR1 '+' MECAR2 ; 'FINS'; *------------ Calcul de la rigidite a la fin du pas ---------------- * Test ci-dessous est inutile, WTAB.'BLOCAGES_MECANIQUES' et ZCLIM0 * correspondent aux memes rigidites ... (cf. ligne 117) 'SI' ( (WTAB.'RECARI' 'OU' IRAUG ) 'OU' 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); 'DETR' HOOKENDO; 'SINON'; 'SI' ((LAG_TOT 'EGA' 1) 'ET' vrai); 'FORM' GEOREF0; 'FORM' GEOM1; 'DETRUI' HOOKRH; 'DETRUI' HOOKRH2; 'SINON'; 'FINSI'; 'FINS'; * RH peut contenir des CL ZCLIM = ZCLIM0 'ET' ZCL ; 'FINSI' ; RRRR = RH 'ET' ZCLIM0 ; * Prise en compte d'eventuelles RIGIDITE_CONSTANTE 'SI' IRCON; RRRR = RRRR 'ET' RIG_CONS; 'FINS'; * ZRAIDNA = RRRR; 'SI' (IRAUG 'ET' ('NON' AUTAUG)); RRRR = RRRR 'ET' RIG_AUG ; 'FINSI'; * 'SI' (AUTAUG 'ET' ('NON' IRAUGLU)); RIG_AUG = 'MASSE' ZMODL MMMM ; *** mess 'actualisation rig_aug'; 'FINSI'; 'SI' (AUTAUG 'ET' IRAUG); * recalcul masse pour l'augmentation RRRR = RRRR 'ET' (RIG_AUG * augauto) 'ET' (RH * augk) ; 'FINS'; 'SI' ('NON' IRAUG); * recalcul masse pour l'augmentation 'SI' IRAUGLU; *** 'MESS' ' augmentation residuelle '; RRRR = RRRR 'ET' RIG_AUG ; 'FINS'; 'FINSI'; * Stockage de la rigidite pour eviter de la recalculer WTAB.'RRRR'=RRRR; 'FINSI'; * 'SINON'; RRRR=WTAB.'RRRR'; 'FINS'; ZRAID=RRRR; *------------ consolidation ou dynamique faut-il recalculer l'operateur? * -----: preparation du pas de temps ------------ 'SI' ( WTAB.'CONSOLIDATION' 'OU' WTAB.'DYNAMIQUE'); DT = TI '-' TEMP0; 'SI' ( '>' (DELTAN '*' 0.9999) DT) ; WTAB . 'RECAOP' = VRAI; 'FINS'; 'SI' ( '<' (DELTAN '*' 1.0001) DT) ; WTAB . 'RECAOP' = VRAI; 'FINS'; 'SI' WTAB.'MATVAR'; WTAB . 'RECAOP' = VRAI; 'FINS'; DELTAN=WTAB.'DT'; 'FINS'; *---------------------- Formation de l operateur ----------------------- 'SI' ( WTAB.'DYNAMIQUE' 'OU' WTAB.'CONSOLIDATION'); 'SI' ('NEG' WTAB.'OPERATEUR' 'INCONNU'); ZRAID = WTAB.'OPERATEUR'; 'SI' (WTAB . 'RECAOP') ; ZRAID = RRRR ; 'FINS'; 'FINS'; 'FINS'; *------------ operateur frequentiel ----------------------- *------------ operateur amortissement en frequentiel 'SI' WTAB.'FREQUENTIEL' ; RR2 = 'CHAN' 'INCO' RRR2 RR3 = 'CHAN' 'INCO' RRR2 RRR2 = RR2 'ET' RR3 ; RR1 = ZRAID ; OMEGI= 2.* PI * TI ; RRR1 = OMEGI * OMEGI * (-1.) * WTAB.'MASSE' ; RR1 = ZRAID 'ET' RRR1 ; ZRAID= RR1 'ET' RR4 ; RR5 = OMEGI '*' RRR2 ; ZRAID= ZRAID 'ET' RR5 ; 'FINS' ; *--------------- et la perméabilité ---------------------------------- 'SI' (WTAB.'CONSOLIDATION') ; 'SI' (WTAB.'GRANDS_DEPLACEMENTS' 'OU' WTAB.'MATVAR') ; WTAB .'RECAOP' = VRAI ; 'FINS'; 'FINS'; *------------- Cas de la consolidation ou de la dynamique ------------- *------------- il faut recalculer l'operateur d'iteration ------------- 'SI' (WTAB . 'RECAOP') ; 'SI' ( WTAB.'DYNAMIQUE'); ZRAID = 4.D0 '/'( DT ** 2) '*' WTAB.'MASSE' 'ET' ZRAID; 'SI' ('NEG' WTAB.'AMORTISSEMENT' 'INCONNU'); ZRAID = WTAB.'AMORTISSEMENT' '*' (2.D0 '/' DT) 'ET' ZRAID; 'FINS'; 'FINS' ; 'SI' ( WTAB.'CONSOLIDATION'); ZRAID =-1.* DT* WTAB.'TETA'* WTAB.'PERMEABILITE' 'ET' ZRAID ; 'FINS' ; WTAB . 'OPERATEUR'= ZRAID ; WTAB . 'RECAOP' = FAUX ; 'FINS'; ZRAIDNA = ZRAID; *-------------- traitement des contacts frottements automatiques ------- CDEP = STAB12.'ZU1' ; BFCONT = FAUX ; BCLIM2 = FAUX ; 'SI' WTAB.'CONTACT'; MODCON= WTAB.'MODCONTA'; ** list resu wtab.'MATCONTA'; 'SINON' ; 'FINS'; ** list modfro; ** list matfro; * 'SI' ('NEG' 0 CJEU) ; 'FINSI' ; 'SI' (WTAB . 'MODAL') ; 'NATURE' 'DISCRETE' ; 'SI' ('EGA' 1 &BCDA) ; CCDA = CHCR ; 'SINON' ; CCDA = CHCR 'ET' CCDA ; 'FINS' ; 'FIN' BCDA ; CDAP = CCDA ; 'FINS' ; * ATTENTION : mettre les conditions de frottement en premier pour * les numeroter en dernier 'SI' ('NEG' CRR 0); BCLIM2 = VRAI ; BFCONT = VRAI ; ZCLIM2 = CRR 'ET' ZCLIM2 ; CCOR = CRR '*' CDEPSLX ; ZFCONT = ZFCONT '+' CCOR ; 'FINS'; 'SI' ( 'NEG' 0 CJEU); BFCONT = VRAI ; ZFCONT = ZFCONT '+' CDAP ; 'FINS'; 'SI' ('NEG' 0 RFROT); BCLIM2 = VRAI ; BFCONT = VRAI ; ZCLIM2 = RFROT 'ET' ZCLIM2 ; CCOR = RFROT '*' CDEPSLX ; ZFCONT = ZFCONT '+' CCOR ; 'FINS'; * * Mise a jour de ZCLIM et ZFCONSTA si necessaire 'SI' BCLIM2; ZCLIM = ZCLIM2 'ET' ZCLIM ; ZRAID = ZCLIM2 'ET' ZRAID; 'FINS'; * 'SI' BFCONT; ZFCONSTA = ZFEXT2 '+' ZFCONT ; 'FINS'; 'FINS'; * ZFCONST1 = ZFCONSTA; zclimp = zclim; * --------------- pilotage automatique ******************************** ISNPB = FAUX ; AL1 = 1. ;COEPI = 1.d0; COEINC=0.d0;COEPI0=1.d0;DAL1=100.D0; CORPREC = 1. ; * CORPREC = 10.; EPS_EPS = 'TEXTE' 'LINEAIRE'; HPP_EPS = VRAI; 'FINS'; * 'SI' IPILOT ; 'SI' ( WTAB.'AUTODEUX' ) ; COEPI = 'ABS' ( STAB12.'AUTOCOEF'); COEPI = COEPI / (1.-COEPI); 'SI' (COEPI > 1.D0) ; COEPI=1.D0;'FINS'; COEPI0=COEPI;STAB12.'AUTOCOEF'=COEPI; 'SINON'; STAB12.'AUTOCOEF' = 1.D0; 'FINS' ; RED1 = 1. ; RED2 = 0 ; 'SI' ('NEG' WTAB.'AUTORED1' 'INCONNU') ; 'SI' (STAB12.'AUTORED1' > 0); STAB12.'AUTORED1' = STAB12.'AUTORED1' - 1; 'SI' (STAB12.'AUTORED1' 'EGA' 0) ; COEPI = 3 * COEPI ; STAB12.'AUTOREDU' = STAB12.'AUTOREDU' / 3.; 'SI' (STAB12.'AUTOREDU' > 1.1 ); * on travaille encore avec un critere reduit STAB12.'AUTORED1' = 4 ; RED1 = 3. ; 'FINS'; 'FINS'; 'FINS'; 'FINS'; 'SI' (COEPI > 1d0); RED1 = RED1 / COEPI ;COEPI =1d0; 'FINS'; 'SI' ( 'NEG' WTAB.'NBPLAS' 'INCONNU') ; 'SI' ( WTAB.'NBPLAS' 'EGA' 0) ; 'SI' (COEPI < 0.) ;COEPI = COEPI * -2.;'FINS'; 'FINS'; 'FINS'; COEPI = 'ABS' COEPI ; COEPI0 = COEPI; * sans pilotage 'SINON'; STAB12.'AUTOCOEF'= 1.D0; 'FINS'; ************************************************************************ *------------- quelques initialisations pour la boucle ETIQ ******** ************************************************************************ URG = FAUX; RED_URG = 0 ; IT = 0 ; c_zdepr = faux;ITACC = 0; ZICONV = VRAI; MMC = 0 ; MMCMAX = 0 ; EPSM = 0.; DPSMAX = xpetit; dpsmaxp = dpsmax; DEPSTDM = 0. ; DEKREAC1 = 0. ;XCONVNOR = 0. ; ITNORM1 = 0 ; DEPSTREF = 100. * (WTAB . 'MAXDEFOR') ; GR_U_K = GR_U_DEB ; DITNORM1 = 0 ;NBCYCLE1 = 0 ; zdeptq = zdept; zdeptp = zdept; zprecnc=1e-5; FDEF = 0. ; ITNV = -5 ; 'SI' IFEFPUL; XUPDA = 1; 'FINS'; *********** Corrections en cas de materiaux variables **************************** DEPST0=0 ; 'SI' WTAB.'MATVAR' ; 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); 'SINON'; 'FINS'; DDEF0 = XXX4 - XXX3; 'DETR' XXX3; 'DETR' XXX4; DEPST0=-1.* DDEF0; 'FINS'; ***** En cas de chargement thermiques ********************************* DMSRT0=0; 'SI' ITHER ; 'SI' (WTAB.'MATVAR' 'ET' IPILOT); 'MESS' 'Le pilotage n est pas possible avec un materiau qui depend de la temperature' ; 'ERREUR' 19 ; 'FINS'; 'SINON'; MCHTETA2 = TETA2; 'FINS'; DTT = ETT '-' ETT0 ; 'SI' WTAB.'POR1'; *------ Cas du milieu poreux avec chargement thermique ---------- * cas isotrope seulement pour le moment * et on ne s'occupe pas du alpha-reference !! DMSRT0= MSRTT '-' MSRTT0 ; 'FINS' ; DEPST0 = DTT '+' DEPST0 ; 'FINS'; * calcul de DSIGT0 et FTHE en cas de chargement si necessaire *********** FTHE = 0 ; 'SI' (DEPST0 'NEG' 0); 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); 'SINON'; 'FINS' ; 'SI' (DMSRT0 'NEG' 0); DSIGT0 = DSIGT0 '+' DMSRT0 ; 'FINS'; 'FINS'; *-------------- deformations imposes ==> contraintes imposees ********************************** 'SI' LOGDEF; 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); 'SINON'; 'FINS'; 'FINS'; *-------------- pression imposee et grands deplacements ?*************** 'SI' (LOGPRE 'ET' IGRD) ; MOP = WTAB.'MOD_PRE' ; * 'SINON' ; 'FINS' ; ZFP0F = 'COPIER' ZFPEXTF ; COEFP=1.D0; 'FINS'; * --------- ktangent et fefp****************************************** 'SI' IKTAN ; 'SI' IFEFP ; IKT_SAUV = VRAI ; 'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ; 'MESS' 'FEFP: Start with LASTKTAN' ; ZRIKTA = STAB12.'LASTKTAN' ; ZRAID = ZCLIM 'ET' ZRIKTA ; 'SINON' ; ZRAID = ZRAID 'ET' ('KSIGMA' ZMODL ZSIG0 ZMAT) ; 'FINS' ; 'SINON' ; 'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ; IKT_SAUV = VRAI ; 'SI' IPERT ; 'MESS' 'Matrice tangente par perturbation - ' 'SINON' ; 'MESS' 'Matrice tangente "coherente" - ' 'FINS' ; ZRIKTA = STAB12.'LASTKTAN' ; ZRAID = ZCLIM 'ET' ZRIKTA ; 'SINON' ; 'SI' IPERT ; IKT_SAUV='NEG' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE'; 'MESS' 'Matrice tangente par perturbation - ' 'SINON' ; IKT_SAUV = ('NEG' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE') 'ET' ('NEG' WTAB.'K_TANGENT_ITER0' 'MAT_TANGENTE') ; 'SI' ('EGA' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE') ; 'MESS' 'Matrice tangente "coherente" - ' 'SINON' ; 'MESS' 'Matrice tangente "coherente" - ' DTTAN = 0. ; ZRIKTA = 'KTAN' ZMODL ZSIG0 ZVAR0 ZMAT 'PREC' ZPREK 'DT ' DTTAN ZKTASYM ; ZRAID = ZCLIM 'ET' ZRIKTA ; 'FINS' ; 'FINS' ; 'FINS' ; 'FINS' ; 'FINS' ; *-------- en grands deplacements option K_SIGMA *********************** * 'SI' (IGRD 'ET' ('NON' HPP_EPS)); 'SI' (IKSIA 'ET' ('NON' IFEFP)) ; KSI1 ='KSIGMA' ZMODL ZSIG0 ZMAT; ZRAIDINI= ZRAID ; ZRAID = ZRAID 'ET' KSI1 ; 'FINS' ; 'FINS' ; * Y a-t-il des forces non conservatives ( forces suiveuses)? *********** ADDISEC0 = FAUX; ADDISEC2 = FAUX; 'SI' WTAB.'PROCEDURE_CHARMECA'; * on ajoute l indice ADDI_MATRICE pour signaler a charmeca qu on * souhaite aussi l operateur linearisé des Forces NL de charmeca PRECED . 'ADDI_MATRICE' = vrai; TFP22 = CHARMECA PRECED WTAB.'T_FINAL'; PRECED . 'ADDI_MATRICE' = faux; * FP22 = F^suiv_n+1 DMZPRES = 0.D0 ; FP22 = TFP22.'ADDI_SECOND' ; FP022 = 'COPIER' FP22 ; DMZPRES = MZPRES ; ADDISEC2= VRAI ; 'FINS'; ZRAID = ZRAID 'ET' TFP22.'ADDI_MATRICE'; 'FINS'; * FP0 = F^suiv_n TFP0 = CHARMECA PRECED TEMPS0 ; FP0 = TFP0.'ADDI_SECOND' ; DMZPRES = DMZPRES '-' MZPRES0 ; ADDISEC0= VRAI ; 'FINS'; * 'SI' ('EXIS' TFP0 'ADDI_MATRICE'); * ZRAID = ZRAID 'ET' TFP0.'ADDI_MATRICE'; * 'FINS'; COEFP = 1.D0 ; 'FINS'; * -----------calcul de la partie constante du second membre ********** * en consolidation * ZFP1 est cense contenir : - B0*SIG0 et * DT*(1-TETA)*FI0 + DT*H*P * * dans ZFCONSTA on met le second membre de u *********************** * en dynamique ******************************************************** * ZFP1 est cense contenir : F0 + 4/DT*M*V0 - B0*SIG0 'SI' IDYN ; UNSURH = 1.D0 '/' STAB12.'DT' ; ZFP1 = WTAB.'FREA1' ; ZDYFEXT = ZFCONSTA 'ENLEVER' 'FLX'; ZFCONSTA= ZFCONSTA '+' ZFP1 ; 'FINS'; ZFEXT = ZFCONSTA 'ENLEVER' 'FLX'; *---------deplacement (ou jeu) e imposer e la fin du pas ************* * on separe les efforts ZFEXT (=F^ext_n+1) *************************** * et deplacement (ou jeu) ZFLX1 (u^imp_n+1) a imposer a la fin du pas 'SI' ('NEG' STAB12.'FFROT' 'INCONNU'); FFROT = STAB12.'FFROT'; 'SINON' ; FFROT=ZFEXT * 0; 'FINS'; FFROTP = FFROT ; * calcul des forces externes deja equilibrees au debut du pas ******** * par B*SIGMA : ZF1 = F^int_n = B*sigma_n + K^cst*u_n 'SI' IRCON; 'FINS'; 'SI' IDYN ; FFDYN = 'COPIER' ZF1; 'FINS'; 'SI' ISOL ; XXXS =((1.- WTAB. 'TETA' )*GRAP0)+ (WTAB. 'TETA' * XXX1); XXX3 = ZF1 ; 'DETR' XXX2 ; 'FINS'; * initialisation des variables forces et deplacement******************* * zzd est le deplacement au pas precedent (=u_n) et ZLX=lambda_n * --------- flxini est la partie des FLX deja realisee au debut du pas **pv FLXINI= ZZD '*' ZCLIM; * FREAP : -1*reactions du pas precedent (=F^reac_n) transportees sur les nou **pv FREAP = ZLX '*' ZCLIM; FZU1 = ZU1 '*' ZCLIM; FEXT0 = ZF1 '+' FREAP; * FEXT0 est le chargement externe (sans reactions vu par la * structure le pas d'avant) = F^ext_n = F^int_n - F^reac_n * * on va calculer le premier residu c'est a dire le desequilibre ******* * entre les forces externes et le calcul B*SIGMA. * le sigma qui sert est celui qui existerait si le champ de * deplacement ne changeait pas (ZU1). faire attention aux FLX * En pilotage on reprend ce residu que l'on multiplie par COEPI * XXX1 = [F^ext_n+1 ; Du^imp] * DFEXT0 = increment des forces et des FLX a imposer (le residu * du pas precedent) = [DF^ext ; Du^imp] XXX1 = ZFCONSTA '-' FLXINI ; DFEXT0 = XXX1 '-' FEXT0 ; * si pression suiveuse dfext0 contient en plus l'increment des forces * de pression du uniquement a la reactualisation de la geometrie (sans * augmentation du module) * mais comme F^int_n equilibre deja F^ext_n + F^suiv_n +... et qu'on est * toujours sur config_n, on doit avoir : DFEXT0= [DF^ext ; Du^imp] * avec DF^ext qui ne contient pas de forces suiveuses (...a vérifier) 'SI' ADDISEC0; DFEXT0 = DFEXT0 '+' FP0 ; 'FINS'; 'SI' (LOGPRE 'ET' IGRD) ; DFEXT0 = DFEXT0 '+' ZFPEXT0 ; 'FINS' ; DFEXT0F = DFEXT0 'ENLEVER' 'FLX'; * DFEXT0 = [DFEXT0F ; DFEXT0L] * = [F^ext_n+1 - (F^int_n - F^reac_n - F^suiv_n) ; Du^imp] * RESIDU = forces exterieures sans reactions * (avec des termes supplementaires le cas echeant p.ex. en * dynamique ou en poreux) - forces interieures * et increment des relations imposees * ---> la resolution fournira dU et dR * RESIDU = [F^ext_n+1 + F^ther_n+1 + F^defi_n+1 - F^int_n ; Du^imp] RESIDU = XXX1 '+'FTHE '+'FDEF '-'ZF1; ZDFINI ='COPIER' DFEXT0 ; ZFPLO = ZF1 '-' FTHE '-' FDEF ; 'SI' (ITHER 'OU' WTAB.'MATVAR'); ZDFINI=ZDFINI '+' FTHE; 'FINS'; 'SI' LOGDEF; ZDFINI= ZDFINI '+' FDEF; 'FINS'; 'SI' (LOGPRE 'ET' IGRD) ; RESIDU = RESIDU '+' ZFPEXTF ; ZDFINI = ZDFINI '+' ZFPEXTF '-' ZFPEXT0 ; 'FINS'; 'SI' ADDISEC0 ; ZDFINI = ZDFINI '-' FP0 ; 'FINS'; 'SI' ADDISEC2 ; * mess ' fp22 ' ; list resu fp22; RESIDU = RESIDU '+' FP22 ; * on tient compte de l'increment des forces suiveuses en direction * et en module ZDFINI = ZDFINI '+' FP22 ; 'FINS'; * ici, on a : * RESIDU = [F^ext_n+1 + F^ther_n+1 + F^defi_n+1 + F^suiv_n+1 * - F^int_n ; Du^imp] * = [ DF^tot ; Du^imp] * ZDFINI = [F^ext_n+1 + F^ther_n+1 + F^defi_n+1 + F^suiv_n+1 * - (F^int_n - F^reac_n) ; Du^imp] * IMPO12= FAUX; IMPO12U = FAUX; IMPO12 = VRAI; IMPO12U = VRAI; DIMPOV = DFEXT0L '-' DIMPO12; 'FINS'; stab12.'SECOND_MEMBRE' = RESIDU '*' 1.D0; ************************************************************************ *** 1ERE RESOLUTION *** ************************************************************************ 'SI' (WTAB.'ADHERENCE') ; 'SINON' ; FADRE = FADHE ; 'FINSI' ; FEXCI = FEXCI 'ET' FADRE ; 'FINSI' ; * 'SI' (WTAB.'CAFROTTE') ; FEXCI = FEXCI 'ET' FFROT ; 'FINSI' ; * 'SI' LOGPIL ; * Modification du residu RESIDU = RESIDU 'ET' (ETA0 '*' ZFPILIN); 'FINSI' ; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12); ZDEP1 BID BID BID BID = 'RESO' ZRAID 'SOUC' RESIDU 'INIB' STAB12.'RIBLO_M' STAB12.'LISEA_M' FEXCI ; 'SINON'; ZDEP1 BID BID BID = 'RESO' ZRAID 'SOUC' RESIDU 'INIB' STAB12.'RIBLO_M' STAB12.'LISEA_M' ; 'FINS'; 'OUBLIER' STAB12 'RIBLO_M'; 'SINON'; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12U); 'SINON'; 'FINS'; 'FINS'; 'SI' LOGPIL ; ZDEPILO = 0. * ZDEPII ; * Calcul de D_eta par appel a la procedure PILOINDI D_ETA = PILOINDI PRECED ZU1 ZDEPILO ZDEP1 ZDEPII DTAU; * Mise à jour ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ; ETA = ETA0 ; 'FINSI' ; 'SI' IDYN; STAB12.'ZRAIDV'= ZRAID; 'FINS'; STAB12.'RIBLO_M' = ZRAID_T. 7 ; STAB12.'LISEA_M' = ZRAID_T. 6 ; 'FINS'; * pour eventuel calcul de augauto * et initialisation zdepl (sert pour xnum) * On sauve le deplacement initial pour la convergence forcee *********** 'SINON'; zdeptini = zdep1; 'FINS'; * * calcul d'une norme pour la convergence************************** * XXX1= ZFEXT; 'SI' ADDISEC2; XXX1=ZFEXT + FP22; 'FINS'; 'SI' (LOGPRE 'ET' IGRD) ; XXX1 = ZFEXT + ZFPEXTF ; 'FINS'; 'SI' LOGPIL ; * Modification du chargement XXX1 = XXX1 '+' ((ETA '+' D_ETA) '*' ZFPILIN); 'FINSI' ; ZDEP1P50 = ZDEP1 ; 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET')) MLPRIM MLDUAL; 'SI' ('VERI' xdeno) ; 'SINON'; xdeno = 1; 'FINSI'; XDENO1 = 'ABS' XDENO + (MZFM * MZDEP1M); MZDEP1M = MZDEP1M + XPETIT; XDENO=XDENO1/MZDEP1M; XDENO = XDENO + MZFM; XDENO = XDENO + XPETIT ; XDENOM=XDENO; 'SI' TSTMOM ; ZDEP1P50 = ZDEP1 + XPETIT; ; XDENOM = XDENOM + XPETIT ; 'FINS' ; 'SI' IPILOT ; 'SI' ( WTAB.'AUTODEUX' ) ; XDENO = STAB12.'XDENO'; XDENOM = STAB12.'XDENOM'; 'FINS'; 'FINS'; * 'MESS' 'non convergence on garde xdeno xdenom ' xdeno xdenom; 'SI' (XDENO < STAB12.'XDENO'); XDENO = STAB12.'XDENO'; 'FINSI'; 'SI' (XDENOM < STAB12.'XDENOM'); XDENOM = STAB12.'XDENOM'; 'FINSI'; 'FINS'; IAFAIR=FAUX; STAB12.'INCREMENT' = 'COPIER' ZDFINI ; 'FINS'; RESIDNOR = 'COPIER' RESIDU ; 'SI' IPILOT; XXX3=DFEXT0F * ( 1-COEPI) ; XXX2 = 1.D0 -COEPI * FREAP;RESIDU = XXX1 + XXX2; 'SI' IMPO12; RESIDU = 1.D0 - COEPI * DIMPO12 + RESIDU; 'FINS'; IAFAIR=VRAI; 'FINS'; * * petite correction du residu pour esperer gagner du temps ************ * INIT = FAUX ; 'SI'(('NEG' STAB12.'FNONL' 'INCONNU') 'ET' (WTAB.'INITIALISATION')); 'SI' IPILOT; 'SI' ( WTAB.'AUTODEUX' 'ET' (COEPI 'NEG' 1.D0)) ; 'MESS' 'Initialisation a partir du pas precedent ' COEPI; IAFAIR=VRAI; INIT = VRAI; XXX1= RED1 * STAB12.'FNONL'; XXX2= XXX1 + RESIDU ; 'FINS' ; 'SINON'; * on fait la correction si le pas precedent a converge sans non convergence * on fait la correction si le pas precedent etait non lineaire. * on enleve le residu du pas precedent pour recuperer l'increment * nominal du second membre e imposer f1. * f2 et l'increment du second membre du pas precedent STAB12.'INCREMENT'=STAB12.'INCREMENT' - STAB12.'RESIDU'; zdeps = WTAB.'ZDEP1' + zdep1; F1 = STAB12.'INCREMENT' ; F2 = INCRPREC ; AMPL = f12/(f22 + XPETIT); DTPREC=1E30; AMPLT = 0; 'SI' ('NEG' STAB12.'DTPREC' 'INCONNU'); DTPREC = STAB12.'DTPREC'; 'SI' (DTPREC > XPETIT); AMPLT=WTAB.'DT_INIT' /DTPREC; 'FINS'; 'FINS'; * Le chargement n'est il pas de fluage ou de thermique ?'; 'SI' (('ABS' XDCOMP) < (ZPREC * XDENO * mzdep1m) 'OU' (('ABS' FFNO) > (('ABS' F22) * 2.e2 ))); DTPREC = STAB12.'DTPREC'; AMPL=AMPLT; * la decharge est-elle significative 'SI'((F12/(f22 + XPETIT)) < -0.05);ampl=0.;'FINS'; STAB12.'INITEMPS'=VRAI;LOGTEMP=FAUX;; 'MESS' 'Pas d increment de charge, initialisation calculee avec le temps'; 'SINON'; AMPL = F12 / (F22 + XPETIT); LOGTEMP=STAB12.'INITEMPS'; STAB12.'INITEMPS'=FAUX; 'FINS'; * changement de modele on n'initialise pas 'SI' ('NEG' ZMODLI ZMODLP); AMPL = 0; 'FINSI'; 'SI' ((AMPL > 0) 'ET' (AMPL < 2e1)) ; 'MESS' 'Initialisation a partir de la solution precedente Coeff' AMPL; XXX1 = AMPL * STAB12.'FNONL'; XXX2 =RESIDU+ XXX1; RESIDU = XXX2; IAFAIR=VRAI; INIT = VRAI; 'FINS'; 'FINS'; 'FINS'; 'FINS'; 'FINS'; WTAB.'ZDEP1'=zdep1; * * initialisation en plastique ***************************************** * 'SI' IPLAVI ; 'SI' ('NEG' WTAB.'MOVA' 'RIEN') ; 'FINS' ; 'FINS' ; * * debut des iterations internes boucle etiquette ******************** * IKT = FAUX ; IPREM = VRAI ; RECA_K = FAUX; RECA_N = 0; DEPSTP = depst ; DEPSTK = DEPSTP ; * * initialisation acceleration de convergence ******************** * iafair = vrai; PASTEST = FAUX; ZITAC= 0 ; * on peut mettre n'importe quoi c'est pour * ne pas faire de tests dans la boucle ACFP1 = 'COPIER' ZFEXT2 *0. ; ACFP2 = ACFP1 ; ZDEPLD = ACFP1 ; ACFP3 = ACFP1 ; ACFEP1 = ACFP1 ; ACFEP2 = ACFP1 ; FCORF = 'COPIER' FREAP; * initialisation du meilleur critere XCONVMIN = 1e20; DPSMREF = 0 ; XCONV = 0. ; XCONVP = 1. ; NSOINCR = 1 ; 'SI' ('NEG' WTAB.'SOUS_INCREMENT' 'INCONNU') ; NSOINCR = WTAB.'SOUS_INCREMENT' ; 'FINS'; NONCONV = FAUX; CORREC=0; PASREINI=VRAI; coefmt=1.; * * initialisation des messages pour l'iteration en cours **************** * 'SI' IKLFFF ; 'MESSAGE' ; 'SINON'; 'MESSAGE' 'Iter'*13 'Nplas'*26 'Fresidu'*39 'Deps.max'*52 'Eps.max'*65 'Mresidu'*78 ; 'FINS'; INSTAB = FAUX; RESIDIN = RESIDU * 1.; hpp_exit = faux; 'SI' hpp_eps; hpp_exit=vrai; 'FINSI'; *======================================================================= *======= DEBUT DE LA BOUCLE DE CONVERGENCE ======= *======================================================================= * trac cach v1 nclk; * trac (surf1 et surf2 et surf3 et _rext et _rint) nclk; * trac cach (engr1 et engr2) nclk; dpsmaxp = 0; IT_RECA = 0; depst = depst * 0; depstp = depst; 'REPETER' ETIQ ; resmul = 2. * resmul; 'SI' (resmul > 1.0); resmul = 1.0; 'SINON'; itacc = 4; 'FINSI'; **resmul = 1.; * sauvegarde des etats sur la configuration debut de sous-increment: geom1 * mise a jour apres la sortie de etiq zdeptp = zdept; depstp = depst; zsig0 = zsig0_1; INSTAB = FAUX; nconvr = faux; *IT est le compteur de ETIQ, ITACC doit etre =< 0 pour qu'on accelere IT= IT + 1 ; ITACC = ITACC - 1; ZITAC = ZITAC + 1 ; * *--------------------------------------------------------------------- * 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 *----------------------------------------------------------------------- * * petits travaux pour acceleration de convergence * CORRECP = CORREC; CORREC = 0; PASTEST=FAUX; ACFEP0 = ACFP0; ACFEP0 = ACFEP0 - CORRECP ; 'SI' IGRD ; 'FORM' GEOM1; * rem : 0.1 est la valeur par defaut de EKREAC (= . 'REAC_GRANDS') 'SI' ((IT > 1) 'ET' (URG 'OU' (ITACC > 3) 'OU' INSTAB 'OU' URG = FAUX; IT_RECA = IT; dekreac1 = 0; INSTAB = FAUX; URG=FAUX; PASREINI=FAUX; * A cause de FORM, ZMATTEMP n'est plus parallele... zdepr = zdept; c_zdepr = faux; RECA_K = VRAI; RECA_N = RECA_N + 1; 'SI' (RECA_N > 20) ; nonconv = vrai; 'FINS'; 'SOUC' 0; 'SI' ((LAG_TOT 'EGA' 1) 'ET' vrai); 'FORM' GEOREF0; 'FORM' GEOR; 'DETRUI' HOOKRH; 'DETRUI' HOOKRH2; 'FORM' GEOM1; 'SINON'; 'FINSI'; PASOK = 'SOUCI'; 'SI' PASOK; 'FORM' GEOM1; 'FINSI'; ZMATTEMP = 0 ; * RH peut contenir des CL ZCLIM = ZCLIM0 'ET' ZCL ; 'SINON' ; ZCLIM = ZCLIM0 ; 'FINSI' ; BFCONT = FAUX ; BCLIM2 = FAUX ; 'SI' WTAB.'CONTACT'; * recalcul raideur de contact et jeu * Maintenant, on veut la configuration finale pour les directions et les jeux 'FORM' GEOM1; 'FORM' ZDEPT ; MODCON = wtab.'MODCONTA'; 'SINON' ; 'FINS' ; 'SI' ('NEG' 0 CJEU) ; 'FINSI' ; 'SI' (WTAB . 'MODAL') ; PBCDA = MCDAP 'POINT' &BCDA ; PCRR ='POINT' MCRR &BCDA ; 'NATURE' 'DISCRETE' ; 'SI' ('EGA' 1 &BCDA) ; CCDA = CHCR ; 'SINON' ; CCDA = CCDA 'ET' CHCR ; 'FINS' ; 'FIN' BCDA ; CDAP = CCDA ; 'FINS' ; * ATTENTION : mettre les conditions de frottement en premier pour * les numeroter en dernier 'SI' ('NEG' CRR 0); BCLIM2 = VRAI ; BFCONT = VRAI ; ZCLIM2 = CRR 'ET' ZCLIM2 ; CCOR = CRR '*' CDEPSLX ; ZFCONT = ZFCONT '+' CCOR ; 'FINS'; 'SI' ( 'NEG' 0 CJEU); BFCONT = VRAI ; ZFCONT = ZFCONT '+' CDAP ; 'FINS'; 'SI' ('NEG' 0 RFROT); BCLIM2 = VRAI ; BFCONT = VRAI ; ZCLIM2 = RFROT 'ET' ZCLIM2 ; CCOR = RFROT '*' CDEPSLX ; ZFCONT = ZFCONT '+' CCOR ; 'FINS'; * Mise a jour de ZCLIM, ZFCONSTA et RESIDU si necessaire 'SI' BCLIM2; ZCLIM = ZCLIM2 'ET' ZCLIM ; * pas de correction a faire sur residu car il est hors reactions pour qu'elles * sortent en total a la resolution incrementale 'FINS'; 'FINS'; 'SI' BFCONT; ZFCONSTA = ZFEXT2 '+' ZFCONT; 'SINO'; ZFCONSTA = ZFEXT2 ; 'FINS'; ZFCONST1 = ZFCONSTA; zclimp = zclim; 'SI' IDYN; ZFCONSTA = ZFCONSTA '+' ZFP1 ; 'FINS'; ZRAID = ZCLIM 'ET' ZRI; * repasser sur la configuration de calcul des raideurs 'FORM' geor; 'SI' iksia ; ksigtc= 'KSIGM' sigttcca zmodl zmat; ZRAID= ZRAID 'ET' ksigtc; 'FINS'; 'SI' IRCON; ZRAID = ZRAID 'ET' RIG_CONS; 'FINS'; * remise a jour impo12 IMPO12= FAUX; IMPO12U = FAUX; IMPO12U = VRAI; 'FINS'; IMPO12 = VRAI; DIMPOV = DFEXT0L '-' DIMPO12; 'FINS'; * pv 'SI' IRAUG; * xksx = xtmx zdep1d ksigtc ; *'SI' ((abs xksx) < xpetit ); * zdep1d = (zu1 + zdept) 'ENLE' 'LX'; * xksx = xtmx zdep1d ksigtc ; *'FINSI'; 'SI' (('VERI' xkx) 'ET' ('VERI' xmx) 'NON'); 'FINSI'; 'SI' faux; xksx1 = xkx - xktx; * 'SI' (xksx1 > xksx); xksx = xksx1 ; 'FINSI'; xksx = xksx1; 'SI' (('VERI' xksx) 'ET' ('VERI' xkx) 'ET' ('VERI' xmx) 'NON'); xksx = xkx; ** xksx = xty zu1l (ksigtc * zu1l) mnprim mndual; 'FINSI'; * si les coef sont negatif, c'est qu'on n'est pas instable sur ce mode 'SI' ((abs xksx) > xpetit ); augauto = xksx / xmx ; augk= xksx /xkx ; 'SINON'; 'FINSI'; 'FINSI'; *'SI' iprem; augauto = xkx / xmx; augk = 1.; *'FINSI'; augauto = abs augauto; augk = abs augk; augmult = augmult * 1.02; 'SI' (augmult > 1D4 ); augmult=1D4 ; 'FINSI'; augauto = augauto * augmult ; augk = augk * augmult; 'SI' (iraug 'ET' autaug) ; 'MESS' 'multiplicateur d augmentation masse' augauto ' raideur' augk ; 'FINSI'; ZRAIDNA= ZRAID ; 'SI' AUTAUG; ZRAID = ZRAID 'ET' (RIG_AUG * augauto) 'ET' (RH * augk) ; 'SINON'; ZRAID = ZRAID 'ET' RIG_AUG ; *** 'MESS' ' augmentation residuelle '; 'FINSI'; *** ZNACCE = FAUX; 'SINON'; augmult = augmult * 0.55 ; 'SI' IRAUGLU; zraid = zraid 'ET' rig_aug; 'FINSI'; 'FINSI'; *** 'FINS'; 'SI' WTAB.'PROCEDURE_CHARMECA'; PRECED . 'ADDI_MATRICE' = vrai; TFP22= CHARMECA PRECED WTAB.'T_FINAL'; PRECED . 'ADDI_MATRICE' = faux; zraid = zraid 'ET' TFP22.'ADDI_MATRICE'; 'FINS'; 'FINS'; 'SI' IDYN; * bp : en toute rigueur, il faudrait aussi recalculer la MASSE ... * et ajouter l'amortissement le cas échéant ... MMA = WTAB.'MASSE' '*' ( 4.D0 '/' WTAB.'DT' '/' WTAB.'DT'); ZRAID= ZRAID 'ET' MMA; 'FINS'; 'FORM' GEOM1; * on impose le recalcul de K a la prochaine iteration si it=2 'SI' (IT 'EGA' 2) ; ITACC = 5 ; 'SINON'; * pas de raison d'interferer avec les accelerations de convergence ***** ITACC = 1 ; 'FINS' ; GR_U_K = GR_U_FIN ; zdeptm = zdept ; 'FINS'; 'FINS'; * * acceleration de convergence effective ------------------------------ * ACCEL = FAUX; 'SI' (('NON' instab) 'ET' (ITACC '<EG' 0) 'ET' (IT '>' 3) 'ET' ZNACCE ); ITACC = 2; ACCEL = VRAI; CORREC = 'ACT3' ACFEP2 ACFEP1 ACFEP0 ACFP3 ACFP2 ACFP1 ACFP0 ; 'SI' (wtab.'STABILITE' 'ET' ('NON' HPP_EPS)); * verif que l'acceleration ne renvoie pas en arriere acc_ref = acc_ref + xpetit; acc_rap = acc_dir/acc_ref; 'SI' (acc_rap '>' 256.) ; * en cas d'acceleration trop grande, on la limite 'FINSI'; * 'SI' ((acc_rap '<EG' 0.) 'OU' (acc_dir '<EG' 0.)) ; 'SI' (acc_rap '<EG' 0.) ; * en cas d'acceleration retrograde, on n'accelere pas 'MESS' 'Annulation acceleration: instable ' ' ' acc_dir; correc = 0.; itacc=2; 'FINSI'; 'FINSI'; RESIDU = RESIDU '-' CORREC; RESIDC = RESIDC '-' CORREC; 'FINS'; ACFP3 = ACFP2 ; ACFP2 = ACFP1 ; ACFP1 = ACFP0 ; ACFEP2 = ACFEP1 ; ACFEP1 = ACFEP0 ; * * Resolution --------------------------------------------------------- * on obtient ZDEP1 = [ du ; Dlx ] 'SI' (resmul < 0.99); **ZDETOT = ZZD '+' ZDEPT ; **XXX1 = ZCLIM '*' ZDETOT; * -1*reactions a l'iteration it **FCORF = 'ENLEVER' XXX1 'FLX'; residu = ((residu -fcorf) * resmul) + fcorf; ** residu = (residu 'ENLE' 'FLX') 'ET' ((RESIDU 'EXCO' 'FLX' 'NOID' 'FLX') * resmul); 'FINSI'; 'SOUCI' 0; 'SI' IAFAIR; * 'SI' (WTAB.'ADHERENCE') ; 'SINON' ; FADRE = FADHE ; 'FINSI' ; FEXCI = FEXCI 'ET' FADRE ; 'FINSI' ; * 'SI' (WTAB.'CAFROTTE') ; FEXCI = FEXCI 'ET' FFROT ; 'FINSI' ; * 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12U); ZDEP1 BID BID BID BID = 'RESO' ZRAID 'SOUC' RESIDU 'INIB' STAB12.'RIBLO_M' STAB12.'LISEA_M' FEXCI ; 'SINON'; ZDEP1 BID BID BID = 'RESO' ZRAID 'SOUC' RESIDU 'INIB' STAB12.'RIBLO_M' STAB12.'LISEA_M' ; 'FINS'; 'SINON'; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12U); 'SINON'; 'SI' (WTAB.'MAN' 'ET' IPREM); ORDRE = WTAB.'ORDRE' ; ZDEP2 IOUT = CORMAN ZRAIDINI ZMODL ZMAT ORDRE ZU1 ZSIG0 RESIDNOR WTAB ; 'SI' (IOUT 'EGA' 1) ; ZDEP1=ZDEP2; 'FINS'; 'FINS'; 'FINS'; 'FINS'; monsouc ='SOUCI'; 'SI' monsouc; 'SI' ('NON' autaug); 'FINSI'; 'FINSI'; 'SI' LOGPIL ; * Calcul de D_eta par appel a la procedure PILOINDI D_ETA = PILOINDI PRECED ZU1 ZDEPT ZDEP1 ZDEPII DTAU ; * Mise à jour ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ; ETA = ETA '+' D_ETA ; 'FINSI' ; tabconv . it = 1; * * verif sens ---------------------------------------------------------- * monsouc = souci; pasunil = monsouc; 'SI' (WTAB.'STABILITE' 'ET' (('NON' HPP_EPS) 'OU' monsouc) 'ET' autaug ); itacc = 3 ; 'SI' monsouc ; resmul = resmul * 0.25; 'SI' autaug; iraug = vrai; 'FINSI'; HPP_EPS = FAUX; pastest = vrai; instab = vrai; urg = vrai; augmult = augmult * 1.4; 'SI' iraug; ZDEPT = zdeptp ; 'ITERER' etiq; 'FINSI'; 'FINS'; 'FINS'; * 'SI' IDYN ; STAB12.'ZRAIDV'= ZRAID; 'FINS'; STAB12.'RIBLO_M' = ZRAID_T.7; STAB12.'LISEA_M' = ZRAID_T.6; 'FINS'; 'FINS'; RECA_K = FAUX; * * 1ere iteration ------------------------------------------------------ 'SI' IPREM ; ZDEPT = 'COPIER' ZDEP1 ; ZDEPTM = ZDEPT ; * ZDEPf = ZDEPT ; ZDELA = 'COPIER' ZDEPT ; 'SI' (WTAB . 'MODAL') ; 'FINS'; 'FINS' ; * iprem est faux: on est apres la premiere operation--------------- * 'SINON'; * zdept est l'increment de deplacement total avec les lagrangiens * de la solution complete XXX1 = ZDEPT 'ENLEVER' 'LX' ; * ZDEPf = (ZDEPT + ZDEPTP) * 0.5; * on cumule les deplacements mais pas les lx * Du^(i) = Du^(i-1) + du ; Dlx^(i) = 0 + Dlx ZDEPT = XXX1 '+' ZDEP1 ; 'DETR' XXX1 ; 'FINS' ; * 'SI' (WTAB . 'MODAL') ; * mettre les point materiels dans zdept 'FINS'; 'FINS' ; * * ATTENTION ON EST SUR LA CONFIGURATION GEOM1? * * *--- CAS 2 ------------------------------------------------------------* * si on part dans les decors on redemarre a 0 'SI' ((XCONV > 1E8) 'ET' PASREINI 'ET' ('NON' IPILOT)) ; 'MESS' 'Reinitialisation du schema'; ZDEPT = zdeptp ; PASREINI=FAUX; ITACC=3; 'ITERER' ETIQ; 'FINS'; * *--- CAS 3 ------------------------------------------------------------* 'SI' ACCEL ; * pilotage automatique 'SI' IPILOT; RED1 = STAB12.'AUTOREDU' ; * reduction du critere de pilotage red2 sert d'incateur si * pas accelere tous les 2 pas * MODIFICATION CB215821 : 18/06/2015 'SI' (RED2 < (IT '/' 20)); STAB12.'AUTOREDU' = STAB12.'AUTOREDU' '*' 3.D0 ; RED2 = RED2 '+' 1 ; ITACC= 3; 'FINS'; RED1 = RED1 '/' (STAB12.'AUTOREDU') ; 'SI' (RED1 '<' 0.9) ; STAB12.AUTORED1 = 10 ; 'FINS' ; 'SI' ( STAB12.'AUTOREDU' '>' 1.D0 ) ; 'MESS' 'On divise le critere de pilotage par ' STAB12.'AUTOREDU'; 'FINS' ; * test si snap back et si refus de l'acceleration 'SI' IPLAVI ; XXX1 = 2.D0 '*' XXX2; XXX3 = ZDEPT '-' XXX1; 'SINON' ; XXX1 = COEPI '*' ZDELA; XXX3 = ZDEPT '-' XXX1 ; XXX1 = XXX2 '*' 2.D0 ; XXX4 = COEPI '*' ZDELA; XXX3 = ZDEPT '-' XXX4 '-' XXX1; 'FINS' ; XXX2 = XXX1 '*' 2.D0 ; XXX3 = ZDEPT '-' XXX2; ISNPB = FAUX; 'SI' (SRT '<' 0) ; ISNPB = VRAI ; 'FINS'; OO = WTAB.AUTOCRIT ; OO = OO '/' STAB12.'AUTOREDU' ; AL1 = OO '/' U1MA ; * al1 est le coefficient de normalisation 'SI' (( al1 '>EG' 1.d0) 'ET' (coepi '>EG' 1.d0)) ; * pour eviter d'aller au dela de alpha=1 on ignore le critere AL1 = 1.d0 ; 'FINS' ; 'SI'((AL1 '>' 1.D0) 'ET' (COEPI '>' 0.D0)); 'SI' ( AL1 '>' (1.D0 '/' COEPI)); AL1 = 1.D0 '/' COEPI; 'FINS'; 'FINS'; * normalisation PASTEST=VRAI; XXX1 = ZLX '*' ( 1.d0 '-' AL1); XXX3 = ZDEPT '*' AL1 ; ZDEPT = XXX3 '+' XXX1; 'DETR' XXX3; 'FINS'; 'SINON'; 'SI' IPILOT ; ** * le calcul d'al1 est fait pour eviter la convergence * s'il est externe a l'interval 0.5 --- 2. OO = WTAB.'AUTOCRIT' ; OO = OO '/' STAB12.'AUTOREDU' ; AL1 = OO '/' U1MA ; 'SI' (( AL1 '>EG' 1.d0) 'ET' (COEPI '>EG' 1.d0)) ; AL1 = 1.d0 ; 'FINS' ; 'FINS' ; 'FINS'; * * garder les reactions pour le test de convergence * ZDEPLP = ZDEPL ; *---------------------------------------------------------------------- * le nouveau champ est fixe on va tester l'equilibre(convergence) * et calculer la force motrice pour l'iteration suivante *---------------------------------------------------------------------- * * calcul de fcorf = lambda * m force de reaction ----------------- * fcoru = m * u depimp (flx) 'DETRUIRE' FCORF; * calcul du deplacement total u zdetotp = zdetot; ZDETOT = ZZD '+' ZDEPT ; XXX1 = ZCLIM '*' ZDETOT; * -1*reactions a l'iteration it FCORF = 'ENLEVER' XXX1 'FLX'; * valeur des relations imposees a l'iteration it 'DETR' XXX1 ; * * dans le cas des modeles endommageables de Lemaitre, on ecoule * en tenant compte, dans les iterations internes, de la variation du * materiau avec la temperature * ZMATT = ZMAT ; 'SI' ('ET' ('ET' ITHER ('OU' IENDOM IVIDOM ) ) WTAB.'MATVAR') ; * on recupere certain materiau avec les parametres fct de la temperatue * voir PAS_mate (il ne faut que la dependance thermique) 'FINS'; 'SI' IFEFP; * Update or total lagrangian ----------------------------------------- 'SI' IFEFPUL; * mess ' update lagrangian ZRIKTA'; ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP' ZMODL ZEPS0 ZVAR0 ZDEPT ZMAT ZPREK NITMA XUPDA; 'SINON'; * mess ' total lagrangian ZRIKTA'; chp_z = ZDEPT + ZU1 ; ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP' ZMODL ZEPS0 ZVAR0 chp_z ZMAT ZPREK NITMA ; chp_z = 1 ; 'FINS'; ZRAID = ZRIKTA 'ET' ZCLIM ; 'FORM' GEOM1 ; ** ACC = 'ABS' ( XXX1 - ACC0 ) ; *'DETR' XXX1 ; MMC = 'MASQUE' ACC 'SUPERIEUR' 1.D-10 'SOMME' ; 'SI' (MMC '>' MMCMAX) ; MMCMAX = MMC ; 'FINS' ; 'SINON'; * cas standard -------------------------------------------------------- ZMAT05=ZMAT; GEOM2 = GEOM1; 'SI' IGRD; zdepm = zdept '*' 0.5; 'FORM' GEOM1 ; 'FORM' GEOM1 ; 'FINS'; * * calcul de l'increment de deformation elast totale DEPST-------- * * Pour les hypotheses de deformations non lineaires, * on calcule l'increment de deformation sur le pas avec EPSI 'LINE' * mais en se placant sur une configuration intermediaire (mi pas) 'SI' (('NEG' HYPDEF 'LINEAIRE') 'ET' IGRD); 'SOUC' 0; 'SI' ITCAR; * A cause de FORM, ZMATTEMP n'est plus parallele... ZMATTEMP = 0 ; 'SINON'; 'FINS'; PASOK = 'SOUCI'; 'SI' ('NON' PASOK); 'FORM' geom1; PASOK = VRAI; 'FINS'; 'FINS'; 'SI' PASOK; resmul = resmul * 0.25; 'SI' autaug; iraug = vrai; 'FINSI'; HPP_EPS = FAUX; pastest = vrai; urg = vrai; itacc = 3; augmult = augmult * 1.4; ZDEPT = zdeptp ; 'SI' IGRD; 'FORM' GEOM1 ; 'FINS'; 'ITERER' etiq; 'FINSI'; 'SINON'; PASOK = vrai; 'FINS'; ** on n'a pas pu calcule DEPST : on utilise les deformations lineaires 'SI' PASOK; * 'MESS' 'deformation lineaire'; 'SI' IGRD; 'FORM' GEOM1; 'FINSI'; 'FINS'; * calcul des increments de deformation et contrainte sur la configuration initiale * en lagrangien total 'SI' IGRD ; 'SI' (LAG_TOT 'EGA' 1); * passage de depst et sig0 sur for0 'FORM' GEOREF0; 'FINSI'; * passage de depst et sig0 sur geom2 'SI' (LAG_TOT 'EGA' 2); 'FORM' geom2; 'FINS'; 'SI' (LAG_TOT 'EGA' 3); 'FORM' geom2m; 'FINS'; 'FINS'; * ICI depst et zsig0 sont sur la configuration initiale en lagrangien total et * il faut soustraire les parties thermique et deformations imposees 'SI' IPLAVI; 'SI' ('OU' ('OU' IENDOM IVIDOM) ICERAM); 'SI' IGRD; 'SINON'; 'FINS'; 'FINS'; * si thermoplastique on enleve e alpha *dt et d'autres termes * si le materiau depend de la temperature 'SI' (ITHER 'OU' WTAB.'MATVAR') ; XXX1 = DEPST0 '*' COEPI; XXX2 = DEPST '-' XXX1 ; * 'DETR' DEPST ; 'DETR' XXX1; DEPST = XXX2 ; 'FINS' ; 'SI' LOGDEF ; XXX1 = DDEFOR0 '*' COEPI; XXX2 = DEPST '-' XXX1 ; * 'DETR' DEPST ; 'DETR' XXX1; DEPST = XXX2 ; 'FINS' ; 'SINON'; * en elasticite, on a directement l'increment de contrainte 'SI' (('NON' ITHER) 'ET' WTAB.'MATVAR') ; XXDEFO = XXXXX4 '-' XXXXX3 ; DSIGT = DSIGT '+' DSIGTV ; * 'DETR' XXXXX3 ; 'DETR' XXXXX4 ; * 'DETR' XXDEFO ; 'DETR' DSIGTV ; 'FINS' ; * si thermoplastique on enleve alpha *dT et d'autres termes * si le materiau depend de la temperature 'SI' ITHER ; XXX1 = DSIGT0 '*' COEPI; XXX2 = DSIGT '-' XXX1; DSIGT = XXX2 ; 'FINS'; 'SI' LOGDEF ; XXX1 = DSI1 '*' COEPI; XXX2 = DSIGT '-' XXX1; * 'DETR' DSIGT ; 'DETR' XXX1; DSIGT = XXX2 ; 'FINS'; 'FINS'; * * Calcul du gradient des deplacements (total) a la fin du pas de ------ * temps par rapport a la configuration de reference GEOREF0 (initiale) * rappel : On a : ZDETOT = ZZD + ZDEPT ; 'SI' igrd; 'SI' ('NEG' GR_U_DEB 'INCONNU'); D_GR_U = GR_U_FIN '-' GR_U_DEB ; 'SINON'; 'FINS'; 'FINS'; * 'SI' ZCCONV ; STAB12.'DT' = DTINI '*' COEPI ; 'SINON' ; STAB12.'DT' = 0. ; 'FINS' ; * * si on est plastique ou viscoplastique on ecoule pour avoir----------- * le nouveau champ de contraintes * MMC = 0 ; 'SI' IPLAVI ; *...cas SSTE... 'SI' ISSTE; ZRIKTA ZSIGF ZVARF ZDEIF = 'SSTE' ZMODL ZSIG0 ZVAR0 DEPST ZMATT ZPREK NSSTE NITMA; *...cas non SSTE... 'SINON'; ZSIG01=ZSIG0; ZVAR01=ZVAR0; ZEPS01=ZEPS0; ** defineto = DEPIN0 ; TABCONT='TABLE'; 'SI' (nsoincr 'NEG' 1); ZSIGM=ZSIG0 '/' 2; 'FINS'; 'REPETER' sousinc nsoincr; tsodeb='FLOT' stab12 .'DT'/nsoincr*( &sousinc - 1)+conti.'TEMPS'; tsofin='FLOT' stab12 .'DT'/nsoincr* &sousinc +conti.'TEMPS'; 'SI' ('NON' ZCCONV) ; tsodeb = tsofin ; 'FINS' ; che1 = che1 che2 = che2 'FINS'; 'SI' ither ; TETDE=STAB12.'TET1' + (DTETD *('FLOT' (&sousinc - 1)/nsoincr)); TETDF=STAB12.'TET1' + (DTETD *('FLOT' &sousinc /nsoincr)); 'FINS' ; che1 = che1 'ET' che3 ; che2 = che2 'ET' che4 ; aa5 = DEPST '*' ( (&sousinc '-' 1.) '/' nsoincr); che5 = DEFT0 '+' aa5; cc5 = DEPST '*' ('FLOT' &sousinc '/' nsoincr); che6 = DEFT0 '+' CC5; *On met les deformations en tete des champs pour COMP ('KTAN' 'PERT') che1 = che5 'ET' che1 ; che2 = che6 'ET' che2 ; * pour les materiaux on garde toujours la valeur a la fin du pas * et au debut du pas * 'SI' ('NEG' ZMAT1 'INCONNU') ; che1 = che1 'ET' ZMAT1 ; * 'FINS' ; * * Modele NON LINEAIRE UTILISATEUR + GRANDES DEFORMATIONS * On ajoute les gradients de deplacements qui seront transformes en * gradients de transformation avant appel a UMAT (dans WKUMA1.ESO) 'SI' (( 'NEG' GR_U_DEB 'INCONNU') 'ET' igrd) ; gru1 = GR_U_DEB + ( D_GR_U * ('FLOT' (&sousinc-1) / nsoincr)); che1 = che1 'ET' gru1 ; gru2 = GR_U_DEB + ( D_GR_U * ('FLOT' &sousinc / nsoincr)); che2 = che2 'ET' gru2 ; 'FINS'; 'SI' LNLOC ; CHE1A = che1 ; CHE2A = che2 ; ISTEP = 1 ; che11 = CHE1A 'ET' chm_z ; che11 = che11 'ET' ZSIG01 'ET' ZVAR01 'ET' ZEPS01 ; che22 = CHE2A 'ET' chm_z ; che22 = che22 'ET' ZMATT ; 'SI' PARTLOCA ; cha11 = 0 ; cha22 = 0 ; cho22 = 0 ; 'SINON'; 'FINS' ; 'TYPE' 'VARIABLES INTERNES' ; * NLOC ne traitant pas des champs //, on reduit tout au modele initial. 'SI' (PARALLEL 'ET' ('NON' PARTLOCA)) ; 'FINS' ; * cas non-local MOYE ou SB 'SI' ('NEG' WTAB.'NON_LOCAL' 'HELM'); 'SI' ('EGA' WTAB.'NON_LOCAL' 'SB') ; MOD_SB = WTAB.'NLOC_MODL' ; ZVARF = ZVARF '+' CONTP '+' WTAB.'NLOC_SB_REGU' ; 'FINS' ; 'TYPE' 'VARIABLES INTERNES' ; 'SINON' ; * cas non-local HELM MOD_HELM = WTAB.'NLOC_MODL' ; 'REPE' BH NHELM ; LEMO = TAHELM . &BH. 'NOM' ; ZVNEW = 'CHAN' 'CHAM' ZVNEW MOD_HELM 'STRESSES' 'VARIABLES INTERNES'; DZVN = ZVNEW '-' ZVAUX ; ZVARN = ZVARN '+' DZVN ; * cas particulier SMAX 'TYPE' 'VARIABLES INTERNES' 'STRESSES' ; ZVARN = ZVARN '+' ZSMAN '-' ZSMAX ; 'FINS' ; 'FIN' BH ; 'FINS' ; * On rend ZVARN // si besoin : 'SI' (PARALLEL 'ET' ('NON' PARTLOCA)) ; 'FINS' ; ISTEP = 2 ; che11 = CHE1A 'ET' chm_z ; che11 = che11 'ET' ZSIG01 'ET' ZVARN 'ET' ZEPS01 ; che22 = CHE2A 'ET' chm_z ; che22 = che22 'ET' ZMATT ; 'SI' PARTLOCA ; cha11 = 0 ; cha22 = 0 ; cho22 = 0 ; 'SINON'; 'FINS' ; * Un peu de menage CHE1A = 1 ; CHE2A = 1 ; ZVARN = 1 ; ZVARF = 1 ; chm_z = 1 ; 'SINON' ; che11 = che1 'ET' ZSIG01 'ET' ZVAR01 'ET' ZEPS01 ; che22 = che2 'ET' ZMATT ; * si ( 'EXIST' PRECED 'ECRIT' ) ; * 'FINS'; 'SI' PARTLOCA; cha11 = 0 ; cha22 = 0 ; cho22 = 0 ; 'SINON'; * 'SI' ( 'EXIS' PRECED 'ECRIT' ); * 'FINS'; 'FINS'; 'FINS' ; 'FINS' ; 'SI' (nsoincr 'NEG' 1); 'SI' (&sousinc 'EGA' nsoincr); aaa1=zsigf/2; aaa2=aaa1+ zsigm; zsigm=aaa2; 'SINON'; aaa2=zsigf+ zsigm; 'DETR' zsigm; zsigm=aaa2; 'FINS'; 'FINS'; zsig01 = zsigf ; ZVAR01 = ZVARF ; ZEPS01 = ZDEIF ; ** defineTO=ZDEIF+defineTO; 'FIN' sousinc; * * Matrice tangente par perturbation evaluee pour la derniere iteration calculee * A voir : cas grand deplacement ZSIGF et ZMAT, cas poreux et thermique ZSIGF 'SI' (IKTAN 'ET' IPERT) ; Z1COMP = che11 ; Z2COMP = che22 'ET' ZSIGF ; 'FINS' ; * 'SI' (nsoincr 'NEG' 1); * pour tenir compte de ce que le travail de la correction est 1/2 FU, * on la multiplie par 2 ZSIGM = ZSIGM*(2. /nsoincr) ; ** ZDEPSPL=defineTO - ZEPS0 ; 'FINS'; * on enleve toutes traces de champs inutiles zsig01=1;ZVAR01=1;ZEPS01=1; ** defineTO=1; tabcont =1; * le fin d'en dessous est le fin de "si isste ... sinon ... finsi" * on est encore dans "si iplavi ... finsi" 'FINS'; *...fin du si ISSTE sinon ... * * cas particulier poreux et thermique * 'SI' ITHER; 'SI' POR1 ; ZSIGF = ZSIGF - ( COEPI * DMSRT0 ) ; 'FINS'; 'FINS'; * * max de epse pendant l'iteration * nbre de points qui ont une evolution non lineaire * 'SI' ('EGA' WTAB.'MOVA' 'RIEN') ; ** EPSM = 0.; DPSMAX=0.; MMC=0; 'SINON'; ** EPSM = 'MAXI' XXX1 ; ACC = 'ABS' ( XXX1 - ACC0 ) ; ** 'DETR' XXX1 ; MMC = 'MASQUE' ACC 'SUPERIEUR' 1.D-10 'SOMME' ; 'SI' (MMC > MMCMAX) ; MMCMAX = MMC ; 'FINS' ; ** DPSMAX = 'MAXI' ACC ; 'FINS'; * * dans le cas SANS plasticite,viscoplas,endommagement * on calcule directement le nouveau champ de contrainte * 'SINON' ; ZSIGF = ZSIG0 + DSIGT ; 'FINS'; * si il y a lieu transport sur configuration debut de pas 'SI' IGRD ; ** ZSIG0 = ZSIG0_1; 'SI' (LAG_TOT 'EGA' 1); ** ZSIG0 = 'PICA' HYPDEF ZSIG0 ZMODL ZU1; 'FORM' GEOM1; 'FINS'; 'SI' (LAG_TOT 'EGA' 2); 'FORM' GEOM1; ** ZSIG0 = 'CAPI' HYPDEF ZSIG0 ZMODL ZDEPT; 'FINS'; 'SI' (LAG_TOT 'EGA' 3); 'FORM' GEOM1; ** ZSIG0 = 'CAPI' HYPDEF ZSIG0 ZMODL ZDEPM; 'FINS'; 'FINS'; * DEPSTI = DEPST; 'SI' (IGRD 'ET' ('NON' HPP_EPS) 'ET' ('EGA' LAG_TOT 1)); 'FORM' georef0; 'FINSI'; EPSMX = DEFT0 '+' DEPSTI; 'DETR' EPSMX; 'SI' (IGRD 'ET' ('NON' HPP_EPS) 'ET' ('EGA' LAG_TOT 1)); 'FORM' geom1; 'DETR' DEPSTI; 'FINSI'; * *--- CAS 4 ------------------------------------------------------------- * on determine les forces resultantes de l'increment------------------ * * en grand deplacements on transforme pi en sigma * on met a jour les coordonnees * ZMAT2=ZMAT; 'SI' (IGRD 'ET' (('NON' HPP_EPS) 'OU' (IT > 0 ))); ref = 50; coefmul = 1.d0 / (DPSMAX + XPETIT); 'SI' (((coefmul < ref) 'ET' AUTAUG) 'OU' PASUNIL); * il faut reduire d'urgence les deplacements ; coefmul = 0.0; 'SI' ((red_urg > 1) 'ET' ('NON' instab)); resmul = resmul * 0.25; 'FINSI'; augmult = augmult * 1.5; urg = vrai; 'SI' autaug; iraug = vrai; 'FINSI'; * on garde zdep1d pour le calcul des coeff d'acceleration 'SI' AUTAUG; * reprise du vieux zdept. Tout le reste est recalcule a chaque pas de etiq ZDEPT = zdeptp ; itacc = 3; 'FINSI'; 'SI' ((dpsmaxp < 5e-4) 'OU' (dpsmaxp > 2e-2 )) ; PASTEST=VRAI; 'SINON'; 'SI' (red_urg > 4); nconvr = vrai; 'FINSI'; tabconv . it = 1; 'FINSI'; 'SI' (RED_URG > 2 ) ; ZMAXIT = IT+5; 'FINS'; nonconv = vrai; 'FINS'; 'SI' (resmul < 0.99); resmul = resmul * 0.25; 'FINSI'; RED_URG = RED_URG + 1; 'SI' ((dpsmaxp > 1.8d-2) 'ET' (augmult > 100)); zprecnc= 2.1d-2; 'SI' (zmaxit > it) ; zmaxit = it - 1; 'FINSI'; 'FINSI'; 'SI' ((zmaxit > it) 'ET' (RED_URG > 4 )); ZMAXIT = IT-1; 'FINS'; 'ITERER' etiq; 'FINS'; 'FINS'; 'FINS'; * au dessus le finsi du si fefp sinon(Update or total lagrangian) ---- * dpsmaxp = dpsmax; * fin de boucle de reduction de zdept *--------------------------------------------------------------------- 'SI' ('NON' IFEFP); 'SI' (IGRD 'ET' ('NON' HPP_EPS)); zdepf = zdept ; * grands deplacements ----------------------------------------------- 'SOUC' 0; 'SI' ITCAR ; 'SINON'; 'FINS'; PASOK = 'SOUCI'; 'SI' PASOK; dekreac1 = 1e30; pastest = vrai; instab = vrai; urg = vrai; 'SI' autaug; iraug = vrai; 'FINSI'; tabconv . it = 1; 'FORM' geom1; 'SI' AUTAUG; * reprise du vieux zdept ZDEPT = zdeptp ; ** dsigt = dsigt * 0.; 'ITERER' etiq; 'FINSI'; 'FINS'; 'SI' IRCON; 'FINS'; 'SI' LOGPRE ; ZDFINI = ZDFINI - ZFPEXTF ; 'DETR' ZFPEXTF ; 'SINON' ; 'FINS' ; ZDFINI = ZDFINI + ZFPEXTF ; 'FINS' ; 'SI' ADDISEC2; ZDFINI = ZDFINI '-' FP22 ; 'FINS'; 'SI' WTAB.'PROCEDURE_CHARMECA';