* UNPAS PROCEDUR JK148537 24/10/29 21:15:10 12056 *----------------------------------------------------------------------* * PROCEDURE UNPAS * * * * Calcul d'un increment de solution en grand deplacement plastique par * * la methode des residus. * * * * Les differentes configurations qui interviennent sont : * * WTAB.'FOR0' : configuration de debut de calcul * * WTAB.'GE0_DEB' : * * -> en GRANDS_DEPLACEMENTS ou FEFP_FORMULATION (1) : * * configuration de debut de pas * * -> dans les autres cas (2) : * * configuration de debut de calcul (= WTAB.'FOR0) * * Dans le cas (1), la configuration WTAB.'GE0_DEB' est actualisee * * dans PASAPAS avec GEOM2. * * Dans le cas (2), il n'y a pas d'actualisation de la configuration * * au cours du calcul. * * * * Pour les calculs DYNAMIQUE, un schema de Newmark (implicite) * * avec gamma = 1/2 et beta = 1/4 est utilise. * * * *----------------------------------------------------------------------* LOG_CNTRL = VRAI; *CB215821 : Recuperation de XPETIT (07/12/2016) * WTAB = PRECED.'WTABLE' ; conti = PRECED.'CONTINUATION'; estim = PRECED.'ESTIMATION' ; * * Liste de composantes utiles 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' 'EPS' ; * Pour MOD_LIAISON : MLPRIM_LIA = MLPRIM ; MVPRIM_LIA = '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'; MLDUAL_LIA = MLDUAL ; *=DEB===== Formulation HHO = Definition des DDLs =============================== * Definition des DDLs primaux : inconnues de deplacements des cellules et des Faces MNPRIM_HHO = 'MOTS' * 'UXC0' 'UXC1' 'UXC2' 'UXC3' 'UXC4' 'UXC5' 'UXC6' 'UXC7' 'UXC8' 'UXC9' * 'UYC0' 'UYC1' 'UYC2' 'UYC3' 'UYC4' 'UYC5' 'UYC6' 'UYC7' 'UYC8' 'UYC9' * 'UZC0' 'UZC1' 'UZC2' 'UZC3' 'UZC4' 'UZC5' 'UZC6' 'UZC7' 'UZC8' 'UZC9' 'UXF0' 'UXF1' 'UXF2' 'UXF3' 'UXF4' 'UXF5' 'UXF6' 'UXF7' 'UXF8' 'UXF9' 'UYF0' 'UYF1' 'UYF2' 'UYF3' 'UYF4' 'UYF5' 'UYF6' 'UYF7' 'UYF8' 'UYF9' 'UZF0' 'UZF1' 'UZF2' 'UZF3' 'UZF4' 'UZF5' 'UZF6' 'UZF7' 'UZF8' 'UZF9' ; MNDUAL_HHO = 'MOTS' * 'FXC0' 'FXC1' 'FXC2' 'FXC3' 'FXC4' 'FXC5' 'FXC6' 'FXC7' 'FXC8' 'FXC9' * 'FYC0' 'FYC1' 'FYC2' 'FYC3' 'FYC4' 'FYC5' 'FYC6' 'FYC7' 'FYC8' 'FYC9' * 'FZC0' 'FZC1' 'FZC2' 'FZC3' 'FZC4' 'FZC5' 'FZC6' 'FZC7' 'FZC8' 'FZC9' 'FXF0' 'FXF1' 'FXF2' 'FXF3' 'FXF4' 'FXF5' 'FXF6' 'FXF7' 'FXF8' 'FXF9' 'FYF0' 'FYF1' 'FYF2' 'FYF3' 'FYF4' 'FYF5' 'FYF6' 'FYF7' 'FYF8' 'FYF9' 'FZF0' 'FZF1' 'FZF2' 'FZF3' 'FZF4' 'FZF5' 'FZF6' 'FZF7' 'FZF8' 'FZF9' ; * Pour l'istant on ne considere que les inconnues (DEPL et FORC) de face pour * le calcul de l'equilibre et l'acceleration de convergence * D'ou les lignes commentees ci-dessus dans la definition de MNPRIM_HHO et MNDUAL_HHO MNPRIM = MNPRIM 'ET' MNPRIM_HHO ; MLPRIM = MLPRIM 'ET' MNPRIM_HHO ; MNDUAL = MNDUAL 'ET' MNDUAL_HHO ; MLDUAL = MLDUAL 'ET' MNDUAL_HHO ; MLDEPL = MLDEPL 'ET' MNPRIM_HHO ; *A Ameliorer : ne faire que si presence HHO et n'ajouter que les DDLs utilises par HHO ! *=FIN===== Formulation HHO ===================================================== * * Pour stocker des informations necessaires aux calculs d'usure *----------------------------------------------------------------------* * Options de pasapas - parametres du calcul * *----------------------------------------------------------------------* * 1- rigidite IKSIA = WTAB.'K_SIGMA'; IKTAN = WTAB.'K_TANGENT' ; ZKTASYM = 'TEXTE' ' ' ; 'SI' WTAB.'K_TANGENT_SYME' ; 'FINS'; IMPLP = WTAB.'LIAISON_PERSISTANTE'; AUTAUG = WTAB.'AUTOAUGM'; IRAUGLU = WTAB.'RAIDAUGM'; IRAUG = FAUX; 'SI' IRAUGLU ; RIG_AUG = WTAB.'RIGIDITE_AUGMENTEE'; 'SI' ('NON' AUTAUG) ; IRAUG = VRAI ; 'FINS'; 'FINS'; IRCON = WTAB.'RAIDCONST'; 'SI' IRCON; RIG_CONS = WTAB.'RIGIDITE_CONSTANTE'; 'FINS'; * declenchement du recalcul de la matrice DLTAIT = WTAB.'DELTAITER' ; ITRCLC = 1 '*' DLTAIT ; 'SI' (ITRCLC < 20) ; ITRCLC = 20; 'FINSI'; * * 2- Type de formulation IFEFP = WTAB.'FEFP_FORMULATION' ; LAG_TOT= WTAB.'LAG_TOT'; ISSTE = WTAB.'SUBSTEPPING'; IFEFPUL= WTAB.'UPDATE_LAGRANGIAN'; 'SI' (LNLOC 'ET' ('EGA' WTAB.'NON_LOCAL' 'HELM')); TAHELM = WTAB.'HELMHOLTZ' ; NHELM = TAHELM . 'N_VARI_NL' ; 'FINSI' ; * * 3- Type de materiau ICERAM = WTAB.'CERAMIQUE' ; IENDOM = WTAB.'ENDOMMAGEMENT'; IPLAVI = WTAB.'IPLAVI'; POR1 = WTAB.'POR1' ; IVIEXT = WTAB.'VISCO_EXTERNE'; IVIDOM = WTAB.'VISCODOMMAGE'; IVISCO = WTAB.'VISCOPLASTIQUE'; * * 4- Critere de convergence/non-convergence ZMAXIT = WTAB.'MAXITERATION' ; NITMA = WTAB.'NITERINTER_MAX'; ZPREK = WTAB.'PRECISINTER' ; ZPREC = WTAB.'PRECISION' ; ZPRECM = WTAB.'PRECFLEX' ; ZPRECHPP = ZPREC ; IFTOL = 'NEG' WTAB.'FTOL' 'INCONNU' ; 'SI' IFTOL ; ZFTOL = 'ABS' WTAB.'FTOL' ; 'FINS'; IMTOL = 'NEG' WTAB.'MTOL' 'INCONNU'; 'SI' IMTOL ; ZMTOL = 'ABS' WTAB.'MTOL' ; 'FINS'; 'SI' ('OU' IVISCO IVIDOM IVIEXT); ZPREK = 5.E-7 ; 'FINS'; 'SI' IENDOM; ZPREK = ZPREC ; 'FINS'; NSOINCRN = WTAB.'SOUS_INCREMENT' ; * * 5- Type de calcul IPILOT = WTAB.'AUTOMATIQUE'; ISOL = WTAB.'CONSOLIDATION'; IDYN = WTAB.'DYNAMIQUE'; IGRD = WTAB.'GRANDS_DEPLACEMENTS'; IJAUMA = 'EGA' HYPDEF 'JAUMANN' ; HPP_EPS = 'EGA' WTAB.'PREDICTEUR' 'HPP'; * * 6- Chargement particulier LOGDEF = WTAB.'CHAR_DEFI'; LOGPRE = WTAB.'CHAR_PRES' ; LOGPIL = WTAB.'CHAR_PILO' ; ITHER = WTAB.'ITHER' ; * * 7- Acceleration de convergence ZNACCE = VRAI; * * 8- Instant de debut et de fin de pas TEMPS0 = WTAB.'TEMPS0'; TI = WTAB.'T_FINAL'; 'SINON'; TEMPS0 = conti . 'TEMPS'; TI = estim . 'TEMPS'; 'FINSI'; * * 9- Blocages mecaniques (BLOM) 'SI' WTAB.'CHAR_BLOM' ; WTAB.'BLOCAGES_MECANIQUES' = BLOM1 ; 'FINSI' ; ZCLIM0 = WTAB.'BLOCAGES_MECANIQUES'; ZCLIM = ZCLIM0 ; *----------------------------------------------------------------------* * Verification de la compatibilite entre options * *----------------------------------------------------------------------* 'SI' (IKTAN 'ET' ('NON' IPLAVI)) ; ' on utilise la rigidite elastique' ; IKTAN = FAUX ; 'FINS' ; IPERT = WTAB.'K_TANGENT_PERT' 'ET' ('NON' LNLOC) 'ET' IPLAVI ; * Matrice tangente : pas d'acceleration en cas de modele FEFP ou SSTE 'SI' IKTAN ; 'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS'; 'FINS' ; * * Configurations de reference et de debut de pas GEOREF0 = WTAB.'FOR0' ; GEOM1 = WTAB.'GE0_DEB' ; *----------------------------------------------------------------------* * Modele complet * * * * Notations utilisees : * * ZMODLI : modeles mecanique + poreux (non //) * * -> aucune parallelisation alors ZMODL = ZMODLI * * -> parallelisation comportement alors ZMODL = ZMODLI * * MODRELOC est // * * -> parallelisation automatique alors ZMODL est // * *----------------------------------------------------------------------* * PAS_MODL : mise a jour des indices de PRECED.WTABLE relatifs aux * modeles si PRECED.WTABLE.MODELE a ete modifie. PAS_MODL PRECED ; 'SI' WTAB.'FOR_MECA' ; PARALLEL = FAUX ; PARTLOCA = FAUX ; ZMODLI = WTAB.'MO_TOT' ; ZMODL = ZMODLI ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'COMPORTEMENT') ; PARALLEL = VRAI ; PARTLOCA = VRAI ; 'FINS' ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'AUTOMATIQUE') ; PARALLEL = VRAI ; PARTLOCA = FAUX ; 'FINS' ; * * Si le modele a change, adapter egalement les MCHAMLS 'SI' ('NEG' ZMODLI WTAB.'MO_TOT_PREC') ; * 'SI' IPLAVI ; 'FINS'; * * 'SI' ITHER ; 'SI' POR1; 'FINS'; 'FINS'; * 'SI' LOGDEF ; 'FINS'; * 'FINS'; 'FINSI'; *----------------------------------------------------------------------* * Champ materiau et caracteristiques a l'instant TEMPS0 * * * * Notations utilisees : * * ZMAT11 : materiau a l'instant TEMPS0 (non //) * * -> aucune parallelisation alors ZMAT1 = ZMAT11 * * -> parallelisation comportement alors ZMAT1 = ZMAT11 * * -> parallelisation automatique alors ZMAT1 est // * *----------------------------------------------------------------------* ZMAT11 = WTAB.'MAT1' ; ZMAT1 = ZMAT11; 'SI' WTAB.'FOR_MECA' ; 'SI' ('EGA' WTAB.'PROCESSEURS' 'AUTOMATIQUE') ; 'FINS'; * * Distinguer les modeles poreux des modeles mecaniques 'SI' POR1; * ** kich ma_por0 intervient si ISOL ** initialement MA_POR0= 'REDU' MO_POR MAT1 ; 'FINS'; 'FINSI'; *----------------------------------------------------------------------* * Champ materiau et caracteristiques a l'instant TI * * * * Notations utilisees : * * ZMAT22 : materiau a l'instant TI (non //) * * -> aucune parallelisation alors ZMAT2 = ZMAT22 * * -> parallelisation comportement alors ZMAT2 = ZMAT22 * * -> parallelisation automatique alors ZMAT2 est // * *----------------------------------------------------------------------* * 'SI' WTAB.'FOR_MECA' ; 'SI' WTAB.'ITCAR'; 'FORM' GEOREF0 ; 'FORM' GEOM1 ; 'FINS'; 'SI' POR1; 'FINS'; 'FINSI' ; *----------------------------------------------------------------------* * Chargement a imposer a l'instant TI (F^ext_n+1) * *----------------------------------------------------------------------* 'SI' ('NEG' TYP_2 'CHPOINT '); 'FINS'; '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' ; * ZFEXT2 = ZFEXT2 + F2_DEP; 'FINS'; * zclim_mail = zclim extrai 'MAILLAGE'; ZFEXT2 = ZFEXT2 + F2_INC + F2_base; ; 'FINS'; * 'SI' LOGDEF ; 'FINS'; * 'SI' WTAB.'FOR_MECA' ; 'SI' ITHER ; 'SI' WTAB.'CHAR_THE'; 'FINSI'; 'SI' WTAB.'FOR_THER'; WTAB.'TET2' = estim.'TEMPERATURES' ; 'FINSI'; * Champs de temperature en debut (ZTEMP1) et fin de pas (ZTEMP2) * Si les maillages mecanique/thermique sont differents, * on projette le champ de temperature sur le modele mecanique 'SI' WTAB.'PROJECTION'; 'FORM' GEOREF0 ; 'SINON' ; 'FINSI' ; 'FORM' GEOM1 ; * Sinon, on vient seulement recuperer les champs de temperature 'SINON' ; ZTEMP1 = WTAB.'TET1'; ZTEMP2 = WTAB.'TET2'; 'FINSI' ; 'SI' POR1; 'FINS'; 'SINON' ; 'FINSI'; 'FINSI'; *----------------------------------------------------------------------* * Dynamique : ajout au second membre de * * F0 + 4/DT*M*V0 - B0*SIG0 * *----------------------------------------------------------------------* * Matrice de masse 'SI' (IDYN 'OU' WTAB.'FREQUENTIEL') ; ICALMAS = FAUX ; ICALMAS = VRAI ; 'FINSI' ; 'SI' (IGRD 'ET' IDYN) ; ICALMAS = VRAI ; 'FINSI' ; 'SI' ICALMAS ; ZMAT22 WTAB.'MO_TOT' ) ; 'SI' (IDYN 'ET' WTAB.'MASSCONST') ; WTAB.'MASSE' = WTAB.'MASSE' 'ET' WTAB.'MASSE_CONSTANTE'; 'FINS'; 'FINS'; 'FINS'; * 'SI' IDYN ; UNSURH = 1.D0 '/' WTAB.'DT' ; 'SI' ('EGA' WTAB.'FREA1' 'INCONNU'); * * Forces exterieures a TEMPS0 F1 = F1 '+' F1M ; 'FINS'; * F1 = F1 '+' F1F ; 'FINS'; * 'SI' WTAB.'PROCEDURE_CHARMECA'; TFF1 = CHARMECA PRECED TEMPS0 ; FF1 = TFF1 .'ADDI_SECOND'; 'FINS'; F1 = F1 '+' FF1; 'FINS'; * 'SI' LOGPRE ; MOP = WTAB.'MOD_PRE' ; 'SINON' ; 'FINS' ; F1 = F1 '+' FF1; 'FINS'; * * Forces interieures a TEMPS0 *HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO) 'SI' IRCON ; LAF0 = LAF0 'ET' 'FINS'; * * FREA1 : (masse*acceleration initiale)+(amortissement*vitesse initiale) WTAB.'FREA1' = F1 '-' LAF0 ; 'FINS'; * 'SI' IMPLP ; * 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'; * FF = WTAB.'MASSE' '*' conti.'VITESSES'; * * partie du second membre qui ne depend que des informations du pas prec WTAB.'FREA1' = FF4 '+' WTAB.'FREA1'; 'FINS'; * * Advection 'SI' WTAB.'NVSTNL' ; LOGADV = VRAI ; 'SINON' ; LOGADV = FAUX ; 'FINSI' ; 'FINSI' ; *----------------------------------------------------------------------* * Consolidation : ajout au second membre de * * -B0*SIG0 + DT*(1-TETA)*FI0 + DT*H*P * *----------------------------------------------------------------------* 'SI' ISOL ; FF4 = 'EXCO' WTAB.'MOT_POR' FF WTAB.'MOT_POR' 'NOID' 'NATURE' 'DISCRET' ; * * ---- traitement des flux si besoin ---- FACFLU = -1. * WTAB.'DT'; FLUXT = ( FACFLU * (1 - WTAB.'TETA') * FLUXT0 ) + ( FACFLU * WTAB.'TETA' *FLUXTI ) ; ZFEXT2 = ZFEXT2 '+' FLUXT ; 'FINS' ; 'FINS'; *----------------------------------------------------------------------* * Pilotage indirect : preparation * *----------------------------------------------------------------------* 'SI' LOGPIL ; ZFPILIN = WTAB.'FORCES_PILOTEES' 'ET' WTAB.'DEPLACEMENTS_PILOTES' ; ETA0 = 'EXTR' PRECED.'COEFFICIENT_DE_PILOTAGE' D_ETA=0.; 'FINS' ; *----------------------------------------------------------------------* * Partie constante du second membre sur le pas de temps * *----------------------------------------------------------------------* ZFCONST1 = ZFEXT2 ; 'SI' IDYN ; ZFEXT2 = ZFEXT2 '+' WTAB.'FREA1'; 'FINSI' ; *----------------------------------------------------------------------* * Recuperation de valeurs de WTAB pour initialisations eventuelles * *----------------------------------------------------------------------* 'SI' WTAB.'FOR_MECA' ; * ZXDENO = WTAB.'XDENO' ; ZXDENOM = WTAB.'XDENOM' ; * XDENOo = WTAB.'XDENO'; * XDENOMo = WTAB.'XDENOM'; ZETAT1 = WTAB.'ETAT1'; 'SI' LOGDEF ; ZDEFOR1 = WTAB.'DEFOR1' ; ZDEFOR2 = WTAB.'DEFOR2' ; 'FINS'; 'SI' ITHER ; ZTET1 = WTAB.'TET1' ; ZTET2 = WTAB.'TET2' ; 'FINS'; ZLASTKTAN = WTAB.'LASTKTAN' ; * 'SI' ('NEG' WTAB.'RIBLO_M' 'INCONNU') ; ZRIBLO_M = WTAB.'RIBLO_M' ; ZLISEA_M = WTAB.'LISEA_M' ; 'FINS' ; * ZFNONL = WTAB.'FNONL' ; 'SI' ('NEG' WTAB.'INCREMENT' 'INCONNU'); INCRPREC = WTAB.'INCREMENT'; 'FINS'; * * teste t'on les moments ? * teste t'on les POREUX ? 'SI' POR1 ; TSTMOM=VRAI ; 'FINS'; * IKLFFF = FAUX ; 'SI'TSTMOM; 'SI' IFTOL; 'SI' IMTOL; IKLFFF=VRAI; 'FINS'; 'FINS'; 'SINON' ; 'SI' IFTOL; IKLFFF=VRAI; 'FINS'; 'FINS'; 'SI' IKLFFF ; 'FINSI' ; * * Jacobien du modele pour ponderer les champs de act3 zjac = zjac **-0.5; * * Initialisation CHArgement SANS T : * * isoucomp : Indicateur SOUci dans COMPortement, FAUX par defaut isoucomp = faux ; PASUNIL = FAUX; ZNITE = 0 ; ISOUSPPP = WTAB.'ISOUSPAS' ; WTAB.'ISOUSPAS' = 0; 'FINSI' ; * RED_URG = 0; resmul = 1; augmult = 0.30000000; augauto = augmult; augk = augmult; XCONV = 0.; DEPSTDM = 0.; dpsmax = 0; KNOCONV = 0 ; *----------------------------------------------------------------------* * Valeurs des champs en debut de pas TEMPS0 * *----------------------------------------------------------------------* 'SI' WTAB.'FOR_MECA' ; ZDEP0 = CONTI.'DEPLACEMENTS' ; 'SI' IPLAVI ; 'SI' ('NEG' WTAB.'MOVA' 'RIEN') ; 'FINS' ; 'SI' ISOL ; 'FINS' ; 'FINS' ; 'SI' ('NEG' WTAB.'MOD_MEC' 'INCONNU') ; 'SI' (IGRD 'ET' WTAB.'UTILISATEUR'); 'FORM' GEOREF0 ; 'FORM' GEOM1 ; 'FINS' ; 'FINS' ; 'FINSI' ; *######################################################################* *----------------------------------------------------------------------* * Boucle de non convergence BONOCONV (Configuration GEOM1) * *----------------------------------------------------------------------* 'REPETER' BONOCONV WTAB.'MAXSOUSPAS' ; * trac cach v1 nclk; * * Initialisation objectif de fin de pas. Forces et deplacements ZFCONSTA = ZFEXT2 ; * 'SI' IGRD ; 'FORM' GEOM1; 'FINSI' ; * excfconv = 0; augmulto = augmult; 'SI' ('EGA' RED_URG 0); 'SI' ((XCONV < ZPREC) 'ET' (DEPSTDM < ZPREC) 'ET' (resmul > 0.99) 'ET' (AUGMULT < 100.)); augmult = 0.30000000; *pv resmul = 1; 'SI' autaug ; DE_CNTRL = FAUX; 'SI' (dpsmax < zprec); iraug = faux; 'FINSI'; 'FINSI'; znacce=VRAI; 'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS'; 'SINON'; augmult = augmult / 1.8; si (augmult < 1d-3); augmult = 1d-3; finsi; 'FINSI'; 'SINON'; augmult = augmult * 1.5; 'SI' autaug; IRAUG = VRAI; 'FINSI'; 'FINSI'; augk = augmult / augmulto * augk; augauto = augmult / augmulto * augauto; * 'SI' WTAB.'FOR_MECA' ; KNOCONV = KNOCONV+1 ; 'SI' (knoconv '>' 1) ; HPP_EPS = FAUX ; 'FINSI' ; * ZDEP0N = -1. '*' ZDEP0 ; ZDEPT = ZDEP0 '-' ZD0SLX ; ZSIGF = ZSIG0 ; * 'SI' ITHER ; DTETD = ZTEMP2 '-' ZTEMP1 ; 'FINS' ; * * -------------------------------------------------------------------- * Reevaluer les materiaux en convergence forcee (TEMPS0 =TI) 'SI' (knoconv '>' 1) ; 'SI' WTAB.'ITCAR'; 'FORM' GEOREF0 ; 'FORM' GEOM1 ; 'FINS'; * 'SI' (IGRD 'ET' ('EGA' LAG_TOT 1)) ; ZMAT2I = ZMAT22 ; 'FINSI' ; 'SI' WTAB.'ITCAR'; 'FORM' GEOREF0 ; 'FORM' GEOM1 ; 'FINS'; 'FINS' ; * * -------------------------------------------------------------------- * Rigidite a la fin du pas * 'SI' ('OU' IENDOM IVIDOM ICERAM); 'DETR' HOOKENDO; 'SINON'; 'SI' (IGRD 'ET' ('EGA' LAG_TOT 1)) ; 'FORM' GEOREF0; 'FORM' GEOM1; 'SINON'; 'FINSI'; 'FINS'; * * RH peut contenir des CL ZCLIM = ZCLIM0 'ET' ZCL ; 'FINSI' ; RRRR = RH 'ET' ZCLIM0 ; * * Option 'RIGIDITE_CONSTANTE' 'SI' IRCON; RRRR = RRRR 'ET' RIG_CONS; 'FINS'; * * Option 'RIGIDITE_AUGMENTEE' 'SI' IRAUGLU ; 'SI' ('NON' (AUTAUG 'ET' IRAUG)) ; RRRR = RRRR 'ET' RIG_AUG ; 'FINSI'; 'FINSI'; * * Option 'AUGMENTATION_AUTOMATIQUE' 'SI' (AUTAUG 'ET' ('NON' IRAUGLU)); RIG_AUG = 'MASSE' ZMODL ZMAT2 ; *** mess 'actualisation rig_aug'; 'FINSI'; 'SI' (AUTAUG 'ET' IRAUG); RRRR = RRRR 'ET' (RIG_AUG * augauto) 'ET' (RH * augk) ; ** 'MESS' 'multiplicateur d augmentation masse' augauto ' raideur' augk ; 'FINS'; * * Stockage de la rigidite pour eviter de la recalculer WTAB.'RRRR'=RRRR; 'FINSI'; * 'FINS'; 'DETR' ZRAID; ZRAID = WTAB.'RRRR' ; 'FINSI' ; *FOR_MECA * 'SI' WTAB.'NVSTNL' ; 'SI' ('NON' ('EXISTE' WTAB 'RNSL')) ; WTAB.'MAT_NSL' = WTAB.'MAT_NSL' 'ET' WTAB.'NSMA' = 'MASSE' WTAB.'MOD_NSL' WTAB.'MAT_NSL' ; WTAB.'RNSL' = RINSL 'ET' ZCLIM0 'ET' ndiv ; 'FINSI' ; 'FINSI' ; * *------------ consolidation ou dynamique faut-il recalculer l'operateur? 'SI' (ISOL 'OU' IDYN); 'SI' ( '>' (WTAB.'DTPREC' '*' 0.9999) WTAB.'DT') ; WTAB.'RECAOP' = VRAI; 'FINS'; 'SI' ( '<' (WTAB.'DTPREC' '*' 1.0001) WTAB.'DT') ; WTAB.'RECAOP' = VRAI; 'FINS'; 'SI' WTAB.'MATVAR'; WTAB.'RECAOP' = VRAI; 'FINS'; 'SI' ('NON' WTAB.'RECAOP') ; 'SI' ('NEG' WTAB.'OPERATEUR' 'INCONNU'); ZRAID = WTAB.'OPERATEUR'; 'FINS'; 'FINS'; 'FINS'; * *------------ 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' ISOL ; 'SI' (IGRD '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' IDYN; ZRMAS = 4.D0 '/' (WTAB.'DT' '**' 2) '*' WTAB.'MASSE' ; ZRAID = ZRMAS 'ET' ZRAID; 'SI' ('NEG' WTAB.'AMORTISSEMENT' 'INCONNU'); ZRAID = WTAB.'AMORTISSEMENT' '*' (2.D0 '/' WTAB.'DT') 'ET' ZRAID; 'FINS'; 'FINS' ; * 'SI' ISOL ; ZRAID =-1.* WTAB.'DT'* WTAB.'TETA'* WTAB.'PERMEABILITE' 'ET' ZRAID ; 'FINS' ; * WTAB.'OPERATEUR'= ZRAID ; 'FINS'; * 'SI' WTAB.'FOR_MECA' ; * * -------------------------------------------------------------------- * Conditions de contact-frottement 'SI' WTAB.'CONTACT'; MODCON = WTAB.'MODCONTA'; 'SINON' ; 'FINS'; * * Separer ce qui concerne l'adherence (FADH) et les jeux (FLX) 'SI' ('NEG' CJEU 0) ; 'FINSI' ; * 'SI' WTAB.'MODAL' ; 'NATURE' 'DISCRETE' ; 'SI' ('EGA' 1 &BCDA) ; CCDA = CHCR ; 'SINON' ; CCDA = CHCR 'ET' CCDA ; 'FINS' ; 'FIN' BCDA ; CDAP = CCDA ; 'FINS' ; * 'SI' ('NEG' CRR 0) ; ZCLIM = CRR 'ET' ZCLIM ; ZRAID = CRR 'ET' ZRAID ; 'FINS'; * 'SI' ('NEG' CJEU 0) ; ZFCONT = ZFCONT '+' CDAP ; 'FINS'; * 'SI' (('NEG' CRR 0) 'OU' ('NEG' CJEU 0)) ; ZFCONSTA = ZFCONSTA '+' ZFCONT ; 'FINS'; * 'FINS'; * * -------------------------------------------------------------------- * Matrice tangente et fefp 'SI' IKTAN ; 'SI' IFEFP ; IKT_SAUV = VRAI ; 'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ; 'MESS' 'FEFP: Start with LASTKTAN' ; ZRIKTA = ZLASTKTAN ; ZRAID = ZCLIM 'ET' ZRIKTA ; 'SINON' ; ZRAID = ZRAID 'ET' ('KSIGMA' ZMODL ZSIG0 ZMAT2) ; 'FINS' ; 'SINON' ; 'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ; IKT_SAUV = VRAI ; 'SI' IPERT ; 'MESS' 'Matrice tangente par perturbation - ' 'SINON' ; 'MESS' 'Matrice tangente "coherente" - ' 'FINS' ; ZRIKTA = ZLASTKTAN ; 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 ZMAT2 'PREC' ZPREK 'DT ' DTTAN ZKTASYM ; ZRAID = ZCLIM 'ET' ZRIKTA ; 'FINS' ; 'FINS' ; 'FINS' ; 'FINS' ; 'FINS' ; 'SI' WTAB.'MAN' ; ZRAIDINI = ZRAID ; 'FINSI' ; * * -------------------------------------------------------------------- * Raideur geometrique (option K_SIGMA) 'SI' (IGRD 'ET' ('NON' HPP_EPS)) ; 'SI' (IKSIA 'ET' ('NON' IFEFP)) ; * Passage de GEOREF0 -> GEOM1 'FORM' GEOM1; KSIGTC = 'KSIGMA' ZMODL ZSIGKS ZMAT2; ZRAID = ZRAID 'ET' KSIGTC; 'FINS' ; 'FINS' ; *----------------------------------------------------------------------* * Pilotage automatique *----------------------------------------------------------------------* COEPI = 1.D0 ; * 'SI' IPILOT ; 'SI' WTAB.'AUTODEUX' ; COEPI = WTAB.'AUTOCOEF' ; 'FINS' ; ZAUTOREDU = 1.D0 ; RED2 = 0 ; 'FINS'; * 'FINSI' ; *FOR_MECA *----------------------------------------------------------------------* * Partie variable du chargement a imposer a TI * * (geometrie, materiaux, chargement en deformations) * *----------------------------------------------------------------------* XDENDEF = 0. ; FDENDEF = 0. ; DEPST0 = 0. ; DMSRT0 = 0. ; DFTHE = 0. ; DFDEF = 0. ; ZFSUIV = 0. ; * * - Corrections en cas de materiaux variables 'SI' WTAB.'MATVAR' ; 'SI' ('OU' IENDOM IVIDOM ICERAM); 'SINON'; 'FINS'; DDEF0 = XXX4 - XXX3; DEPST0 = -1.* DDEF0; 'FINS'; * * - En cas de chargement thermiques 'SI' ITHER ; 'SINON'; MCHTETA2 = ZTEMP2; 'FINS'; DTT = ETT '-' ETT0 ; XDENDEF = XDENDEF '+' ETT ; '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'; * * - Contraintes et forces equivalentes associees a ces deformations 'SI' (DEPST0 'NEG' 0); 'SI' ('OU' IENDOM IVIDOM ICERAM); 'SI' (XDENDEF 'NEG' 0) ; 'FINSI' ; 'SINON'; 'SI' (XDENDEF 'NEG' 0) ; 'FINSI' ; 'FINS' ; 'SI' (DMSRT0 'NEG' 0); DSIGT0 = DSIGT0 '+' DMSRT0 ; XDENDEF = XDENDEF '+' MSRTT ; 'FINS'; 'SI' (XDENDEF 'NEG' 0) ; 'FINSI' ; 'FINS'; * * - Chargement en deformations imposees 'SI' LOGDEF; 'SI' ('OU' IENDOM IVIDOM ICERAM); 'SINON'; 'FINS'; 'FINS'; * * - Chargement de pression suiveuse en grands deplacements 'SI' (LOGPRE 'ET' IGRD) ; MOP = WTAB.'MOD_PRE' ; 'SINON' ; 'FINS' ; ZFSUIV = ZFSUIV '+' ZFPEXTF ; 'FINS'; * * - Procedure UTILISATEUR : y a-t-il des forces non conservatives '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 TI ; PRECED.'ADDI_MATRICE' = faux; * * FP22 = F^suiv_n+1 'SI' ADDISEC2 ; FP22 = TFP22.'ADDI_SECOND' ; FP022 = 'COPIER' FP22 ; ZFSUIV = ZFSUIV '+' FP22 ; 'FINS'; ZRAID = ZRAID 'ET' TFP22.'ADDI_MATRICE'; 'FINS'; 'FINS'; * * A-t-on des C.L. unilaterales ? *----------------------------------------------------------------------* * Second membre RESIDU * * * * Calcul du premier residu : 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 (ZDEP0). * * * * | forces exterieures sans reactions - forces interieures | * * | F^ext_n+1 + DF^suiv + DF^ther + DF^defi - F^int_n | * * | |-------------------------| | * * | -1.*ZFPLO | * * | (composantes de forces FX FY FZ ...) | * * RESIDU = | | * * | increment des relations imposees Du^imp | * * | (composantes de depl. FLX) | * * * * A F^ext peuvent s'ajouter des termes supplementaires p.ex. en * * dynamique ou en poreux. * *----------------------------------------------------------------------* * ZFCONSTA = [ F^ext_n+1 ; u^imp_n+1 ] * 'SI' WTAB.'FOR_MECA' ; * * Forces externes deja equilibrees au debut du pas par B*sigma * ZF1 = F^int_n = B*sigma_n + K^cst*u_n *HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO) * transport de zsig0 sur config debut de pas ZSIGBS = ZSIG0; * Passage de GEOREF0 -> GEOM1 'SI' IGRD ; 'FORM' GEOM1; 'FINSI'; 'SI' IRCON; 'FINS'; * 'SI' IDYN ; FFDYN = 'COPIER' ZF1; 'FINS'; * 'SI' ISOL ; XXXS =((1.- WTAB. 'TETA' )*GRAP0)+ (WTAB. 'TETA' * XXX1); XXX3 = ZF1 ; ZF1 = XXX3 - XXX2; 'FINS'; * * [ FLXINI ; FREAP ] = [ u^imp_n ; -F^reac_n ] FZDEP0 = ZDEP0 '*' ZCLIM; * * XXX1 = [ F^ext_n+1 ; u^imp_n+1 - u^imp_n ] XXX1 = ZFCONSTA '-' FLXINI ; * * FEXT0 : chargement externe (sans reactions) au pas precedent * FEXT0 = F^int_n - F^reac_n = F^ext_n FEXT0 = ZF1 '+' FREAP; * * 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 * * RESIDU = [DF^tot ; Du^imp] * RESIDU = [F^ext_n+1 + F^suiv_n+1 + DF^ther + DF^defi - F^int_n ; Du^imp] RESIDU = XXX1 '+' DFTHE '+' DFDEF '-' ZF1; RESIDU = RESIDU '+' ZFSUIV ; * * ZFPLO = -1.* [- F^int_n + DF^ther + DF^defi ] ZFPLO = ZF1 '-' DFTHE '-' DFDEF ; * * ZDFINI = [F^ext_n+1 + DF^ther + DF^defi - (F^int_n - F^reac_n) ; Du^imp] ZDFINI = XXX1 '-' FEXT0 ; * * Inc. de forces et de deplacements en distinguant CL unil et autres 'SI' IPILOT ; 'SI' IMPO12 ; 'FINS'; 'FINSI' ; * 'SI' (ITHER 'OU' WTAB.'MATVAR'); ZDFINI = ZDFINI '+' DFTHE; 'FINS'; 'SI' LOGDEF; ZDFINI = ZDFINI '+' DFDEF; 'FINS'; * ZSDMBR = RESIDU '*' 1.D0; * 'FINSI' ; * FOR_MECA * 'SI' WTAB.'NVSTNL' ; ZVIFL = conti.'VITESSES_FLUIDE' ; ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' (WTAB.'NSMA' '*' 1.5) ; fmass = WTAB.'NSMA' '*' ((2. '*' un) '-' (0.5 '*' unm)) ; 'SINON' ; ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' WTAB.'NSMA' ; fmass = WTAB.'NSMA' * un ; 'FINSI' ; ftnsl = ZFCONSTA 'ET' fmass ; 'SI' LOGADV ; ZRNSL = ZRNSL 'ET' nuadv ; fadv = nugrad * UADV ; ftnsl = ftnsl '-' fadv ; 'FINSI' ; resnsl = ftnsl '-' (ZRNSL * ZVIFL) ; * resolution ZVIFL = ZVIFL '+' ZVIF1 ; 'FINSI' ; *----------------------------------------------------------------------* * 1ere Resolution - ZDEP1 = [ Du^0 ; LX_n+1 ] * * * * -> estimation a TI a partir de l'instant TEMPS0 * * * * K*Du^0 + At*LX_n+1 = F^ext_n+1 - Bsigma_n * * A*Du^0 = d1 - A*u_n * * |----------------| |-------------------| * * ZRAID * ZDEP1 = RESIDU * *----------------------------------------------------------------------* 'SI' WTAB.'FOR_MECA' ; 'SI' WTAB.'ADHERENCE' ; 'SINON' ; FADRE = FADHE ; 'FINSI' ; FEXCI = FEXCI 'ET' FADRE ; 'FINSI' ; * * Modification du residu si automatique 'SI' LOGPIL ; RESIDU = RESIDU 'ET' (ETA0 '*' ZFPILIN); 'FINSI' ; * * Force limite de frottement 'SI' (WTAB.'CAFROTTE' 'ET' IMPO12); 'SI' WTAB.'FROCABL' ; FEXCI = FEXCI '+' FFROT ; 'FINSI'; 'SI' WTAB.'FROCOUL' ; excfconv = excfconv + 1 ; FEXCI = FEXCI '+' FFROT ; 'FINSI'; 'FINS'; * 'SI' ('NEG' ZRIBLO_M 'INCONNU') ; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12); ZRIBLO_M ZLISEA_M FEXCI RFNS RTRSF ; 'SINON'; ZRIBLO_M ZLISEA_M ; 'FINS'; 'SINON'; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12); 'SINON'; 'FINS'; 'FINS'; * 'SI' LOGPIL ; ZDEPILO = 0. * ZDEPII ; * Calcul de D_eta par appel a la procedure PILOINDI * Mise à jour ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ; ETA = ETA0 ; 'FINSI' ; * * Pour initialiser le statut des CL unil lors du prochain appel a RESO * Conserver les conditions unilaterales pour VITEUNIL 'SI' IDYN; WTAB.'ZRAIDV' = ZRAID; 'FINS'; ZRIBLO_M = ZRAID_T. 7 ; ZLISEA_M = ZRAID_T. 6 ; 'FINS'; * * Initialisation zdepl (sert pour xnum) *----------------------------------------------------------------------* * Calcul d'une norme pour la convergence * *----------------------------------------------------------------------* XXX1 = ZFEXT '+' ZFSUIV ; 'SI' LOGPIL ; XXX1 = XXX1 '+' ((ETA '+' D_ETA) '*' ZFPILIN); 'FINSI' ; * 'SI' ('EGA' knoconv 1) ; ZDEP1P50 = ZDEP1 ; 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET')) MLPRIM MLDUAL; XDENO1 = 'ABS' XDENO + (MZFM * MZDEP1M); MZDEP1M = MZDEP1M + XPETIT; XDENO=XDENO1/MZDEP1M; XDENO = XDENO + MZFM; 'SI' (XDENO < XPETIT); XDENO = 1.; 'FINSI'; XDENOM=XDENO; 'SI' TSTMOM ; ZDEP1P50 = ZDEP1 + XPETIT; ; XDENOM = XDENOM + XPETIT ; 'FINS' ; 'FINS' ; * ** mess 'xdeno propose' ' ' xdeno ' xdeno precedent' ' ' xdenoo; ** XDENO = MAXI XDENOo XDENO; ** XDENOo = XDENO; ** XDENOM = MAXI XDENOMo XDENOM; ** XDENOMo = XDENOM; * ZINCREMENT = ZDFINI '+' ZFSUIV ; 'FINS'; * RESIDNOR = 'COPIER' RESIDU ; * 'SI' IPILOT; * Objectif non atteint : conserver le XDENO 'SI' WTAB.'AUTODEUX' ; XDENO = ZXDENO ; XDENOM = ZXDENOM ; 'FINS'; * * Forces/deplacements en fin de pas DFEXT = COEPI '*' DFEXT0F ; ZFEXT = DFEXT '+' FEXT0 ; 'SI' IMPO12; DUIMP = (COEPI '*' DUIMPO) '+' DUUNIL ; 'SINON' ; DUIMP = COEPI '*' DFEXT0L ; 'FINS'; ZFLX1 = DUIMP '+' FLXINI ; * XFORC = ZFEXT '+' (COEPI '*' ZFSUIV) '-' ZF1 ; RESIDU = XFORC '+' DUIMP ; 'FINS'; *----------------------------------------------------------------------* * Corriger le residu a partir du pas precedent * *----------------------------------------------------------------------* * petite correction du residu pour esperer gagner du temps *********** INIT = FAUX ; 'SI'(('NEG' ZFNONL 'INCONNU') 'ET' WTAB.'INITIALISATION'); 'SI' IPILOT; 'SI' (WTAB.'AUTODEUX' 'ET' (COEPI 'NEG' 1.D0)) ; 'MESS' 'Initialisation a partir du pas precedent '; INIT = VRAI ; RESIDU = RESIDU '+' ZFNONL ; 'FINS' ; 'SINON'; * Faire la correction si le pas precedent a converge sans non convergence * Faire la correction si le pas precedent etait non lineaire. * on enleve le residu du pas precedent pour recuperer l'increment nominal du * second membre a imposer et l'increment du second membre du pas precedent ZINCREMENT = ZINCREMENT - WTAB.'RESIDU' ; zdeps = WTAB.'ZDEP1' + zdep1; AMPL = f12/(f22 + XPETIT); AMPLT = 0; 'SI' ('NEG' WTAB.'DTPREC' 0); DTPREC = WTAB.'DTPREC' ; 'SI' (DTPREC > XPETIT); AMPLT = WTAB.'DT' '/' 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 = WTAB.'DTPREC' ; AMPL=AMPLT; * la decharge est-elle significative 'SI'((F12/(f22 + XPETIT)) < -0.05);ampl=0.;'FINS'; ** 'MESS' 'F12 F22_' f12 f22; 'MESS' 'Pas d increment de charge, initialisation calculee avec le temps' ; 'SINON'; AMPL = F12 / (F22 + XPETIT); 'FINS'; * changement de modele on n'initialise pas 'SI' ('NEG' ZMODLI WTAB.'MO_TOT_PREC'); AMPL = 0; 'FINSI'; 'SI' ((AMPL > 0) 'ET' (AMPL < 2e1)) ; 'MESS' 'Initialisation a partir de la solution precedente Coeff'AMPL; XXX1 = AMPL * ZFNONL ; XXX2 =RESIDU+ XXX1; RESIDU = XXX2; INIT = VRAI; 'FINS'; 'FINS'; 'FINS'; 'FINS'; 'FINS'; WTAB.'ZDEP1'=zdep1; * * Initialisation pour l'acceleration de convergence (on peut mettre * n'importe quoi c'est pour ne pas faire de tests dans la boucle) ACFP1 = 'COPIER' ZFEXT2 *0. ; ACFP2 = ACFP1 ; ACFP3 = ACFP1 ; ACFEP1 = ACFP1 ; ACFEP2 = ACFP1 ; FCORF = 'COPIER' FREAP ; CORREC = 0; * *----------------------------------------------------------------------* * Initialisation des messages pour le pas de temps courant * *----------------------------------------------------------------------* 'SI' IPILOT ; 'MESSAGE' 'Iter'*13 'Nplas'*26 mocrit*39 'Deps.max'*52 'Eps.max'*65 moflex*78 'Alpha'*91; 'SINON' ; 'MESSAGE' 'Iter'*13 'Nplas'*26 mocrit*39 'Deps.max'*52 'Eps.max'*65 moflex*78 ; 'FINSI' ; 'FINSI' ; *FOR_MECA * NONCONV = FAUX; ZICONV = VRAI; PASTEST = FAUX; PASREINI = VRAI; ITACC = 0 ; RECA_N = 0 ; IT_RECA = 0 ; RED_URG = 0 ; DPSMAXP = 1; DEPSTDM = 0. ; MMCMAX = 0 ; IT = 0 ; MMC = 0 ; XCONV = 0. ; DPSMAX = xpetit ; EPSM = 0. ; IPREM = VRAI ; IRATE = FAUX; URG = FAUX ; *######################################################################* *----------------------------------------------------------------------* * Boucle de convergence ETIQ (Configuration GEOM1) * *----------------------------------------------------------------------* * preparation pour la gestion du pas avant etiq depstp zdeptp zsigfp fcorfp correcp = depst zdept zsigf fcorf 0; * 'REPETER' ETIQ ; * * IT est le compteur de ETIQ IT = IT + 1 ; * 'SI' WTAB.'FOR_MECA' ; * depst zdept zsigf fcorf correc = depstp zdeptp zsigfp fcorfp correcp; * * zdep1d doit etre coherent avec residu 'SI' IPREM; 'SINON' ; 'SI' ('NON' IRATE) ; 'FINSI'; 'FINSI'; * * Cas ou l'iteration precedente ne s'est pas faite entierement * -> partir du bon etat (mis a jour si le sous increment est acceptable) 'SI' IRATE; HPP_EPS = FAUX ; URG = VRAI ; ITACC = 4 ; 'SI' AUTAUG ; IRAUG = VRAI ; 'FINSI' ; 'FINSI' ; IRATE = VRAI; * * Champ materiau final ZMAT2F = ZMAT2 ; * 'SI' autaug; resmul = 2. * resmul; 'SI' ((dpsmax < 1d-3) 'ET' (resmul > 1D-20) 'ET' ('NON' PASTEST)); resmul = resmul / (dpsmax + xpetit) * 1d-3; 'FINSI'; 'SI' (resmul > 1.0); resmul = 1.0; 'SINON'; itacc = 4; 'FINSI'; 'SI' (resmul < 1.d-40); 'ERREUR' 996; 'FINSI'; 'FINSI'; * * * Pas de temps courant (WTAB.'DT') eventuellement sous-decoupe via COEPI 'SI' (knoconv '>' 1) ; ZDT = 0. ; 'SINON' ; ZDT = WTAB.'DT' '*' COEPI ; 'FINS' ; * * ITACC doit etre =< 0 pour qu'on accelere ITACC = ITACC - 1; * PASTEST = FAUX; RECA_K = FAUX; * * Pour conserver les criteres de convergence tabconv. it = 1; *----------------------------------------------------------------------* * Recalcul de la rigidite si necessaire * *----------------------------------------------------------------------* 'SI' IGRD ; 'FORM' GEOM1; * * URG = FAUX ; IT_RECA = IT; RECA_K = VRAI; RECA_N = RECA_N + 1; HPP_EPS = FAUX; 'SI' (RECA_N > 20) ; nonconv = vrai; 'FINS'; PASREINI = FAUX; * * A cause de FORM, ZMATTEMP n'est plus parallele... * ZDEPTR = ZDEP0 '+' ZDEPT ; * * ---------------------------------------------------------------- * Recalcul de la rigidite a la fin du pas 'SI' ('EGA' LAG_TOT 1) ; 'FORM' GEOREF0; 'FORM' GEOR; 'SINON'; 'FINSI'; 'FORM' GEOM1; ZMATTEMP = 0 ; 'SI' (SOUCI) ; 'FINSI'; * * RITC peut contenir des CL ZCLIM = ZCLIM0 'ET' ZCL ; 'SINON' ; ZCLIM = ZCLIM0 ; 'FINSI' ; * * ---------------------------------------------------------------- * Option 'AUGMENTATION_AUTOMATIQUE' : actualiser la rig augm 'SI' (AUTAUG 'ET' ('NON' IRAUGLU)); 'FORM' GEOR ; RIG_AUG = 'MASSE' ZMODL ZMAT2F ; RH = RITC; 'FORM' GEOM1; *** mess 'actualisation rig_aug'; 'FINSI'; * * ---------------------------------------------------------------- * Ajout de relations issues de la procedure DEFO_IMP ZFCONT = 0 ; IACTURES = FAUX ; ** DE_CNTRL = FAUX; 'SI' (AUTAUG 'ET' DE_CNTRL); IACTURES = VRAI ; ZCLIM = RE_CNTRL 'ET' ZCLIM ; 'FINSI'; * * ---------------------------------------------------------------- * Recalcul des conditions de contact-frottement ds la conf finale 'SI' WTAB.'CONTACT'; 'FORM' GEOR ; * * MODCON = wtab.'MODCONTA'; 'SINON' ; 'FINS' ; * * Separer ce qui concerne l'adherence (FADH) et les jeux (FLX) 'SI' ('NEG' CJEU 0) ; * * Ajout du glissement deja parcouru sur les conditions de frot. 'SI' WTAB.'CAFROTTE'; CDAP = CDAP ET (-1*CGLI) ; 'FINSI'; 'FINSI'; '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' ; * 'SI' ('NEG' CRR 0) ; IACTURES = VRAI ; 'SI' ('NEG' ZFCONT 0) ; ZFCONT = ZFCONT '+' CHPZ ; 'SINON' ; ZFCONT = CHPZ ; 'FINSI' ; ZCLIM = CRR 'ET' ZCLIM ; 'FINS'; * 'SI' ('NEG' CJEU 0) ; IACTURES = VRAI ; ZFCONT = ZFCONT '+' CDAP ; 'FINS'; 'FINS'; * * Mise a jour des jeux dans le residu (partie force inchangee) ZFCONSTA = ZFEXT2 ; ZFLX1 = ZFLXB ; 'SI' IACTURES ; ZFCONSTA = ZFCONSTA '+' ZFCONT ; RJEUX = ZFCONT '+' ZFLXB '-' FLXPRE ; RESIDU = RFORCE 'ET' RJEUX ; ZFLX1 = ZFLX1 'ET' ZFCONT ; 'FINS'; * 'SI' (IPILOT 'ET' IMPO12) ; 'FINS'; * 'DETR' ZRAID; ZRAID = ZCLIM 'ET' ZRI; * Mise a jour de FREAP (faut-il faire cette actualisation?) * * ---------------------------------------------------------------- * Ajout de la partie dynamique 'SI' IDYN; * bp : en toute rigueur, il faudrait aussi recalculer la MASSE ... * et ajouter l'amortissement le cas échéant ... ZRAID = ZRAID 'ET' ZRMAS ; 'FINS'; * * ---------------------------------------------------------------- * Recalcul de la raideur geometrique (option K_SIGMA) 'SI' IKSIA ; 'FORM' GEOR; KSIGTC = 'KSIGMA' ZMODL ZSIGKS ZMAT2; ZRAID = ZRAID 'ET' ksigtc ; 'FORM' GEOM1; 'FINS'; * * ---------------------------------------------------------------- * Prise en compte d'eventuelles 'RIGIDITE_CONSTANTE' 'SI' IRCON; ZRAID = ZRAID 'ET' RIG_CONS; 'FINS'; * * ---------------------------------------------------------------- * Traitement en de la rigidite augmentee 'SI' IRAUG; 'SI' AUTAUG ; 'FINSI'; augauto = xkx / xmx; augk = 1.d-4 ; augk = augk * (1. + (wtab.'ISOUSPAS' * 2.) ); augauto = abs augauto; augk = abs augk; augmult = augmult * 1.02; 'SI' (augmult > 1D2 ); augmult=1D2 ; 'FINSI'; augauto = augauto * augmult ; augk = augk * augmult; ZRAID = ZRAID 'ET' (RIG_AUG * augauto) 'ET' (RH * augk) ; 'MESS' 'multiplicateur d augmentation masse' augauto ' raideur' augk ; 'SINON'; ZRAID = ZRAID 'ET' RIG_AUG ; 'FINSI'; *** ZNACCE = FAUX; 'SINON'; augmult = augmult * 0.55 ; si (augmult < 1d-3); augmult = 1d-3; finsi; 'FINSI'; 'FORM' GEOM1; * * ---------------------------------------------------------------- * Procedure UTILISATEUR 'SI' WTAB.'PROCEDURE_CHARMECA'; PRECED.'ADDI_MATRICE' = vrai; TFP22= CHARMECA PRECED TI ; PRECED.'ADDI_MATRICE' = faux; zraid = zraid 'ET' TFP22.'ADDI_MATRICE'; 'FINS'; 'FINS'; * * * on impose le recalcul de K a la prochaine iteration si it=2 (pq ?) 'SI' (IT 'EGA' 2) ; ITACC = 4 ; 'FINS' ; 'FINS'; 'FINS'; *----------------------------------------------------------------------* * Evaluation de la matrice tangente si demandee * *----------------------------------------------------------------------* 'SI' IKTAN ; 'SI' ('NON' IFEFP) ; * IKT = VRAI ; 'SI' (IGRD 'ET' ('NON' HPP_EPS)); 'FORM' GEOR ; 'FINS' ; 'SI' ('NON' ISSTE) ; 'SI' IPERT ; 'SI' PARTLOCA ; 'C1' WTAB.'K_TANG_PERT_C1' 'C2' WTAB.'K_TANG_PERT_C2' ZKTASYM ; ZRIKTA = 'ET' zktap ; 'SINON'; 'C1' WTAB.'K_TANG_PERT_C1' 'C2' WTAB.'K_TANG_PERT_C2' ZKTASYM ; 'FINS'; 'SINON' ; 'SI' ('NON' (IVISCO 'OU' IVIDOM)) ; DTTAN = 0. ; 'SINON' ; DTTAN = ZDT ; 'FINS' ; ZRIKTA = 'KTAN' ZMODL ZSIGF ZVARF ZMAT2F 'PREC' ZPREK 'DT ' DTTAN ZKTASYM ; 'FINS' ; 'FINS' ; 'SI' IKSIA ; ZKSIG = 'KSIGMA' ZMODL ZSIGF ZMAT2F ; ZRIKTA = ZRIKTA 'ET' ZKSIG ; 'FINS' ; ZRAID = ZRIKTA 'ET' ZCLIM ; 'FINS' ; 'FINS' ; 'FINS' ; 'FINSI' ; *FOR_MECA * 'SI' WTAB.'NVSTNL' ; ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' (WTAB.'NSMA' '*' 1.5) ; 'SINON' ; ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' WTAB.'NSMA' ; 'FINSI' ; ftnsl = ZFCONSTA 'ET' fmass ; 'SI' LOGADV ; ZRNSL = ZRNSL 'ET' NUADV ; ftnsl = ftnsl '-' fadv ; 'FINSI' ; resnsl = ftnsl '-' (ZRNSL * ZVIFL) ; 'FINSI' ; *----------------------------------------------------------------------* * Acceleration de convergence effective * *----------------------------------------------------------------------* 'SI' WTAB.'FOR_MECA' ; CORRECA = CORREC ; CORREC = 0; ACFEP0 = RESIDU - FCORF; ACFEPT = ACFEP0; ACFP0 = ACFEP0 * ZJAC ; ACFEP0 = ACFEP0 - CORRECA; 'SI' ((ITACC '<EG' 0) 'ET' (IT '>' 3) 'ET' ZNACCE); ITACC = 2; CORREC = 'ACT3' ACFEP2 ACFEP1 ACFEP0 ACFP3 ACFP2 ACFP1 ACFP0 ; * * verif que l'acceleration ne renvoie pas en arriere 'SI' (WTAB.'STABILITE' 'OU' IPILOT) ; acc_ref = acc_ref + xpetit; acc_rap = acc_dir/acc_ref; acc_lim = 1E5 ; 'SI' (acc_ref < 0.); acc_lim = 1.; 'FINSI'; 'SI' (acc_rap '>' acc_lim) ; * en cas d'acceleration trop grande, on la limite 'FINSI'; 'SI' (acc_rap '<EG' 0.) ; * en cas d'acceleration retrograde, on n'accelere pas 'MESS' 'Annulation acceleration: retrograde' ' ' acc_rap; correc = 0.; 'FINSI'; 'FINSI'; * RESIDU = RESIDU '-' CORREC; 'FINS'; * ACFP3 = ACFP2 ; ACFP2 = ACFP1 ; ACFP1 = ACFP0 ; ACFEP2 = ACFEP1 ; ACFEP1 = ACFEP0 ; * 'SI' (resmul < 0.99); residu = ((residu -fcorf) * resmul) + fcorf; ACFP0 = (RESIDU - FCORF) ; 'FINSI'; *----------------------------------------------------------------------* * Systeme a resoudre a l'iteration (i) * * * * | K At | | du | | Fsur + Fvol - Bsigma^(i) | * * | | | | = | | * * | A 0 | | LX(i) | | d1 - A*u^(i) | * *----------------------------------------------------------------------* 'SOUCI' 0; 'SI' WTAB.'ADHERENCE' ; 'SINON' ; FADRE = FADHE ; 'FINSI' ; FEXCI = FEXCI 'ET' FADRE ; 'FINSI' ; * * Force limite de frottement 'SI' (WTAB.'CAFROTTE' 'ET' IMPO12); 'SI' WTAB.'FROCABL' ; FEXCI = FEXCI '+' FFROT ; 'FINSI'; 'SI' WTAB.'FROCOUL' ; excfconv = excfconv + 1 ; FEXCI = FEXCI '+' FFROT ; 'FINSI'; 'FINS'; * 'SI' ( ('NEG' ZRIBLO_M 'INCONNU') 'ET' (IPREM 'OU' RECA_K) ) ; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12); ZRIBLO_M ZLISEA_M FEXCI RFNS RTRSF; 'SINON'; ZRIBLO_M ZLISEA_M ; 'FINS'; 'SINON'; 'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12); 'SINON'; 'SI' (WTAB.'MAN' 'ET' IPREM); ORDRE = WTAB.'ORDRE' ; ZDEP2 IOUT = CORMAN ZRAIDINI ZMODL ZMAT2F ORDRE ZDEP0 ZSIG0 RESIDNOR WTAB ; 'SI' (IOUT 'EGA' 1) ; ZDEP1=ZDEP2; 'FINS'; 'FINSI' ; 'FINS'; 'FINS'; * * verif stabilite matrice elastique 'SI' (wtab.'STABILITE' 'ET' ('NON' HPP_EPS) 'ET' AUTAUG); ** diag refactorise la matrise. On evite donc provisoirement de l'appeler ** nbng = 'DIAG' ZRAID; nbng = 0; 'SI' ((acc_ref < 0.) 'OU' (nbng > 0)); 'MESS' 'Raideur negative' ' ' acc_ref ; ** resmul = resmul * 0.25; ** augmult = augmult * 1.4; 'SI' LOG_CNTRL; * controler la deformation max RE_CNTRL2*'RIGIDITE' = DEFO_IMP DEPST (zdept + zdep1) ZRAID ; DE_CNTRL = VRAI; ** cmul = MINI 1. ( 1d-4 / DPSMAX); ** ZFEXT2 = ZFEXT2 + (ZFCONT * cmul); ZFCONSTA = ZFEXT2 ; RE_CNTRL = RE_CNTRL 'ET' RE_CNTRL2; 'FINSI'; ** 'ITERER' ETIQ; ** essai en inversant zdep1 pour avoir la positivite du travail 'FINSI'; 'FINSI'; * 'SI' ('SOUCI') ; 'SI' AUTAUG ; resmul = resmul * 0.25; augmult = augmult * 1.4; 'FINSI'; pastest = VRAI; HPP_EPS = FAUX; URG = VRAI; ITACC = 4; *** 'ITERER' ETIQ; 'FINSI'; * 'SI' LOGPIL ; * Calcul de D_eta par appel a la procedure PILOINDI * Mise à jour ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ; ETA = ETA '+' D_ETA ; 'FINSI' ; * 'SI' (resmul < 1d-2) ; souci 0; 'FINSI'; pasunil = 'SOUCI' ; * * Conserver les conditions unilaterales pour VITEUNIL 'SI' IDYN ; WTAB.'ZRAIDV' = ZRAID; 'FINS'; ZRIBLO_M = ZRAID_T.7; ZLISEA_M = ZRAID_T.6; 'FINS'; 'FINSI' ; *FOR_MECA 'SI' WTAB.'NVSTNL' ; ZVIFL = ZVIFL '+' ZVIF1 ; 'FINSI' ; *----------------------------------------------------------------------* * Increment de deplacements et mult. de Lagrange * * * * ZDEP1 est un increment d'increment de deplacements, note du, et les * * multiplicateurs de Lagrange sont resolus "en total". * * La solution de la boucle ETIQ a l'iteration i a fourni : * * ZDEP1 = [ du ; LX(i) ] * * Le cumul est realise de maniere a avoir : * * ZDEPT = [ Du^(i) ; LX(i) ] * * avec Du^(i) = Du^(i-1) + du * * Le champ de deplacements total s'ecrit : * * ZDETOT = ZDEP0 '+' ZDEPT * *----------------------------------------------------------------------* 'SI' WTAB.'FOR_MECA' ; 'SI' IPREM ; ZDEPT = 'COPIER' ZDEP1 ; ZDELA = 'COPIER' ZDEPT ; 'SINON'; XXX1 = ZDEPT 'ENLEVER' 'LX' ; ZDEPT = XXX1 '+' ZDEP1 ; 'DETR' XXX1 ; 'FINS' ; * 'SI' WTAB.'MODAL' ; 'SI' IPREM ; 'SINON' ; * mettre les point materiels dans zdept 'FINS'; 'FINS'; 'FINS' ; * * Option automatique : determiner le coef de normalisation 'SI' IPILOT; PASTEST = '<' IT WTAB.'AUTORECA' ; 'SI' ACCEL ; * * MODIFICATION CB215821 : 18/06/2015 'SI' (RED2 '<' (IT '/' 20)) ; RED2 = RED2 '+' 1 ; ZAUTOREDU = ZAUTOREDU '*' 3.D0 ; 'FINS'; 'SI' (ZAUTOREDU '>' 1.D0 ) ; 'MESS' 'On divise le critere de pilotage par 'ZAUTOREDU; 'FINS' ; * OO = WTAB.'AUTOCRIT' '/' ZAUTOREDU ; AL1 = OO '/' U1MA ; * * AL1 : coefficient de normalisation 'SI'((AL1 '>' 1.D0) 'ET' (COEPI '>' 0.D0)); 'SI' (AL1 '>' (1.D0 '/' COEPI)); AL1 = 1.D0 '/' COEPI; 'FINS'; 'FINS'; * * Normalisation XXX1 = (1.D0 '-' AL1) '*' ZLX ; XXX3 = AL1 '*' ZDEPT ; ZDEPT = XXX3 '+' XXX1; 'DETR' XXX3; 'FINSI' ; * 'SINON' ; * On part dans le decor : redemarrer a 0 'SI' ((XCONV > 1E8) 'ET' PASREINI) ; 'MESS' 'Reinitialisation du schema'; PASREINI=FAUX; 'ITERER' ETIQ; '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 *---------------------------------------------------------------------- 'DETRUIRE' FCORF; * * Champ de deplacement total ZDETOT = ZD0SLX '+' ZDEPT ; ZDETOTN = -1. '*' ZDETOT ; * * Gradient des deplacements (total) a la fin du pas de temps * par rapport a la configuration de reference GEOREF0 'SI' IGRD; 'SI' ('NEG' ZGRDU0 'INCONNU'); 'FORM' GEOREF0 ; D_GR_U = ZGRDUF '-' ZGRDU0 ; 'FORM' GEOM1 ; 'FINS'; 'FINS'; * * XXX1 = [ FCORF ; FCORU ] = [ -Freac_it ; u^imp_it ] XXX1 = ZCLIM '*' ZDETOT; 'DETR' XXX1 ; *----------------------------------------------------------------------* * Increment de deformations totales DEPST * *----------------------------------------------------------------------* * Update or total lagrangian --------------------------------------- 'SI' IFEFP; 'SI' IFEFPUL ; * mess ' update lagrangian ZRIKTA'; ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP' ZMODL ZDEI0 ZVAR0 ZDEPT ZMAT2F ZPREK NITMA 1 ; 'SINON'; * mess ' total lagrangian ZRIKTA'; chp_z = ZDEPT + ZDEP0 ; ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP' ZMODL ZDEI0 ZVAR0 chp_z ZMAT2F ZPREK NITMA ; chp_z = 1 ; 'FINS'; ZRAID = ZRIKTA 'ET' ZCLIM ; 'FORM' GEOM1 ; * ACC = 'ABS' ( XXX1 - ACC0 ) ; MMC = 'MASQUE' ACC 'SUPERIEUR' 1.D-10 'SOMME' ; 'SI' (MMC '>' MMCMAX) ; MMCMAX = MMC ; 'FINS' ; ZDEFF = ZDEF0 '+' DEPST ; * * cas standard ----------------------------------------------------- 'SINON'; * ZDEF0a = ZDEF0 ; ZSIG0a = ZSIG0 ; ZVAR0a = ZVAR0 ; ZDEI0a = ZDEI0 ; DEPST = 0. ; DSIGT = 0. ; * nsoincr = nsoincrn ; ZDEPTI = ZDEPT '/' nsoincr ; * * --------------------------------------------------------------- * Sous-incrementation du comportement 'REPETER' sousinc nsoincr; * codeb = 'FLOT' (&sousinc '-' 1) '/' nsoincr ; cofin = 'FLOT' &sousinc '/' nsoincr ; * 'SI' IGRD; ZDEPF = ZDEPT '*' cofin ; ZDEPM = ZDEPT '*' ((&sousinc - 0.5) '/' nsoincr) ; 'FORM' GEOM1 ; 'SI' WTAB.'ITCAR'; * A cause de FORM, ZMATTEMP n'est plus parallele... 'SINON'; ZMATTEMP = ZMAT2F; 'FINS'; 'FORM' GEOM1 ; 'FINS'; * * -------------------------------------------------------------- * Calcul de l'increment de deformations 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) iDEFOK = FAUX ; 'SI' (('NEG' HYPDEF 'LINEAIRE') 'ET' IGRD); * 'FORM' GEOM2M ; 'SI' ('OU' IENDOM IVIDOM ICERAM); 'DETR' HOOKENDO; 'SINON'; 'FINSI'; 'FORM' GEOM1; iDEFOK = 'NON' ('SOUCI') ; DPSMAX=1001 ; 'SI' iDEFOK ; 'FINSI' ; * 'SI' (('NON' iDEFOK) 'OU' (DPSMAX > 1000)) ; 'SI' autaug; augmult = augmult * 1.4; resmul = resmul * 0.25; 'ITERER' ETIQ; 'FINSI'; pastest = vrai; 'SINON'; * Increment de deformations transporte sur GEOREF0 'SI' IGRD ; 'SI' IJAUMA ; 'FORM' GEOM1; 'FORM'GEOREF0 ; 'FORM' GEOM1; 'SINON' ; 'FORM' GEOREF0; 'FORM' GEOM1 ; 'FINSI' ; 'FINSI'; iDEFOK = VRAI ; 'FINS'; 'FINS'; * * Calcul des deformations lineaires (demande ou pas reussi) 'SI' ('NON' iDEFOK) ; 'SI' ('OU' IENDOM IVIDOM ICERAM); 'DETR' HOOKENDO; 'SINON'; 'FINSI'; 'SI' (('NEG' HYPDEF 'LINEAIRE') 'ET' IGRD); 'MESS' 'Attention utilisation des deformations lineaires'; 'FINS'; * * Increment de deformations transporte sur GEOREF0 'SI' IGRD ; 'FORM' GEOREF0; 'FORM' GEOM1 ; 'FINSI'; 'FINS'; * ZDEFF = ZDEF0a '+' DEPSTa ; DEPSTZ = DEPST '+' DEPSTa ; 'DETR' DEPST ; DEPST = DEPSTZ ; * * -------------------------------------------------------------- * Transport des quantites sur la configuration adequate * 'SI' IGRD ; * Lagrangien total : sur GEOREF0 'SI' ('EGA' LAG_TOT 1); 'FORM' GEOREF0; 'FINSI'; * Lagrangien fin pas : GEOREF0 -> GEOM2 'SI' ('EGA' LAG_TOT 2); 'FORM' GEOM2; ZDEFIN = ZDEP0 '+' ZDEPF ; ZDEFINN = -1. '*' ZDEFIN ; 'SI' ('NEG' &sousinc 1) ; 'FINSI' ; ZSIG0a ZDEF0a DEPSTa = ZSIG0Z ZDEF0Z DEPSTZ ; 'FINS'; * Lagrangien mi pas : GEOREF0 -> GEOM2M 'SI' ('EGA' LAG_TOT 3); 'FORM' GEOM2M; ZDEMIL = ZDEP0 '+' ZDEPM ; ZDEMILN = -1. '*' ZDEMIL ; 'SI' ('NEG' &sousinc 1) ; 'FINSI' ; ZSIG0a ZDEF0a DEPSTa = ZSIG0Z ZDEF0Z DEPSTZ ; 'FINS'; 'FINS'; * 'SI' IDYN ; 'SI' (('NEG' nsoincr 1) 'ET' ('EGA' &sousinc 1)); ZSIGM0 = ZSIG0a ; ZSIGM = ZSIGM0 '/' 2; 'FINS'; 'FINS'; * * Verifier si les deformations ne sont pas trop importantes 'SI' IGRD ; ref = 20; * 'SI' ((DPSMAX > 1D-2) 'ET' (DPSMAXP < 1D-3)); REF = 100; 'FINSI'; 'SI' (DPSMAXP * 3 < DPSMAX); ref = 100; DPSMAXP = 1.1 * DPSMAXP ;'FINSI'; 'SI' WTAB.'RECALCUL_SOUSINC' ; ** mess 'sous decoupage pour epsi : ' nsoincrn; 'FINSI' ; coefmul = 1.D0 / (DPSMAX + XPETIT); 'SI' (((coefmul < ref) 'OU' PASUNIL) 'ET' AUTAUG); si LOG_CNTRL; RE_CNTRL2*'RIGIDITE' = DEFO_IMP DEPSTa (zdept + zdep1) ZRAID ; DE_CNTRL = VRAI ; ** cmul = MINI 1. ( 1d-4 / DPSMAX); ** ZFEXT2 = ZFEXT2 + (ZFCONT * cmul); ZFCONSTA = ZFEXT2 ; RE_CNTRL = RE_CNTRL 'ET' RE_CNTRL2; finsi; * il faut reduire d'urgence les deplacements 'SI' (red_urg > 2); 'FINSI'; augmult = augmult * 1.5; 'SI' (RED_URG > 2) ; 'SI' ('NON' nonconv); ZMAXIT = IT+5; 'FINS'; nonconv = vrai; 'FINS'; ** 'SI' (resmul < 0.99); ** resmul = resmul * ('MINI' 0.25 (1d-3 / dpsmax)); ** 'FINSI'; '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 > 3)); ZMAXIT = IT-1; 'FINS'; RED_URG = RED_URG + 1; 'ITERER' etiq; 'SINON'; ** RE_CNTRL = 'VIDE' 'RIGIDITE' / 'RIGIDITE'; ** DE_CNTRL = FAUX; 'FINS'; 'FINS'; * * -------------------------------------------------------------- * Loi de comportement sur la configuration definie par LAG_TOT * (Operateur ELAS ou COMP) * * - Comportement elastique lineaire 'SI' ('NON' IPLAVI) ; * * Increment de contraintes * * Soustraire les parties thermique et deformations imposees 'SI' (ITHER 'OU' WTAB.'MATVAR') ; XXX1 = DSIGT0 '*' (COEPI '/' nsoincr) ; XXX2 = DSIGTa '-' XXX1; DSIGTa = XXX2 ; 'FINS'; * 'SI' LOGDEF ; XXX1 = DSI1 '*' (COEPI '/' nsoincr) ; XXX2 = DSIGTa '-' XXX1; DSIGTa = XXX2 ; 'FINS'; * * - Comportement non lineaire : integration avec l'operateur COMP 'SINON'; * * 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 * * Soustraire les parties thermique et deformations imposees 'SI' (ITHER 'OU' WTAB.'MATVAR') ; XXX1 = DEPST0 '*' (COEPI '/' nsoincr) ; DEPST = DEPST '-' XXX1 ; DEPSTa = DEPSTa '-' XXX1 ; 'FINS' ; * 'SI' LOGDEF ; XXX1 = DDEFOR0 '*' (COEPI '/' nsoincr) ; DEPST = DEPST '-' XXX1 ; DEPSTa = DEPSTa '-' XXX1 ; 'FINS' ; * ZMATT = ZMAT2F ; 'SI' ('ET' ('ET' ITHER ('OU' IENDOM IVIDOM )) WTAB.'MATVAR') ; * on recupere certain materiau avec les parametres fct de la temperature * voir PAS_mate (il ne faut que la dependance thermique) 'FINS'; * * ...cas SSTE... 'SI' ISSTE; ZRIKTA ZSIGFa ZVARF ZDEIF = 'SSTE' ZMODL ZSIG0a ZVAR0 DEPSTa ZMATT ZPREK WTAB.'NMAXSUBSTEPS' NITMA; * * ...cas non SSTE... 'SINON'; * tsodeb = TEMPS0 '+' (ZDT '*' codeb) ; tsofin = TEMPS0 '+' (ZDT '*' cofin) ; 'SI' (knoconv '>' 1) ; tsodeb = tsofin ; 'FINS' ; * che1 = che1 che2 = che2 'FINS'; 'SI' ither ; TETDE = ZTET1 + (DTETD '*' codeb) ; TETDF = ZTET1 + (DTETD '*' cofin) ; 'FINS' ; che1 = che1 'ET' che3 ; che2 = che2 'ET' che4 ; *On met les deformations en tete des champs pour COMP ('KTAN' 'PERT') che1 = ZDEF0a 'ET' che1 ; che6 = ZDEF0a '+' DEPSTa ; che2 = che6 'ET' che2 ; * pour les materiaux on garde toujours la valeur a la fin du pas * et au debut du pas 'SI' ('EXISTE' PRECED 'MODELE_STATIONNAIRE') ; 'SINON' ; che1 = che1 'ET' ZMAT1 ; 'FINSI' ; * * 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' ZGRDU0 'INCONNU') 'ET' igrd) ; gru1 = ZGRDU0 + (D_GR_U '*' codeb) ; che1 = che1 'ET' gru1 ; gru2 = ZGRDU0 + (D_GR_U '*' cofin) ; che2 = che2 'ET' gru2 ; 'FINS'; * 'SI' LNLOC ; ISTEP = 1 ; Che1a = che1 ; che1 = che1 'ET' chm_z ; Che2a = che2 ; che2 = che2 'ET' chm_z ; 'FINS'; * che11 = che1 'ET' ZSIG0a 'ET' ZVAR0a 'ET' ZDEI0a ; che22 = che2 'ET' ZMATT ; 'SI' PARTLOCA; souci 0 ; isoucomp = souci ; cha11 = 0 ; cha22 = 0 ; cho22 = 0 ; 'SINON'; souci 0 ; isoucomp = souci ; 'FINS'; * *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * DEBUT MODELE NON LOCAL 'SI' LNLOC ; * NLOC ne traitant pas des champs //, on reduit tout au modele initial. 'SI' (PARALLEL 'ET' ('NON' PARTLOCA)) ; 'FINS' ; ********** boucle sous-iterations Helmholtz ************** * donnees pour acceleration Helmholtz * periode de recours a l acceleration ACT3 pour Helmholtz ZNACCEHL = 3; * NBR min de sous-iterations d Helmholtz a faire avant d appeler ACT3 pour la 1ere fois ITDEPHL = 4; 'REPE' BIHL WTAB.'MAXSOUSITERATION' ; 'SI' ('NEG' WTAB.'NON_LOCAL' 'HELM'); * cas non-local MOYE ou SB 'SI' ('EGA' WTAB.'NON_LOCAL' 'SB') ; MOD_SB = WTAB.'NLOC_MODL' ; ZVARF = ZVARF '+' CONTP '+' WTAB.'NLOC_SB_REGU' ; 'FINS' ; LOGHLIN = VRAI ; 'SINON' ; * cas non-local HELM Sellier Millard MOD_HELM = WTAB.'NLOC_MODL' ; * (re) lecture du materiau si init ou evolutif * (re) calcul des matrices de rigidite de Helmholtz Sellier PAS_HELM PRECED ZMATHL; * variable pour la liste des residus de convergence sur Helmholtz * ----- debut boucle sur les formulations de Helmholtz --- * on reste en lineaire si toutes les vari sont lineaires (sinon on aura des istep a 3) LOGHLIN = VRAI ; 'REPE' BH NHELM ; LOGHLIN=LOGHLIN 'ET' TAHELM . &BH. 'LINEAIRE' ; * nom de la variable a delocaliser LEMO = TAHELM . &BH. 'NOM' ; * recuperation de la variable sur le modele mecanique * projection sur le modele de Helmholtz associe * TRAC zvaux2 TAHELM . &BH . 'H_MODELE' titre ('CHAI' 'Vari a diffuser:'LEMO ); * la variable devient la source * trac fsour (extr TAHELM . &BH . 'H_MODELE' 'MAIL') titre 'Source'; * trac (TAHELM . &BH . 'H_OPER') titre 'H_Oper'; * *** pour istep=1 *** 1ere etape non locale *** 'SI' ('EGA' &BIHL 1) ; * &BIHL est le compteur de sous iteration de Helmholtz * il vaut forcement 1 si istep vaut 1, on fait la resolution classique * sauvegarde pour possible acceleration ulterieure TAHELM . &BH . 'XN' = ZVNEW ; 'SINON'; * recuperation de la solution de la sous iteration de Helmholtz precedente ZVNEW = TAHELM . &BH . 'XN' ; 'FINSI'; * *** traitement non lineaire de la formulation *** 'SI' ( 'NON' TAHELM . &BH . 'LINEAIRE' ) ; * calcul du residu sur l evolution de la source RESN = FSOUR - (TAHELM . &BH . 'H_OPER' '*' TAHELM . &BH . 'XN' ); * stockage du residu pour la mesure d erreur * normalisation des residus par l inverse de la capacite * (calculee dans pas_helm.procedur) ZERR0 = RESN0 '*' (TAHELM . &BH . 'INV_CAPA' ); * erreur residu * 'SI' ('EGA' &BIHL 1); * TAHELM . &BH . 'RES1' = ZERR2; * 'FINSI'; * carre de norme du residu * ZERR1 = 'XTX' ZERR0; * Residu * ZERR2 = (ZERR1)**0.5; * normalisation de la source par la capacite * calcul erreur * carre de la norme de la source * ZERR5 = 'XTX' ZERR3; * norme de la source * ZERR6 = ZERR5 **0.5; * valeur maximale du champs * ZERR3='MAXI' (TAHELM . &BH . 'XN') 'ABS'; * mess 'variable Helmholtz', LEMO; * mess 'residu', ZERR2; * mess 'norme de la source ', ZERR5; * mess 'norme Source', ZERR3; * on teste si la valeur maximale est capable de normaliser 'SI' (ZERR6 '>' 0.); ZERR4 = ZERR2 '/' ZERR6; 'SINON'; ZERR4 = ZERR2; 'FINS'; * preparation acceleration 'SI' ( 'EGA' &BIHL 1 ) ; TAHELM . &BH . 'H_RESO' = TAHELM . &BH . 'H_OPER' ; FWOR = FSOUR * 0. ; TAHELM . &BH . 'ACFP1' = FWOR ; TAHELM . &BH . 'ACFP2' = FWOR ; TAHELM . &BH . 'RESNP' = FWOR ; TAHELM . &BH . 'ACFP3' = FWOR ; TAHELM . &BH . 'ACFEP1' = FWOR ; TAHELM . &BH . 'ACFEP2' = FWOR ; TAHELM . &BH . 'CORREC' = 0.; 'SINON' ; * acceleration de convergence act3 CORRECPHL = TAHELM . &BH . 'CORREC'; CORRECHL = 0; ACFP0HL = RESN - TAHELM . &BH . 'RESNP' ; ACFEP0HL = ACFP0HL - CORRECPHL ; 'SI' (&BIHL > ITDEPHL); CORRECHL = 'ACT3' TAHELM . &BH . 'ACFEP2' TAHELM . &BH . 'ACFEP1' ACFEP0HL TAHELM . &BH . 'ACFP3' TAHELM . &BH . 'ACFP2' TAHELM . &BH . 'ACFP1' ACFP0HL ; RESN = RESN - CORRECHL ; 'FINS'; 'FINS'; TAHELM . &BH . 'RESNP' = RESN ; TAHELM . &BH . 'ACFP3' = TAHELM . &BH . 'ACFP2' ; TAHELM . &BH . 'ACFP2' = TAHELM . &BH . 'ACFP1' ; TAHELM . &BH . 'ACFP1' = ACFP0HL ; TAHELM . &BH . 'ACFEP2' = TAHELM . &BH . 'ACFEP1' ; TAHELM . &BH . 'ACFEP1' = ACFEP0HL ; * calcul correction delta-sigma ZVNEW = ZVNEW + DZVN ; TAHELM . &BH . 'XN' = ZVNEW ; 'FINSI' ; 'SINON'; * cette variable est traitee de facon lineaire ZERR4=0.; 'FINSI'; * trac zvnew mod_helm titre 'ZVNEW'; DZVN = ZVNEW '-' ZVAUX ; ZVARN = ZVARN '+' DZVN ; * mess 'Erreur retenue', ZERR4; * stockage de l erreur normalisee dans la liste des erreurs * cas particulier TWL2 Sellier theorie WL2 'SI' ('EGA' LEMO 'TWL2') ; ZVARN = ZVARN '+' ZSMAN '-' ZSMAX ; 'FINS' ; * ---- cas particulier SER Sellier renfort avec gradient * calcul GRADIENT des contraintes sur les renforts (SERi) * stockage dans la VARI HELM1(&BH,2) du gradient selon l absc curviligne 'SI' ('EGA' MOPRE 'SHR'); 'SI' ( 'EGA' NDIM1 3 ); * valeur SERi du pas precedent dans Helm(&bh,1) * valeur SERi du pas actuel * recuperation vecteur direction de projection * produit scalaire ZVAUX3=(NX11*GRA11)+(NX22*GRA22)+(NX33*GRA33); * stockage du gradient de epse dans YU2j * trac zvaux6 MOD_HELM titre '3D H2'; DZVN = ZVAUX6 '-' ZVAUX7; ZVARN = ZVARN '+' DZVN ; 'SINON'; * valeur SERi du pas precedent dans Helm(&bh,1) * recuperation vecteur direction de projection * produit scalaire ZVAUX3=(NX11*GRA11)+(NX22*GRA22); * stockage du gradient de epse dans YU2j * trac zvaux6 MOD_HELM titr '2DH2'; DZVN = ZVAUX6 '-' ZVAUX7; ZVARN = ZVARN '+' DZVN ; 'FINSI'; 'FINSI'; * ------fin du calcul des gradients et Hessiens-------- 'FIN' BH ; * ------fin boucle sur les formulations de Helmholtz --- 'FINS' ; * On rend ZVARN // si besoin : 'SI' (PARALLEL 'ET' ('NON' PARTLOCA)) ;