* G_CALCUL PROCEDUR JB251061 22/09/05 21:15:01 11444 * ============================================================================= * PROCEDURE DE CALCUL DES INTEGRALES NECESSAIRES A G_THETA * -------------------------------------------------------- * * DESCRIPTION : EFFECTUE LE CALCUL DE L'INTEGRALE SPECIFIEE DANS L'INDICE * 'OBJECTIF' DE LA TABLE SUPTAB. * ON BOUCLE SUR LES PAS DE CALCUL, SUR LES INTEGRALES A CALCULER * ET ENFIN SUR LES NOEUDS DU FRONT DE FISSURE. * ============================================================================= * QUELQUES OBJETS UTILES MOD_MEC_R = OBJUTI.'MOD_MEC_R'; SI BOOL.'DECOUPLAGE'; CH_AUX = SUPTAB.'CH_AUX'; FINSI; * CHAMPS NULS SI (EGA GMODE 'PLANGENE'); 3 'FZ' 0. 'MX' 0. 'MY' 0. 3 'UZ' 0. 'RX' 0. 'RY' 0. FINSI; * EN CAS D'UN MODELE DE PRESSION SI BOOL.'MODE_PRES'; FINSI; FINSI; SI (BOOL.'COQ' ET BOOL.'EL_QUA'); MAT2 = MAT_INST; SINON; FINSI; SI BOOL.'COQ'; SINON; FINSI; NBG = OBJUTI.'NBG'; * FRONT DE FISSURE EN POI1 FF = SUPTAB.'FRONT_FISSURE'; FINSI; FINSI; * QUELQUES MOTS POUR SIMPLIFIER L'ECRITURE NEXTR = GDIME; SI BOOL.'COQ'; NEXTR = 6; FINSI; SI (EGA GMODE 'AXIS'); FINSI; *BP: PETIT AJOUT DE NOM DE COMPOSANTE DE CHPOINT SI (BOOL.'J' OU BOOL.'J_DYNA') ; MOCOMP = 'J' ; FINSI; SI (BOOL.'C*' OU BOOL.'C*H') ; MOCOMP = 'C*' ; FINSI; SI BOOL.'DJ/DA' ; MOCOMP = 'DJDA' ; FINSI; * VALEURS INITIALES S10 = 0. ; S20 = 0. ; S30 = 0. ; S40 = 0. ; S41 = 0. ; S50 = 0.; S60 = 0. ; S70 = 0. ; S80 = 0. ; S90 = 0. ; S100 = 0.; S110 = 0. ; S120 = 0. ; S130 = 0. ; S140 = 0. ; S141 = 0.; S150 = 0. ; S160 = 0. ; S170 = 0. ; S180 = 0.; S190 = 0. ; S200 = 0. ; S210 = 0. ; S220 = 0.; * DETERMINATION DE LA TAILLE DE LA 1ERE COLONNE A AFFICHER SI ((EGA GDIME 3) ET (NON BOOL.'COQ')); FIN INO; SINON; FINSI; * TITRES A AFFICHER SELON LE PROBLEME TRAITER * ******************************************* * TITRE GENERAL NMESS = OBJUTI.'NMESS'; SI (BOOL.'J_DYNA'); SINON; SI (BOOL.'DJ/DA'); FINSI; FINSI; FINSI; SI BOOL.'THER'; SINON; FINSI; SI BOOL.'DECOUPLAGE'; SINON; SI (BOOL.'J' OU BOOL.'DJ/DA'); SINON; SI BOOL.'J_DYNA'; SINON; FINSI; FINSI; SI (BOOL.'C*H'); FINSI; SINON; FINSI; FINSI; * EN-TETE DES RESULTATS NCTITR = 0; NPMAX = NBG + NBDEP; FINSI; SI BOOL.'DECOUPLAGE'; NDDEC = NCTITR + 5; FINSI; SI BOOL.'COQ'; NDCOQ = NCTITR + 7; FINSI; SI ((EGA GDIME 3) ET (NON BOOL.'COQ')); SINON; FINSI; NDFRO = 6; FIN INO; NDFRO = NDFRO + NCTITR + 1; FINSI; NFLOT = 11 + 2; NDTHE = NDMEC + NFLOT; NDVOL = NDTHE + NFLOT; NDINT = NDVOL + NFLOT; SAUT 'LIGNE'; FIN ICHAR; SAUT 'LIGNE'; NDCHA = (NDINT / 2) + (45 / 2); MESS; MESS CTITR; MESS; *********************************************** *********************************************** ********* BOUCLE SUR LE PAS DE CALCUL ********* *********************************************** *********************************************** FINSI; REPE BOUCEXT NBDEP; IABC = NBG + &BOUCEXT; * DECLARATION DES TABLES STOCKANT LES RESULTATS POUR LE PAS DE TEMPS IABC SI (EGA GDIME 3); SI BOOL.'DECOUPLAGE'; REPE IRUPT GDIME; FIN IRUPT; SINON; FINSI; FINSI; FINSI; *************************************************** ** DEPLACEMENTS,CONTRAINTES ... A L INSTANT INST ** *************************************************** * SOLUTION_PASAPAS SI BOOL.'PERSO1'; * CAS PERSO1 INST = ESTIM . 'TEMPS'; SI (BOOL.'GRANDS_DEP' ET (NON BOOL.'ROT_RIG')); FINSI; SI BOOL.'ROT_RIG'; MESS 'ERREUR : Le deplacement du a une rotation'; MESS ' rigidifiante au pas ' IABC ' n est pas' MESS ' donne'; FINSI; DEPINT = DEPINT - FINSI; SINON; FINSI; SI BOOL.'C*'; SI (EGA IABC 0); DELTAT = INST + 1.E+30; DEPINT = ESTIM . 'DEPLACEMENTS'; VITDFI = ESTIM . 'DEFORMATIONS_INELASTIQUES'; SIG1 = SIGF * 1.; FINSI; SI (IABC '>' 0); DELTAT= INST - (ESTIM . 'TEMPS'); DEPINT= (ESTIM.'DEPLACEMENTS') - (ESTIM.'DEPLACEMENTS'); VITDFI= (ESTIM.'DEFORMATIONS_INELASTIQUES') - (ESTIM.'DEFORMATIONS_INELASTIQUES'); FINSI; FINSI; SI BOOL.'J_DYNA'; FINSI; SINON; * CAS OU ON APPELLE G_THETA APRES PASAPAS INST = SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'.IABC; DEPINT = (SUPTAB.'SOLUTION_PASAPAS'.'DEPLACEMENTS'.IABC) REDU ELTETA; SIGF = (SUPTAB.'SOLUTION_PASAPAS'.'CONTRAINTES'.IABC) REDU MOD_MEC_R; SI (BOOL.'GRANDS_DEP' ET (NON BOOL.'ROT_RIG')); FINSI; SI BOOL.'ROT_RIG'; MESS 'ERREUR : Le deplacement du a une rotation'; MESS ' rigidifiante au pas ' IABC ' n est pas' MESS ' donne'; FINSI; DEPINT = DEPINT - FINSI; VARF = (SUPTAB.'SOLUTION_PASAPAS'.'VARIABLES_INTERNES'.IABC) REDU MOD_MEC_R; SINON; FINSI; SI BOOL.'C*'; SI (EGA IABC 0); DELTAT = INST + 1.E+30; DEPINT = SUPTAB.'SOLUTION_PASAPAS'.'DEPLACEMENTS'.IABC; VITDFI = SUPTAB.'SOLUTION_PASAPAS'. 'DEFORMATIONS_INELASTIQUES'.IABC; SIG1 = SIGF * 1.; FINSI; SI (IABC > 0); DELTAT = INST - (SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'. (IABC - 1)); DEPINT = (SUPTAB.'SOLUTION_PASAPAS'.'DEPLACEMENTS'.IABC) - (SUPTAB.'SOLUTION_PASAPAS'.'DEPLACEMENTS'.(IABC - 1)); VITDFI = (SUPTAB.'SOLUTION_PASAPAS'. 'DEFORMATIONS_INELASTIQUES'.IABC) - (SUPTAB.'SOLUTION_PASAPAS'. 'DEFORMATIONS_INELASTIQUES'.(IABC - 1)); FINSI; FINSI; SI BOOL.'J_DYNA'; VITF = (SUPTAB.'SOLUTION_PASAPAS'.'VITESSES'.IABC) REDU ELTETA; ACCF = (SUPTAB.'SOLUTION_PASAPAS'.'ACCELERATIONS'.IABC) REDU ELTETA; FINSI; FINSI; SINON; * SOLUTION_RESO FINSI; * ON CHANGE LE DEPLACEMENT DEPINT EN MCHAML AU NOEUD **************************************************** * MODIFICATION DES CHAMPS THETA ET PI SI GRANDE ROT **************************************************** OTER SUPTAB 'TAB_THETA'; SI BOOL.'DJ/DA'; OTER SUPTAB 'TAB_PI'; FINSI; CH_THETA SUPTAB OBJUTI BOOL; FINSI; ********************************************************************** * CALCUL DES GRADIENTS ET DE LA DIVERGENCE DES CHAMPS THETA ********************************************************************** TTHETA = SUPTAB.'TAB_THETA'; REPE BCNOEU NBOU; IND0 = INDTH.&BCNOEU; CHTHETA = TTHETA.IND0 + DEP000; SUPTAB.'GRTHETA'.IND0 = GRTHETA; SUPTAB.'DIVTHETA'.IND0 = DIVTHETA; FIN BCNOEU; FINSI; *************************************************** ********** TEMPERATURES A L INSTANT INST ********** *************************************************** SI BOOL.'THER'; SI BOOL.'THER_DECO'; SINON; * bp, 2014-11-13 : ajout distinction cas BOOL.'THER' et BOOL.'THER_DECO' TEPINT = REDU ELTETA SUPTAB.'SOLUTION_PASAPAS'.'TEMPERATURES'. IABC; FINSI; SINON; TEPINT = SUPTAB.'TEMPERATURES'; FINSI; FINSI; *************************************************** ********** DEF IMPOSEE A L INSTANT INST ********** *************************************************** SI BOOL.'DEF_IMP'; SINON; DEFI = SUPTAB.'DEFORMATIONS_IMPOSEES'; FINSI; FINSI; *************************************************** ********** CONTACT FROTTANT ********** *************************************************** * ...TODO SINON; * A PRIORI DEPLACEMENT_FISSURE PAS TRES UTILE ... WDEP = SUPTAB . 'DEPLACEMENT_FISSURE'; SINON; FINSI; FINSI; SIGCON = SUPTAB . 'PRESSION_FISSURE'; SINON; FINSI; * PEUT ETRE FAIRE UN TEST SUR SIGCON ... FINSI; FINSI; ***************************************************** * CONTRAINTE RECALCULEE SI BOOL.'THER' = VRAI et NON BOOL.'PASAPAS' * ***************************************************** FINSI; *************************************************** ************ MATERIAU A L INSTANT INST ************ *************************************************** SINON; MAT_INST = MAT_MEC; * bp: ne devrait on pas ecrire ci dessous..? inutile? * MAT_INST = REDU MOD_MEC_R MAT_MEC FINSI; OBJUTI.'MAT_INST' = MAT_INST; *************************************************** ********* RIGIDITE TOTALE A L INSTANT INST ******** *************************************************** * bp: utilite? AUX_MECA = FAUX; SI BOOL.'DECOUPLAGE'; AUX_MECA = EGA SUPTAB.'METH_AUX' 'MECA'; FINSI; SI (BOOL.'DJ/DA' OU AUX_MECA); SI (BOOL.'GRADPROP' ET BOOL.'THER'); (SUPTAB.'SOLUTION_PASAPAS'.'CARACTERISTIQUES'); SINON; M1 = SUPTAB.'SOLUTION_PASAPAS'.'CARACTERISTIQUES'; FINSI; SINON; FINSI; FINSI; ************************************ *** RIGTOT + BLOCAGES MECANIQUES *** ************************************ SI BOOL.'DJ/DA'; RIGTOT = RIGTOT ET (SUPTAB.'SOLUTION_PASAPAS'.'BLOCAGES_MECANIQUES'); * rem bp: on suppose ceux ci identique a ceux de wtab... SINON; RIGTOT = RIGTOT ET SUPTAB.'BLOCAGES_MECANIQUES'; FINSI; FINSI; OBJUTI.'RIGTOT' = RIGTOT; ******************************************************************** ****** CHARGEMENT MECANIQUE (FORCES NODALES) A L INSTANT INST ****** ******************************************************************** PREINT = FOR000; *** SOLUTION_PASAPAS ****************************** FINSI; SI BOOL.'MODE_PRES'; * Dans le cas ou il y a un modele de pression * --> on isole la partie du MMODEL de pression qui est appliquee * sur la fissure SI BOOL.'PRES_FISS'; FINSI; * --> on calcule les forces nodales equivalentes aux pressions * appliquees hors de la fissure FINSI; FINSI; SINON; *** SOLUTION_RESO ********************************** FINSI; SI BOOL.'MODE_PRES'; * Dans le cas ou il y a un modele de pression * --> on isole la partie du MMODEL de pression qui est appliquee sur * la fissure SI BOOL.'PRES_FISS'; FINSI; * --> on calcule les forces nodales equivalentes aux pressions * appliquees hors de la fissure FINSI; FINSI; FINSI; ********************************************************************** ******** POUR LES PRESSIONS APPLIQUEES SUR LA FISSURE * on calcule le MCHAML vectoriel de [pression * normale] ********************************************************************** SI BOOL.'PRES_FISS'; * champ de normale unitaire au maillage de la fissure NF = NF2 / (XX**0.5); * champ de pression normale au maillage de la fissure PNF = PF * NF LCP LCOMP LCOMP; FINSI; **************************************************** ***** APPEL A G_AUX POUR LES CHAMPS AUXILIAIRES **** **************************************************** * G_AUX SI AUX_MECA; FINSI; SUPTAB.'CH_AUX' = CH_AUX; MESS; FINSI; *debut du cas contact frottant BOOL.'FROT' (btrolle 19/02/2013) OBJCON = OBJUTI.'OBJCON'; OBJCON2 = OBJUTI.'OBJCON2'; FINSI; *fin du cas contact frottant BOOL.'FROT' (btrolle 19/02/2013) *********************************************** *********************************************** **** BOUCLE SUR LES INTEGRALES A CALCULER ***** *********************************************** *********************************************** * NBMIXT = nbre d integrale a calculer (=1 si J, =2 si K1 K2, =3 * si K1 K2 K3) NBMIXT = 1; SI BOOL.'DECOUPLAGE'; C_MATE = OBJUTI.'C_MATE'; MU_1 = OBJUTI.'MU_1'; NBMIXT = 2; SI (EGA GDIME 3); NBMIXT = 3; FINSI; FINSI; REPE BOUCMIX NBMIXT; *|=====================================================================| *|======= I. OBJETS NECESSAIRES =====================================| SI BOOL.'DECOUPLAGE'; FINSI; **************************************************** **************** EN CAS DE DECOUPLAGE ************** **************************************************** SI BOOL.'DECOUPLAGE'; MOTMIX = CH_AUX.&BOUCMIX.'MOTMIX'; MOTMIA = CH_AUX.&BOUCMIX.'MOTMIA'; A_PREI = CH_AUX.&BOUCMIX.'A_PREI'; A_DEPI = CH_AUX.&BOUCMIX.'A_DEPI'; A_SIGF = CH_AUX.&BOUCMIX.'A_SIGF'; A_DEPGR = CH_AUX.&BOUCMIX.'A_DEPGR'; *debut du cas contact frottant BOOL.'FROT' (btrolle 19/02/2013) B_DEPGR = CH_AUX.&BOUCMIX.'B_DEPGR'; FINSI; *fin du cas contact frottant BOOL.'FROT' (btrolle 19/02/2013) FINSI; **************************************************** ******* EN CAS DE CALCUL EN VISCO_PLASTICITE ******* **************************************************** * INITIALISATION DE FACT1 FACT1 = 1.; SI (BOOL.'C*' OU BOOL.'C*H'); SI (EGA ITYPEF 1); FINSI; FINSI; SI BOOL.'C*'; SI (EGA ITYPEF 1); * DENSITE D'ENERGIE POUR LES FLUAGES DONT ON A UNE * EXPRESSION EXPLICITE DE L'INTEGRATION SUR LE TEMPS ENERM1 = (CHAF2*((CHAF2 + 1.)**(-1.)))*CHAF1*(VMI1**COE1); ENERM2 = (CHAF4*((CHAF4 + 1.)**(-1.)))*CHAF3*(VMI1**COE2); ENERM3 = (CHAF6*((CHAF6 + 1.)**(-1.)))*CHAF5*(VMI1**COE3); ENERM = ENERM1 + ENERM2 + ENERM3; SINON; SI ((EGA INST 0. 1.E-10) ET ('<'(COE2 - 1) 0.)); V1 = 0.; SINON; V1 = INST**(COE2 - 1); FINSI; ENERM = (CHAF2*((CHAF2 + 1.)**(-1.)))* CHAF1*(VMI1**COE1)*CHAF3*V1; FINSI; SI (BOOL.'COQ' ET BOOL.'EL_LIN') ; ENERM = ENERM*MOD_MEC_R EPAICH ; FINSI; FINSI; SI (EGA ITYPEF 2); * ON N'A PAS UNE EXPRESSION EXPLICITE DE * L'INTEGRATION DU FLUAGE SUR LE TEMPS SIGMOY = 0.5*(SIG1 + SIGF); SI ((EGA IABC 0) ET (NON BOOL.'REPRI')); SINON; FINSI; SIG11 = SIG1 ; SIG1 = SIGF ; VDI1 = VITDFI; FINSI; FINSI; * FACTEUR AFFECTE A J SI ON CALCULE C*(H) SI BOOL.'C*H'; N1 = COE2 / COE3; SI (EGA INST 0. 1.E-10) ; INST = 1.E-10 ; FINSI; V1 = 0.; FIN BF2; FIN BF1; SI (EGA FACT1 0. 1.E-10) ; FACT1 = 1.E+30 ; FINSI; FACT1 = ((V1 ** N1) / FACT1) ** COE3; FINSI; FINSI; ******************************************************** * SI LA COURBE DE TRACTION DEPEND DE LA TEMPERATURE ON * * CALCULE LA VARIATION DE CONTRAINTES DE VON-MISES LORS* * D'UNE AUGMENTATION (DETATE) DE LA TEMPERATURE A INST * ******************************************************** *bp 11/08/2011 : on calcule plutot la variation de la limite d elasticite * ce qui est moins faux si presence de decharge... * => on construit VM1 comm VM2 ! DETATE = 1.; EPS1 = MSQ1 * EPS1; * rem : on pourrait utiliser BORN dans le futur * VMI1 = (CHAN ('VMIS' MOD_MEC_R SIGF MAT_INST) * 'TYPE' 'SCALAIRE')*MSQ1; * DETAVM = VMI1 * 0.; DETAVM = 0.; REPETER BCMOD2 NBOBJ; * VM1 = REDU VMI1 MODI; SI (EGA MODPLA.&BCMOD2 1); 'ECRO' TABTRA.&BCMOD2) TEPINT; EPS2 'STRESSES' 'SCALAIRE'; 'ECRO' TABTRA.&BCMOD2) TEP1; * VM2 = 'VARI' 'NUAG' MODI (EXCO 'ECRO' MA2 SIGM) * (REDU EPS1 MODI) 'STRESSES' 'SCALAIRE'; * VM2 = ((EXCO VM2 SIGM 'SCAL')CHAN 'TYPE' 'SCALAIRE') * * (REDU MSQ1 MODI); EPS2 'STRESSES' 'SCALAIRE'; FINSI; SI ((EGA MODPLA.&BCMOD2 2) OU (EGA MODPLA.&BCMOD2 3)); * EPS2 = REDU EPS1 MODI; SI (EGA MODPLA.&BCMOD2 2); FINSI; * VM2 = VM2 * (REDU MSQ1 MODI); FINSI; DETAVM = ((VM2 - VM1) / DETATE); SINON; DETAVM = DETAVM + ((VM2 - VM1) / DETATE); FINSI; FINSI; FIN BCMOD2; FINSI; ******************************************************* **** ENERGIE DE DEFORMATION ELASTIQUE ET PLASTIQUE **** ******************************************************* *** *** DENSITE D'ENERGIE EN ELASTO OU THERMO-ELASTO-PLASTICITE ET *** DENSITE D'ENERGIE LIEE A LA VARIATION DE COURBE DE TRACTION *** SI (NON BOOL.'C*'); * SI (BOOL.'PASAPAS' ET (NON BOOL.'DECOUPLAGE')); SI (EGA IABC 0); SI (BOOL.'COQ' ET BOOL.'EL_LIN') ; VMI1 = VMI1*MOD_MEC_R EPAICH ; FINSI; FINSI; SINON; SI (BOOL.'COQ' ET BOOL.'EL_LIN') ; VMI1 = VMI1*MOD_MEC_R EPAICH ; FINSI; WVMIS = WVMIS + ((0.5*(DETAV1 + DETAVM))* FINSI; FINSI; ENERM = WELAS + WPLAS; SIG11 = SIG1 ; SIG1 = SIGF*1. ; VAR11 = VAR1 ; VAR1 = VARF*1.; DETAV1 = DETAVM; FINSI; SINON; ENERM = WELAS; FINSI; FINSI; *|=====================================================================| *|======= II. CRITERES DE DECHARGES =================================| ************************************************************ ****** CRITERE GLOBAL 1 DE DECHARGE DES CONTRAINTES ******** ************************************************************ *** *** EVALUATION EN ELASTO-PLASTICITE OU THERMO-ELASTO-PLASTICITE *** PAR COMPARAISON AVEC UN CALCUL EN ELASTICITE NON-LINEAIRE *** EN PRENANT POUR REFERENCE LES CONTRAINTES DE VON-MISES DU *** CALCUL RENDU PAR PASAPAS ET LES CONTRAINTES EQUIVALENTES LUES *** SUR LES COURBES DE TRACTION *** OUICRIT = FAUX; OUICRIT = SUPTAB.'CALCUL_CRITERE'; FINSI; * list modpla; FINSI; XCRIT = 1.; EXISCRIT = FAUX; * VMIS1 = contraintes de Von-Mises en elastoplastique TYPE 'SCALAIRE'; SI (EGA MODPLA.&BCMOD0 1); * modele est plastique isotrope. on utilise la * courbe de traction TABTRA.&BCMOD0 EXISCRIT = VRAI; 'ECRO' TABTRA.&BCMOD0) TEPINT; EPSTOT = EPELAST + EPSET; * VMIS2 = contraintes deduites des deformations * equivalentes par la courbe de traction EPSTOT 'STRESSES' 'SCALAIRE'; FINSI; SI (EGA MODPLA.&BCMOD0 2); * modele est plastique cinematique, cas ou * SIGY et H sont des evolutions EXISCRIT = VRAI; FINSI; SI (EGA MODPLA.&BCMOD0 3); * modele est plastique parfait EXISCRIT = VRAI; FINSI; SINON; *** modele elastoplastique et modele plastique cinematique *** avec caracteristiques independantes de la temperature * extraction de la courbe de traction du modele MODI EXISCRIT = VRAI; EPSTOT = EPELAST + EPSET; * VMIS2 = contraintes deduites des deformations * equivalentes par la courbe de traction * ? VMIS2 = VMIS2 CHAN TYPE 'SCALAIRE'; FINSI; * modele est plastique cinematique traite ici separement du cas * ou SIGY et H sont des evolutions c'est inutile car MAT_INST est * deja instantie en temperature mais il faudra regler le pb de * MODPLA dans son ensemble avec ZONE et une identification * automatique de chacun des sous-modes EXISCRIT = VRAI; FINSI; * modele est plastique parfait EXISCRIT = VRAI; FINSI; FINSI; SI EXISCRIT; VMIS1P = VMIS1 * MSQ; VMIS2P = VMIS2 * MSQ; SI ('NEG' VMIS1S 0.); XCRIT = XCRIT + ((VMIS2S/VMIS1S) - 1.); FINSI; FINSI; FIN BCMOD0; SUPTAB.'CRIT_DECHA_GLOBAL1'.IABC = XCRIT; FINSI; ********************************************************* ****** CRITERE LOCAL 1 DE DECHARGE DES CONTRAINTES ****** ********************************************************* FINSI; EXISCRIT = FAUX; * initialisation INST1 = SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'.(IABC-1); * TEPIN1 : temperature absolue au pas precedent * modele plastique isotrope SI (EGA MODPLA.&BCMOD0 1); 'ECRO' TABTRA.&BCMOD0) TEPINT; * VMISF : contrainte de Von Mises au pas courant 'ECRO' TABTRA.&BCMOD0) TEPIN1; * VMIS1 : contrainte de Von Mises au pas precedent FINSI; * modele plastique cinematique ou plastique parfait SI ((EGA MODPLA.&BCMOD0 2) OU (EGA MODPLA.&BCMOD0 3)); FINSI; * pour ces 3 modeles SI ((EGA MODPLA.&BCMOD0 1) OU (EGA MODPLA.&BCMOD0 2) OU (EGA MODPLA.&BCMOD0 3)); EXISCRIT = VRAI; VMISFP = VMISF * MSQ; VMIS1P = VMIS1 * MSQ; * VMISFP : contrainte de Von Mises aux points d'integration plastifies ** au pas courant * VMIS1P : contrainte de Von Mises aux points d'integration plastifies ** au pas precedent * on considere la plastification lorsque EPSET > 1E-6 VMISF1P = VMISFP - VMIS1P; * VMISFP : difference entre la contrainte equivalente au temps de * calcul t et celle au temps precedent (t-1) VMITOTF = VMITOTF + VMISFP; VMITOT1 = VMITOT1 + VMISF1P; FINSI; SINON; * modele elastoplastique avec ** caracteristiques independantes de la temperature TYPE 'SCALAIRE'; TYPE 'SCALAIRE'; FINSI; * modele plastique cinematique et plastique parfait avec ** caracteristiques independantes de la temperature TYPE 'SCALAIRE'; TYPE 'SCALAIRE'; FINSI; EXISCRIT = VRAI; VMISFP = VMISF * MSQ; VMIS1P = VMIS1 * MSQ; VMISF1P = VMISFP - VMIS1P; VMITOTF = VMITOTF + VMISFP; VMITOT1 = VMITOT1 + VMISF1P; FINSI; FINSI; FIN BCMOD0; SI EXISCRIT; * MSQO = (CHAN 'CHPO' MOD_MEC_R (VMITOT1*((VMITOTF)**(-1)))) * 'MASQ' 'INFERIEUR' 0.; * SUPTAB.'CRIT_DECHA_LOCAL1'.IABC = * (CHAN 'CHPO' MOD_MEC_R ((VMITOT1*((VMITOTF)**(-1)))) * MSQO); * #MC 04/11/98 : pour le MASQ, un produit suffit (pas de risque de /0) * #MC 04/11/98 : ou il y a des 0, on met 1.E-10 SUPTAB.'CRIT_DECHA_LOCAL1'.IABC = )**(-1)))) * MSQO); *lorsque la valeur du critere est positive, c'est-a-dire qu'il n'y a **pas de decharge, celle-ci est ramenee a zero FINSI; FINSI; ************************************************************ ****** CRITERE GLOBAL 2 DE DECHARGE DES CONTRAINTES ******** ************************************************************ FINSI; SI BOOL.'PERSO1'; EPPLASF = ESTIM . 'DEFORMATIONS_INELASTIQUES'; SINON; EPPLASF = SUPTAB.'SOLUTION_PASAPAS'.'DEFORMATIONS_INELASTIQUES' .IABC; FINSI; XCRIT = 1.; PRO1 = 0.; PRO2 = 0.; * initialisation * glob2 n'est pas calcule si reprise car on ne connait pas i-1 * VMISF : contrainte de Von Mises au pas courant VMISF = VMISF + (MQSIG * ((SUPTAB.VMISMAX) - VMISF)); FINSI; SUPTAB.VMISMAX = VMISF; EXISCRIT = FAUX; * EPPLASI = (SUPTAB.'SOLUTION_PASAPAS'. * 'DEFORMATIONS_INELASTIQUES'.IABC) REDU MODI; SI ((EGA MODPLA.&BCMOD0 1) OU (EGA MODPLA.&BCMOD0 2) OU (EGA MODPLA.&BCMOD0 3)); * modele plastique isotrope * modele plastique cinematique * modele plastique parfait * PRO1 = PRO1 + (INTG MODI (ENER MODI SIGFI EPPLASI)); SI (PRO2 'NEG' 0.D0); XCRIT = XCRIT + (1. - (PRO1/PRO2)); FINSI; FINSI; SINON; * modele elastoplastique et * modele plastique cinematique et plastique parfait ** avec caracteristiques independantes de la temperature * PRO1 = PRO1 + (INTG MODI (ENER MODI SIGFI EPPLASI)); * NUMERA = NUMERA + (ENER MODI SIGFI EPPLASI); DENOMI = DENOMI + (VMISFI * EPSEI); * petit traitement pour eviter de diviser par zero SI (AMP < 1.E-15); AMP = 1.E-15; FINSI; SI (MIN0 < (1.E-15 * AMP)); DENOMI = DENOMI + (1.E-15 * AMP); FINSI; CRITLOCA = CRITLOCA + (1. - (NUMERA / DENOMI)); FINSI; FINSI; FIN BCMOD0; * Le maxi sert a corriger le probleme castem sur les EPSE * (ENER MODI SIGFI EPPLASI) n est pas toujours egal a * (VMISFI * EPSEI) lorsque le chargement est proportionnel * A supprimer quand le pb sera resolu SUPTAB.'CRIT_DECHA_GLOBAL2_L'.IABC= CRITLOCA; * mess 'PRO1' , PRO1 , 'PRO2' , PRO2 , 'XCRIT' , XCRIT; SINON; MESS 'ATTENTION : LES DEFORMATIONS INELASTIQUES N''ONT'; MESS ' ETRE CALCULE.'; FINSI; FINSI; FINSI; ********************************************************* ****** CRITERE LOCAL 2 DE DECHARGE DES CONTRAINTES ****** ********************************************************* *initialisation EXISCRIT = FAUX; TYPE 'CARACTERISTIQUES' 'RIGIDITE'; EXISCRIT = VRAI; NOR1 = SIGFTRA**(1./2); DNOR = DSIGTRA**(1./2); DENO = NOR1 * DNOR; CRIT0 = CRIT0 + (INV1 * (DENO**(-1))); FINSI; FIN BCMOD0; SI EXISCRIT; FINSI; FINSI; *|=====================================================================| *|======= III. TRAITEMENT DES CHAMPS MECANIQUES ====================| **************************************************** * ON CHANGE LA TEMPERATURE EN MCHAML AU NOEUD =TEPINT **************************************************** SINON; FINSI; FINSI; SI (BOOL.'THER' ET (NON BOOL.'COQ')); EPT1 = EPT1 - OBJUTI.'EPTREF'; FINSI; *************************************************** ***** GRADIENT DU DEPLACEMENT DEPINT =GRADEP ****** *************************************************** SI (BOOL.'COQ' ET BOOL.'EL_LIN'); FINSI; *************************************************** ***** CAS DU CONTACT FROTTANT : ON CALCULE WSAUTGR = GRADIENT DU SAUT ****** *************************************************** * WDEP = deplacement de la fissure de composante UX ... AX ... * [grad(w)] = grad [w] car linearite (le saut [w] est donne par AX ...) * on utilise un modele liee a la geo de la fissure (cohesif ou contact) * SIGCON = sigma*n de composantes SMX... homogene a une contrainte * supporte sur la levre superieure FINSI; * Vsigcon = 'VECTEUR' sigcon objcon (1E-8) ('MOTS' 'SMX' 'SMY' 'SMZ'); * 'TRACER' vsigcon ('EXTRAIRE' objcon 'MAIL'); SI (EGA GDIME 3); FINSI; FINSI; *************************************************** ***** GRADIENT DE LA VITESSE VITF =GRAVIT ****** *************************************************** SI BOOL.'J_DYNA'; FINSI; *************************************************** ** GRADIENT TEMPERATURE Grad T pour tous elements =TEPEGR ** et Grad Grad T pour les elements massifs seuls (si DJ/DA) =DEPDTGR *************************************************** SI BOOL.'THER'; SI (NON BOOL.'COQ'); SI BOOL.'DJ/DA'; SI (GDIME EGA 3); FINSI; * DEPDTGR = GRAD MOD_MEC_R MAT_INST (CHAN 'CHPO' MOD_MEC_R DEPDT); FINSI; SINON; + CMD000); FINSI; FINSI; *************************************************** **GRADIENT DE PRESSIONS Grad P (si DJ/DA) =PRESGR *************************************************** SI BOOL.'DJ/DA'; FIN B1; FINSI; FINSI; *************************************************** * PROFIL DE L'ENERGIE DANS L'EPAISSEUR DE LA COQUE *************************************************** SI BOOL.'COQ'; MCOU1 = OBJUTI.'MODCOU'.&NBJ7; SI ('<' (E1 / EPAITO) 1.E-4); FINSI; FIN NBJ7; SI (EGA SOM1 0.); SOM1 = SOM1 + 1.E-10; FINSI; EVENR = EVENR * EPAITO / SOM1; FINSI; **************************************************** * RECUPERATION DU NOMBRE DE CALCULS A FAIRE **************************************************** TTHETA = SUPTAB.'TAB_THETA'; SI BOOL.'DJ/DA'; TPI = SUPTAB.'TAB_PI'; FINSI; ******************************************************** * DENSITE DE MATERIAU EN CAS DE CALCUL DE J DYNAMIQUE ** ******************************************************** SI BOOL.'J_DYNA'; FINSI; *|=====================================================================| *|======= BOUCLE SUR LES NOEUDS A AVANCER VIRTUELLEMENT =============> REPE BCNOEU NBOU; ************************************************** * RECUP DU CHAMP THETA ASSOCIE AU NOEUD PM DE NUMERO NUNOE * ET DU CHAMP PI (SI DJ/DA) ************************************************** IND0 = INDTH.&BCNOEU; CHTHETA = TTHETA.IND0 + DEP000; SI BOOL.'DJ/DA'; CHPI = TPI.(INDPI.&BCNOEU); FINSI; ************************************************** * GRADIENT, DIVERGENCE DU CHAMP THETA =TETAGR, DIVTETA * CHAMELEM THETA =(TETX,TETY,TETZ) ************************************************** TETAGR = SUPTAB.'GRTHETA'.IND0; DIVTETA = SUPTAB.'DIVTHETA'.IND0; * btrolle 2013 : ajout de TETAX, TETAY et TETAZ SI (GDIME EGA 3); FINSI; ************************************************** * GRADIENT, DIVERGENCE DU CHAMP PI =PIGR, DIVPI * CHAMELEM PI =(PIX,PIY,PIZ) * ... ************************************************** SI BOOL.'DJ/DA'; SI (EGA GDIME 3); FINSI; CHPI = CHPI + DEP000; PITAGR = (MOD_MEC_R PIGR * TETAGR); ADJ = (DIVPITA - (DIVPI * DIVTETA)); SI (GDIME EGA 3); FINSI; FINSI; ******************************************************** *********** SI LE MATERIAU N'EST PAS CONSTANT ********** * YOU1 = YOUGR * THETA = (grad E) * THETA * SIGPRIM = (grad D)*THETA*(EPSI elas) ou grad D =D(grad E,nu) =DMAT * ALF1 = ALFGR*THETA = (grad Alpha)*THETA * DMAT = D(1,nu,ALF1) ******************************************************** SI BOOL.'GRADPROP'; *** Le coefficient de Poisson est constant *** * NU1 = REDU (EXCO MAT_INST 'NU' 'SCAL') MOD_MEC_R; * NU1 = EXCO NU1 'SCAL' 'NU'; ****** Gradient de Module d'young ****** SI BOOL.'GRADYOUN'; YOUGR = GRA000 ; I = 0; I = I + 1; * YO1 = CHAN 'CHPO' TABMOD.I * (REDU (EXCO MAT_INST 'YOUN' 'SCAL') TABMOD.I); * YOUGR = YOUGR + (GRAD TABMOD.I MAT_INST ((NOMC MU1 YO1) * + (REDU DEP000 (EXTR YO1 'MAIL')))); FIN NBJ2; MOD_MEC_R TETX) + MOD_MEC_R TETY); SI (GDIME EGA 3); FINSI; FINSI; ****** Gradient de Coefficient de Dilatation ****** SI (BOOL.'THER' ET BOOL.'GRADALPH'); ALFGR = GRA000 ; I = 0; I = I + 1; *! AL1 = CHAN 'CHPO' TABMOD.I *! (REDU (EXCO MAT_INST 'ALPH' 'SCAL') TABMOD.I); *! ALFGR = ALFGR + (GRAD TABMOD.I MAT_INST ((NOMC MU1 AL1) *! + (REDU DEP000 (EXTR AL1 'MAIL')))); FIN NBJ3; MOD_MEC_R TETX) + MOD_MEC_R TETY); SI (GDIME EGA 3); MOD_MEC_R TETZ); FINSI; FINSI; FINSI; *********************************************** * TEMU = (Grad T)*THETA * TEMU1 = (Grad T)*PI * TEMU2 = (Grad (Grad T))*PI*THETA * EPSTU = D^-1 sigma^th * avec sigma^th = (-3lambda-2mu)*alpha*(grad T)*THETA*I * GR_EPTH_THET = (grad epsilon_th)*THETA *********************************************** SI BOOL.'THER'; SI (NON BOOL.'COQ'); 'TALP' 0. 'TREF' 0.; SINON; MOD_MEC_R TETX) + MOD_MEC_R TETY) + MOD_MEC_R TETZ); MOD_MEC_R TETX) + MOD_MEC_R TETY) + MOD_MEC_R TETZ); MOD_MEC_R TETX) + MOD_MEC_R TETY) + MOD_MEC_R TETZ); FINSI; SI BOOL.'DJ/DA'; SI (GDIME EGA 3); FINSI; SI BOOL.'COQ'; SINON; 'TEMPERATURES'; FINSI; MAT_INST; MOD_MEC_R PIX; MOD_MEC_R PIY; MOD_MEC_R PIX; MOD_MEC_R PIY; TEMU2 = (TETX*MOD_MEC_R (TXXPIX + TYXPIY)) + (TETY*MOD_MEC_R (TXYPIX + TYYPIY)); SI (GDIME EGA 3); TEMU2 = (TETX*MOD_MEC_R (TXXPIX + TYXPIY + TZXPIZ)) + (TETY*MOD_MEC_R (TXYPIX + TYYPIY + TZYPIZ)) + (TETZ*MOD_MEC_R (TXZPIX + TYZPIY + TZZPIZ)); FINSI; SI BOOL.'COQ'; SINON; FINSI; FINSI; FINSI; * modif sm *********************************************** * CUMDEFI = SIG_{IJ} * (deps_{ij} / dx_{,k}) * thet_{k} *********************************************** SI (BOOL.'DEF_IMP' ET (NON BOOL.'COQ')); REPETER BDEFI NIDEF; GREIJ = GRAD MOD_MEC_R MAT_INST CUMDEFI = CUMDEFI + DEIJX + DEIJY; SI (GDIME EGA 3); CUMDEFI = CUMDEFI + DEIJZ; FINSI; FIN BDEFI; FINSI; *********************************************** * PEMU1 = (Grad P)*TETA * PEMU2 = (Grad P)*(Grad PI)*TETA *********************************************** SI BOOL.'DJ/DA'; * NOUVELLE IMPLEMENTATION JB : * PEMU = VIDE 'CHPOINT'/'DISCRET'; * REPE ICOMP GDIME; * IDEB = 1 + (3*(&ICOMP - 1)); * IFIN = IDEB + GDIME - 1; * LMO1 = EXTR MGI (LECT IDEB PAS 1 IFIN); * LMO2 = EXTR MUI (LECT 1 PAS GDIME); * FI = PSCA PRESGR PI LMO1 LMO2; * FI = NOMC (EXTR MUI &ICOMP) FI; * PEMU = PEMU + (CHAN 'CHPO' MOD_MEC_R FI); * FIN ICOMP; * CETTE SECTION N'EST PAS TESTEE ACTUELLEMENT DANS LA BASE DE CAS-TESTS * IL FAUT DONC AJOUTER UN CAS-TEST AVANT DE POUVOIR REMPLACER CE QUI SUIT * PAR LES LIGNES COMMENTEES CI-DESSUS MOD_MEC_R PIX) + MOD_MEC_R PIY); MOD_MEC_R PIX) + MOD_MEC_R PIY); SI (EGA GDIME 3); F1 = F1 + MOD_MEC_R PIZ); F2 = F2 + MOD_MEC_R PIZ); MOD_MEC_R PIX) + MOD_MEC_R PIY) + MOD_MEC_R PIZ); FINSI; MOD_MEC_R TETX) + MOD_MEC_R TETY); MOD_MEC_R TETX) + MOD_MEC_R TETY); SI (EGA GDIME 3); F1 = F1 + MOD_MEC_R TETZ); F2 = F2 + MOD_MEC_R TETZ); MOD_MEC_R TETX) + MOD_MEC_R TETY) + MOD_MEC_R TETZ); FINSI; MOD_MEC_R TETX) + MOD_MEC_R TETY); MOD_MEC_R TETX) + MOD_MEC_R TETY); SI (EGA GDIME 3); F1 = F1 + MOD_MEC_R TETZ); F2 = F2 + MOD_MEC_R TETZ); MOD_MEC_R TETX) + MOD_MEC_R TETY) + MOD_MEC_R TETZ); FINSI; MOD_MEC_R F1) + MOD_MEC_R F2); MOD_MEC_R F1) + MOD_MEC_R F2); SI (EGA GDIME 3); F4 = F4 + MOD_MEC_R F3); F5 = F5 + MOD_MEC_R F3); MOD_MEC_R F1) + MOD_MEC_R F2) + MOD_MEC_R F3); FINSI; FINSI; FINSI; *********************************************** * DANS LE CAS DE CALCUL DE DJ/DA * U^aux=A_DEPI et SIGMA^aux=A_SIGF *********************************************** SI BOOL.'DJ/DA'; GRAD11 = (MOD_MEC_R GRADEP * PIGR); SI (GDIME EGA 3); EPSIA1 = MANU 'CHML' MOD_MEC_R EP1 EPXX1 EP2 EPYY1 EP3 EPZZ1 EP4 GAXY1 EP5 GAXZ1 EP6 GAYZ1 TYPE 'DEFORMATIONS' 'STRESSES'; SINON; EPSIA1 = MANU 'CHML' MOD_MEC_R EP1 EPXX1 EP2 EPYY1 EP3 EPZZ1 EP4 GAXY1 TYPE 'DEFORMATIONS' 'STRESSES'; FINSI; SI BOOL.'THER' ; EPSIA1 = EPSIA1 + EPSTU1 ; FINSI; FINSI; A_PREI = FOR000; SI (BOOL.'COQ' ET BOOL.'EL_LIN'); SINON; FINSI; FINSI; *********************************************** * DEP0 = (grad U)*THETA *********************************************** LO1 = ((BOOL.'J') OU BOOL.'C*' OU BOOL.'C*H' OU BOOL.'DJ/DA' OU BOOL.'J_DYNA' ) LO2 = FAUX; SI BOOL.'DECOUPLAGE'; FINSI; SI (LO1 OU LO2); * mess 'gcalcul : CHAN CHPO grad(U)'; *bp : on se ramene a un chpoint pour faire XTY avec preint ensuite * on moyenne avce CHAN CHPO, pourrait on l'eviter? SI (GDIME EGA 3); DEPLX = DEPLX + DEPLY = DEPLY + FINSI; SI (BOOL.'COQ' ET BOOL.'EL_LIN'); DEP0 = DEP0 + DEPFX + DEPFY; FINSI; FINSI; *************************************************************** * S'il y a une pression sur la fissure, on calcule * GRUTE = (grad U)*THETA sur les elements du modele de pression * et dans la zone d'etude (ou le champ TETA est definit) *************************************************************** SI BOOL.'PRES_FISS'; * on isole les elements de la fissure ou la pression s'applique et ou le * champ teta est definit * on interpole le gradient de deplacement (GRADEP) sur les elements du * modele de pression touchant la levre inferieure MODPFT = MODPFT ET MODPFI; GRADEPF = GRADEPF ET GRUI; FINSI; * on interpole le gradient de deplacement (GRADEP) sur les elements du * modele de pression touchant la levre superieure MODPFT = MODPFT ET MODPFS; GRADEPF = GRADEPF ET GRUS; FINSI; * on interpole le champ theta (THETA) sur les elements du modele de * pression ==> TETAF * on fait ensuite le produit (GRADEPF * TETAF) ==> GRUTE SI (EGA GDIME 2); FINSI; SI (EGA GDIME 3); FINSI; FINSI; *********************************************** * DEP1 = (grad A_DEPI)*THETA *********************************************** LO2 = FAUX; SI BOOL.'DECOUPLAGE'; FINSI; SI (LO1 OU LO2); * mess 'gcalcul : CHAN CHPO grad(U^aux)'; *bp : on se ramene a un chpoint pour faire XTY avec preint ensuite * on moyenne avec CHAN CHPO, pourrait on l'eviter? SI (GDIME EGA 3); DEPLX = DEPLX + DEPLY = DEPLY + FINSI; SI (BOOL.'COQ' ET BOOL.'EL_LIN'); DEP1 = DEP1 + DEPFX + DEPFY; FINSI; FINSI; *********************************************** * Termes supplemenetaire en Dynamique * VCARE = V*V * GRUTV = (grad U)*THETA*V (Attigui) * GRVTV = V^T * (grad V)*THETA (BP) * GRWTU = W^T * (grad U)*THETA (BP) *********************************************** SI BOOL.'J_DYNA'; * NOUVELLE IMPLEMENTATION JB : * CHML_TETA = CHAN 'CHAM' MOD_MEC_R TETA 'STRESSES'; * GRUTV = MANU 'CHML' MOD_MEC_R 'SCAL' 0. * 'TYPE' 'DEPLACEM' 'STRESSES'; * GRVTV = COPI GRUTV; * GRWTU = COPI GRUTV; * REPE ICOMP GDIME; * CALCUL DE DEPFI ET VITFI * IDEB = 1 + (3*(&ICOMP - 1)); * IFIN = IDEB + GDIME - 1; * LMO1 = EXTR MGI (LECT IDEB PAS 1 IFIN); * LMO2 = MUI; * DEPFI = PSCA GRADEP CHML_TETA LMO1 LMO2; * VITFI = PSCA GRAVIT CHML_TETA LMO1 LMO2; * CALCUL DE GRUTV * MO1 = EXTR MUI &ICOMP; * GRUTV = GRUTV + ((EXCO MO1 VITF 'SCAL') * DEPFI); * CALCUL DE GRVTV * GRVTV = GRVTV + ((EXCO MO1 VITF 'SCAL') * VITFI); * CALCUL DE GRWTU * GRWTU = GRWTU + ((EXCO MO1 ACCF 'SCAL') * DEPFI); * FIN ICOMP; * CETTE SECTION N'EST PAS ASSEZ TESTEE ACTUELLEMENT DANS LA BASE DE CAS-TESTS * IL FAUT DONC AJOUTER DES CAS-TESTS AVANT DE POUVOIR REMPLACER CE QUI SUIT * PAR LES LIGNES COMMENTEES CI-DESSUS SI (EGA GDIME 3); SINON; FINSI; FINSI; *|=====================================================================| *|======== IV. CALCUL DE J,C*,C*(h) ou J_DYNA =======================| *********************************************** * S10 = w*(div THETA) * S20 = sigma*(grad U)*(grad THETA) * S40 = Tr(sigma)*Alpha*(Grad T)*THETA * S50 = F*(grad U)*THETA * S60 = 0.5*(grad D)*THETA*(EPSI elas)*(EPSI elas) * S70 = Tr(sigma)*T*(Grad Alpha)*THETA * S100 = (W(b) - W(a))*THETA(x) * S110 = SIGF*n*[ (dU/dX)(b) - (dU/dX)(a) ]*THETA(x) * S120 = wvmis*(Grad T)*THETA * ou wvmis = SOME d(Von_mises)/d(Temperature) d(EPSE) *********************************************** * Termes supplementaires pour le J dynamique (Formulation de Attigui) * S130 = 0.5*RHO*(V^2)*(div THETA) * S140 = [Delta (RHO*V*(grad U)*THETA)]/(Delta Temps) * Termes supplementaires pour le J dynamique (Formulation de BP) * S130 = -0.5*RHO*(V^2)*(div THETA) * S140 = RHO*W*(grad U)*THETA * S141 = -RHO*V*(grad V)*THETA *********************************************** * Termes supplementaires pour le J en grands-deplacements * S150 = sigma*((grad U)t)*(grad U)*(grad THETA) *********************************************** *********************************************** * Termes supplementaires pour un chaergement en deformation imposee (modif sm) * S160 = sigma*(deps/dx)*(THETA) * = SIG_{IJ} * (deps^{imp}_{ij} / dx_{,k}) * thet_{k} *********************************************** * Termes supplementaires pour le contact frottant (modif BP,BT) * S111 = SIGF*n*[ (dU/dX)(b) - (dU/dX)(a) ]*THETA(x) *********************************************** SI ((BOOL.'J') OU BOOL.'C*' OU BOOL.'C*H' OU BOOL.'J_DYNA' ); GMCANI = 0.; GTERMI = 0.; GPRESS = 0.; GMCANI = GMCANI - S10; SI (BOOL.'COQ' ET BOOL.'EL_LIN'); (GRADEP*MOD_MEC_R TETAGR) (GRADEF*MOD_MEC_R TETAGR)); SINON; (GRADEP*MOD_MEC_R TETAGR)) MAT2; FINSI; GMCANI = GMCANI + S20; SI BOOL.'GRANDS_DEP'; SI (BOOL.'COQ' ET BOOL.'EL_LIN'); SINON; FINSI; GMCANI = GMCANI + S150; FINSI; SI (BOOL.'GRADPROP' ET BOOL.'GRADYOUN'); GMCANI = GMCANI - S60; FINSI; I = I + 1; LE1 = IND1.&NBJ4; L1 = LINTER.LE1; *** ((ENRMBL - ENRMAL) * OBJINT TETXLC) MAT2); *** 'STRESSES'; 'STRESSES'; 'STRESSES'; 'STRESSES'; 'STRESSES'; 'STRESSES'; AAAA1 = (SIXXL * OBJINT (UXXBL - UXXAL)) + (SIXYL * OBJINT (UYXBL - UYXAL)); FIN NBJ4; GMCANI = GMCANI - S100 + S110; FINSI; * Termes supplementaire en Dynamique ********** SI BOOL.'J_DYNA'; GMCANI = GMCANI - S130 + S140 - S141; FINSI; * Termes supplementaire pour le contact frottant (btrolle 19/02/2013) * teta SI (EGA GDIME 3); FINSI; * [grad(w)] = grad [w] = WSAUTGR * sigma*n SI (EGA GDIME 3); AAAAX = AAAAX AAAAY = AAAAY AAAA111 = (AAAAX * OBJCON2 TETXC) + (AAAAY * OBJCON2 TETYC) + (AAAAZ * OBJCON2 TETZC); SINON; AAAA111 = (AAAAX * OBJCON2 TETXC) + (AAAAY * OBJCON2 TETYC); FINSI; * AAAA111 = 'MANUEL' CHML objcon 'SCAL' 1 TYPE 'SCALAIRE' 'STRESSES'; GMCANI = GMCANI - S111; SINON; S111 = 0.0; FINSI; * fin du cas contact frottant (btrolle 19/02/2013) S50 = 0.; FINSI; SI BOOL.'PRES_FISS'; FIN B1; FINSI; GPRESS = GPRESS - S50; * Termes supplementaire pour un chargement en deformation imposee (SM) ********* S160 = 0.; SI (BOOL.'J' ET BOOL.'DEF_IMP'); GPRESS = GPRESS + S160; FINSI; SI BOOL.'THER'; SI BOOL.'COQ'; GTERMI = GTERMI + S40; SI BOOL.'GRADALPH'; GTERMI = GTERMI + S70; FINSI ; SINON; GTERMI = GTERMI + S41; FINSI; GTERMI = GTERMI - S120 ; FINSI; FINSI; * somme des termes *************************** GTOTA = (GMCANI + GTERMI + GPRESS)*FACT1; FINSI; *|=====================================================================| *|======= V. CALCUL DE DJ/DA ========================================| *********************************************** * S10 = -SIGF*(Grad U)*(Grad PI)*(Grad TETA) * S20 = -SIGF*(Grad U)*(Grad TETA)*(Grad PI) * S30 = SIGF*(Grad U)*(Grad TETA)*(Div PI) * S40 = SIGF*(Grad U)*(Grad PI)*(Div TETA) * S50 = ENEGIE*(ADJ TETA*PI) * S60 = ALPH*SIGF*(Grad (Grad T))*PI*TETA * S70 = ALPH*SIGF*((Grad T)*PI)*(Div TETA) * S80 = ALPH*SIGF*((Grad T)*TETA)*(Div PI) * S90 = SIG11*(Grad U)*(Grad TETA) * S100= SIGF*(Grad U1)*(Grad TETA) * S110= -SIGF*(Grad U1)*(Div TETA) * S120= ALPH*SIG11*((Grad T)*TETA) * S130= -PRESSION*(grad A_DEPI)*THETA * S140= -(Grad PRESSION)*(Grad PI)*TETA*U * S150= (Grad PRESSION)*TETA*U*(Div PI) * S160= -(Grad PRESSION)*PI*(Grad U)*TETA *********************************************** SI BOOL.'DJ/DA'; (MOD_MEC_R GRADEP * (MOD_MEC_R PIGR * TETAGR)))); (MOD_MEC_R GRADEP * (MOD_MEC_R TETAGR * PIGR)))); SIGF (MOD_MEC_R GRADEP * TETAGR)) * DIVPI); SIGF (MOD_MEC_R GRADEP * PIGR)) * DIVTETA); SI BOOL.'THER'; FINSI; FINSI; MOD_MEC_R SIGF A_DEPGR) * DIVTETA)); GMCANI = S10 + S20 + S30 + S40 + S50 + S90 + S100 + S110; GTERMI = S60 + S70 + S80 + S120; GPRESS = S150 - S130 - S140 - S160; GTOTA = GMCANI + GTERMI + GPRESS; FINSI; *|=====================================================================| *|====== VI. CALCUL DES FACTEURS D INTENSITE DES CONTRAINTES =======| *********************************************** * En pratique on les deduit de l integrale d interaction * M(u,u^aux) = 2/E^* (K1*K1^aux + K2*K2^aux) * On calcule donc : * J1 = 1/E^* K1^2 = M(u,u^aux_1)**2 / 4*J(u^aux_1,u^aux_1) * J2 = 1/E^* K2^2 = M(u,u^aux_2)**2 / 4*J(u^aux_2,u^aux_2) *********************************************** * l'integrale J(u^aux, u^aux) est la somme de : * S10 = 0.5*Tr(A_sigf*(grad A_depi))*(div THETA) * S20 = Tr[A_sigf*(grad A_depi)*(grad THETA)] * * S50 = A_prei*(grad A_depi)*THETA *********************************************** * l'integrale M(u, u^aux) et la somme de : * S60 = Tr[A_sigf*(Grad U)*(Grad THETA)] * S70 = Tr[SIGF*(Grad A_depi)*(Grad THETA)] * S80 = Tr[SIGF*(Grad A_depi)*(Div THETA)] * S90 = PREINT*(grad A_depi)*THETA * S100= A_prei*(grad U)*THETA * S110= ALPH*A_sigf*((Grad T)*TETA) *********************************************** * Termes supplementaires pour le contact frottant (modif BP,BT) * S111 = SIGF*n*[ (dU^aux/dX)(b) - (dU^aux/dX)(a) ]*TETA(x) *********************************************** SI BOOL.'DECOUPLAGE'; * si(BOOL.'XFEM'); * * REM: pour les elements XFEM on a contruit les champs aux. tq * * (K_I^aux, K_II^aux) = (1,0) et (0,1) * * => inutile de calculer J(u^aux,u^aux) (= C_MATE) *bp 07/07/2011 : cela est faux dans le cas d'une fissure courbee * COE_M2K = 2. * C_MATE; * SINON; S20 = INTG MOD_MEC_R FINSI; JETOILE = 'ABS' (S20 - S10 - S50); SI ((EGA &BOUCMIX 1) OU (EGA &BOUCMIX 2)); FACTK = 0.5 / ((C_MATE * JETOILE)**0.5); SINON; FACTK = 0.5 * ((2. * MU_1 / JETOILE)**0.5); FINSI; * FINSI; * integrale M(u, u^aux) *rem: on pourrait calculer M avec les 2 lignes ci dessous: * S81 = INTG MOD_MEC_R (MOD_MEC_R (WORK MOD_MEC_R A_SIGF GRADEP) * DIVTETA); * S80 = 0.5*(S80+S81); * Termes supplementaire pour le contact frottant (btrolle 19/02/2013) * on calcule WSAUTGR^aux = GRADIENT DU SAUT^aux * WDEP = deplacement de la fissure de composante UX ... AX ... * [grad(w)] = grad [w] car linearite (le saut [w] est donne par AX ...) * on utilise un modele liee a la geo de la fissure (cohesif ou contact) * mess 'I avec contact frottant'; * SIGCON = sigma*n de composantes SMX... homogene a une contrainte * supporte sur la levre superieure FINSI; SI (EGA GDIME 3); FINSI; * teta SI (EGA GDIME 3); FINSI; * [grad(w)] = grad [w] = WSAUTGR * sigma*n SI (EGA GDIME 3); AAAAX = AAAAX AAAAY = AAAAY AAAA111 = (AAAAX * OBJCON2 TETXC) + (AAAAY * OBJCON2 TETYC) + (AAAAZ * OBJCON2 TETZC); SINON; AAAA111 = (AAAAX * OBJCON2 TETXC) + (AAAAY * OBJCON2 TETYC); FINSI; SINON; S111 = 0.0; FINSI; * fin du cas contact frottant (btrolle 19/02/2013) FINSI; FINSI; SI BOOL.'THER'; FINSI; * facteurs d intensite des contraintes GMCANI = (S60 + S70 - S80 + S111) * FACTK; GTERMI = (S110) * FACTK; GPRESS = (-1.*(S100 + S90)) * FACTK; GTOTA = GMCANI + GTERMI + GPRESS; FINSI; *|=====================================================================| *|======= VII. STOCKAGE DES RESULTATS ET IMPRESSIONS ===============| SI ((NBOU > 1) ET (EGA &BCNOEU NBOU)); SINON; SI ((EGA GDIME 2) OU BOOL.'COQ'); SINON; NUNOE = &BCNOEU; SINON; FINSI; FINSI; FINSI; *************************************************** ** PROFILE DE G DANS L EPAISSEUR EN CAS DE COQUE ** *************************************************** SI BOOL.'COQ'; SUPTAB.'EPAISSEUR_RESULTATS'.IABC = EVENR * GTOTA; SINON; SUPTAB.'EPAISSEUR_RESULTATS' = EVENR * GTOTA; FINSI; FINSI; **************************** ** AFFICHAGE DES RESULTATS** **************************** FINSI; SI BOOL.'DECOUPLAGE'; FINSI; SI ((EGA GDIME 3) ET (NON BOOL.'COQ')); FINSI; SI BOOL.'COQ'; REPE IVAL 4; (GPRESS*XVAL)*NDVOL (GTOTA*XVAL)*NDINT; FIN IVAL; SINON; GTOTA*NDINT; FINSI; *************************************************** ** STOCKAGE DES RESULTATS DANS SUPTAB.'RESULTATS'** *************************************************** * CAS PASAPAS SI (EGA GDIME 2); SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX.IABC = GTOTA; SINON; TABRES.IABC = GTOTA; FINSI; FINSI; SI ((EGA GDIME 3) ET (NON BOOL.'COQ')); SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX.IABC.NUNOE = GTOTA; SINON; TABRES.IABC.NUNOE = GTOTA; FINSI; FINSI; SI BOOL.'COQ'; SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX = GTOTA; SINON; TABRES.IABC.'SUPERI' = GTOTA*V_SUPE; TABRES.IABC.'MEDIAN' = GTOTA*V_MOYE; TABRES.IABC.'INFERI' = GTOTA*V_INFE; TABRES.IABC.'GLOBAL' = GTOTA; FINSI; FINSI; * CAS RESOU SINON; SI (EGA GDIME 2); SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX = GTOTA; SINON; SUPTAB.'RESULTATS' = GTOTA; FINSI; FINSI; SI ((EGA GDIME 3) ET (NON BOOL.'COQ')); SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX.NUNOE = GTOTA; SINON; TABRES.NUNOE = GTOTA; FINSI; FINSI; SI BOOL.'COQ'; SI BOOL.'DECOUPLAGE'; TABRES.MOTMIX = GTOTA; SINON; TABRES.'SUPERI' = GTOTA*V_SUPE; TABRES.'MEDIAN' = GTOTA*V_MOYE; TABRES.'INFERI' = GTOTA*V_INFE; TABRES.'GLOBAL' = GTOTA; FINSI; FINSI; FINSI; **************************** ** AJOUT AU CHPO OU CHAML ** **************************** SI ((&BCNOEU < NBOU) OU (EGA NBOU 1)); SINON; FINSI; FINSI; FIN BCNOEU; *<== FIN DE BOUCLE SUR LES NOEUDS A AVANCER VIRTUELLEMENT =============| *|=====================================================================| * EN XFEM 3D ON A LE RESULTAT UNIQUEMENT SOUS FORME DE CHAMELEM * ON LE TRANSFORME DONC EN CHPO SINON; * EN STD ON REDUIT LE CHAMP RESULTAT AU FRONT DE FISSURE FINSI; * ON AJOUTE LES COMPOSANTES DANS L'EPAISSEUR DANS LE CAS DES COQUES SI BOOL.'COQ'; GCHPO1 = GCHPO2 ET GCHPO3 ET GCHPO4 ET GCHPO1; FINSI; ********************************************************************* ** STOCKAGE DES CHAMPS RESULTATS (+ logique qu'une table a priori) ** ********************************************************************* *CAS PASAPAS SUPTAB.'CHPO_RESULTATS'.IABC = SUPTAB.'CHPO_RESULTATS'.IABC ET GCHPO1; SUPTAB.'CHAM_RESULTATS'.IABC = SUPTAB.'CHAM_RESULTATS'.IABC ET GCHAM1; FINSI; *CAS RESOU SINON; SUPTAB.'CHPO_RESULTATS' = SUPTAB.'CHPO_RESULTATS' ET GCHPO1; SUPTAB.'CHAM_RESULTATS' = SUPTAB.'CHAM_RESULTATS' ET GCHAM1; FINSI; FINSI; * ON SE PREPARE A PARTIR... MENA; MESS; FINSI; MESS; FINSI; SI (EGA ((&BOUCEXT/10)*10) &BOUCEXT); 'SAUT' 'LIGNE'; FINSI; FINSI; FIN BOUCMIX; * FIN DE BOUCLE SUR LES INTEGRALES A CALCULER ========================* **************************************************** ******* FIN DE BOUCLE SUR LES PAS DE CALCUL ******** **************************************************** FIN BOUCEXT; **************************************************** ** STOCKAGE DES RESULTATS DANS L OBJET EVOLUTIONS ** **************************************************** * en plus de ce qui existe dans SUPTAB.'RESULTATS' * et dans CHPO_RESULTATS * TITRE DES EVOLUTIONS SI BOOL.'DECOUPLAGE'; SINON; CHAINT = MOBJ; SI BOOL.'J_DYNA'; FINSI; FINSI; TITR CHATIT; TINST = SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'; CHPO1 = SUPTAB.'CHPO_RESULTATS'.(&IINST - 1); FIN IINST; * MCOMP EST UN LISTMOT, ON LE TRANSFORME EN MOT S'IL Y EN A QU'UN FINSI; SI ((EGA GDIME 2) OU BOOL.'COQ'); * EN 2D ET 3D COQUE L'INDICE EVOLUTION_RESULTATS EST UNE EVOLUTION SUPTAB.'EVOLUTION_RESULTATS' = EVO1; SINON; * SINON C'EST UNE TABLE INDICEE PAR LES POINTS DU FRONT SUPTAB.'EVOLUTION_RESULTATS'.POIN1 = EVO1; FIN INOEU; * ET IL Y A LE RESULTAT GLOBAL G0 = SUPTAB.'RESULTATS'.(&IINST - 1).'GLOBAL' ; FIN IINST ; SUPTAB.'EVOLUTION_RESULTATS'.'GLOBAL' = EVO1 ; FINSI; FINSI; ******************************************************** **** STOCKAGE POUR UNE EVENTUELLE REPRISE DE CALCUL **** ******************************************************** SUPTAB.'IABC' = IABC; SUPTAB.'MAT_INST' = MAT_INST; SUPTAB.'END1' = ENERM; SUPTAB.'ENV1' = WVMIS; FINSI; SI (EGA ITYPEF 2); SUPTAB.'VDI1' = VDI1; FINSI; FINSI; ******************************************************** ************* OBJETS POUR REPRISE ET XFEM ************** ******************************************************** SUPTAB.'COU1' = SUPTAB.'COUCHE'; FINSI; *bp: soyons coherent... CHAMP_THETA est deja renseigne SUPTAB.'CHAMP_THET1' = SUPTAB.'CHAMP_THETA'; FINSI; *OTER SUPTAB 'CHAMP_THETA'; OTER SUPTAB 'GRTHETA'; OTER SUPTAB 'DIVTHETA'; FINP;
© Cast3M 2003 - Tous droits réservés.
Mentions légales