* fichier : rupt25.dgibi ************************************************************************ * Section : Mecanique Endommagement ************************************************************************ ******************************************************* * Test rupt25.dgibi: Jeux de données * * --------------------------------- * ******************************************************* * CAS TEST DU 15/12/15 PROVENANCE : TEST *Cas test de validation pour le calcul de J sous plusieurs chargement *avec les procedures g_theta.procedur et g_calcul.procedur * *- chargement en traction *- chargement avec pression sur levres *- chargement thermique *Calcul en dimension 3 avec des elements CU20 sur un maillage complet *symetrique opti dime 3 elem cu20 echo 1 ; ************************ *Données paramètriques : ************************ * a : profondeur de la fissure * * t : epaisseur du tube * * ri, re : rayon interne/externe * * h : hauteur du tube * h = 1. ; t = 60.e-3 ; a = t / 5. ; ri = t * 5.; re = ri + t; *POINTS POUR L'AXE DE REVOLUTION p0 = (0. 0. 0.); py = (0. 1. 0.); *COORDONNEE DE LA POINTE DE LA FISSURE pf = (a + ri) 0. 0.; *NOMBRE D'ELEMENTS AUTOUR DE LA POINTE DE LA FISSURE (1 et 2 COUT) nfiss = 10; *TAILLE D'UN ELEMENT DE LA 1ERE ET 2EME COUTURE* tel = 200e-6 ; tel2 = 400e-6 ; *Facteur d'agrandissement de la taille du derafinement ttel2 = 4.*tel2 ; *LONGUEUR DE LA 1ERE ET 2EME COUTURE* lc1 = nfiss * tel ; lc2 = tel2 * nfiss; *NIVEAU DE CHARGEMENT p0T = -400. ; p0P = 400. ; dt0 = 300.; *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * DEBUT DU MAILLAGE *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ************************************************************** ********************** 1ERE COUTURE ************************ ************ (Autour de la pointe de la fissure) *********** ************************************************************** p1cbd = pf plus (lc1 0. 0.) ; p1chd = pf plus (lc1 lc1 0.) ; pf1 = pf plus (0. lc1 0.) ; p1chg = pf1 moin (lc1 0. 0.) ; p1cbg = pf moin (lc1 0. 0.) ; d1ch = droi (2*nfiss) p1chg p1chd; d1cg = droi (nfiss) p1cbg p1chg ; d1cd = droi (nfiss) p1cbd p1chd ; d1cbg = droi (nfiss) p1cbg pf ; d1cbd = droi (nfiss) pf p1cbd ; cout1 = regl nfiss d1ch (d1cbg et d1cbd) ; *cout1 = coul jaun cout1 ; ************************************************************** ********************** 2EME COUTURE ************************ ************ (Autour de la pointe de la fissure) *********** ************** Partie au-dessus de la fissure ************** ************************************************************** p2cbd = pf plus (lc2 0. 0.) ; p2chd = pf plus (lc2 lc2 0.) ; pf2 = pf plus (0. lc2 0.) ; p2chg = pf2 moin (lc2 0. 0.) ; p2cbg = pf moin (lc2 0. 0.) ; d2ch = droi (2*nfiss) p2chg p2chd; d2cg = droi (nfiss) p2cbg p2chg ; d2cd = droi (nfiss) p2cbd p2chd ; cout2 = regl nfiss (d1cd et d1ch et d1cg) (d2cd et d2ch et d2cg) ; *cout2 = coul vert cout2 ; cout1et2 = cout1 et cout2; ************************************************************** ***************** DERAFINEMENT DES COUTURES **************** ************** Partie au-dessus de la fissure ************** ************************************************************** *------------------( DERAF A 4 ELEMENT )------------------- pid1 = p2chg moin (0. ttel2 0.) ; pid2 = pid1 plus (0. tel2 0.) ; pid3 = pid2 plus (0. tel2 0.) ; pid4 = pid3 plus (0. tel2 0.) ; pid5 = pid4 plus (0. tel2 0.) ; pid6 = pid2 moin (tel2 0. 0.) ; pid7 = pid3 moin (tel2 0. 0.) ; pid8 = pid4 moin (tel2 0. 0.) ; pid9 = pid1 moin (ttel2 0. 0.) ; pid10 = pid3 moin (ttel2 0. 0.) ; pid11 = pid5 moin (ttel2 0. 0.) ; did1 = droi 1 pid1 pid2 ; did2 = droi 1 pid2 pid3 ; did3 = droi 1 pid3 pid4 ; did4 = droi 1 pid4 pid5 ; did5 = droi 1 pid9 pid6 ; did6 = droi 1 pid6 pid7 ; did7 = droi 1 pid7 pid8 ; did8 = droi 1 pid11 pid8 ; did9 = droi 1 pid10 pid7 ; si1 = (regl 1 did1 did5) et (regl 1 did2 did6) et (regl 1 did3 did7) et (regl 1 did4 (inve did8)) et (regl 1 did8 did9) et (regl 1 did9 did5) ; elim si1 1.e-5 ; *------------------( DERAF A 3 ELEMENT )------------------- pad1 = pf moin (lc2 0. 0.) ; pad2 = pad1 plus (0. tel2 0.) ; pad3 = pad2 plus (0. tel2 0.) ; pad4 = pad3 plus (0. tel2 0.) ; pad5 = pad2 moin (tel2 0. 0.) ; pad6 = pad3 moin (tel2 0. 0.) ; pad7 = pad1 moin (ttel2 0. 0.) ; pad8 = pad4 moin (ttel2 0. 0.) ; dad1 = droi 1 pad1 pad2 ; dad2 = droi 1 pad2 pad3 ; dad3 = droi 1 pad3 pad4 ; dad4 = droi 1 pad7 pad5 ; dad5 = droi 1 pad5 pad6 ; dad6 = droi 1 pad6 pad8 ; sa1 = (regl 1 dad1 dad4) et (regl 1 dad2 dad5) et (regl 1 dad3 dad6) et (regl 1 dad4 (inve dad6)); saa1 = sa1 ; repe i0 1 ; ssa1 = sa1 plus ( 0. (3.*&i0*tel2) 0.) ; fin i0 ; sa1 = sa1 et ssa1 ; elim sa1 1.e-5 ; *---------------------- PARTIE GAUCHE ----------------------- sig = sa1 et si1 ; elim sig 1.e-5 ; *---------------------- PARTIE DROITE ----------------------- sid = sig syme droi ((coor 1 pf) 0. 0.) ((coor 1 pf) lc2 0.) ; elim sid 1.e-5 ; *---------------------- PARTIE HAUTE ----------------------- *lignes diagonales pour la symetrie p_diagod = p2chd plus (lc1 lc1 0.); p_diago = p2chg moin (lc1 0. 0.); p_diagog = p_diago plus (0. lc1 0.); d_diagog = droi 1 p1chg p_diagog; d_diagod = droi 1 p1chd p_diagod; sihg = sig syme droi p1chg p_diagog ; elim sihg 1.e-5 ; sihd = sid syme droi p1chd p_diagod ; elim sihd 1.e-5 ; sih = sihd et sihg ; elim sih 1.e-5 ; *---------------------- PARTIE COIN ----------------------- dg = droi 1 pid11 p2chg; dcg = dg tran 1 (0. ttel2 0.); dcd = dcg syme droi ((coor 1 pf) 0. 0.) ((coor 1 pf) lc2 0.) ; sic = dcd et dcg; elim sic 1.e-5; cout3 = sig et sid et sih et sic ; elim cout3 1.e-5 ; couttot = cout1 et cout2 et cout3; elim d1cbd couttot 1.e-5 ; ************************************************************** ********************* RESTE DU MAILLAGE ******************** ************** Partie au-dessus de la fissure ************** ************************************************************** *Partie de gauche *---------------- pt1 = mini (coor 1 couttot) 0. 0.; pt2 = (mini (coor 1 couttot)) (maxi (coor 2 couttot)) 0.; ptpartg = cout3 poin droit pt1 pt2 1.e-5 ; d_partg = (cont couttot) elem appuye strictement ptpartg ; pri = ri 0. 0.; pg = d_partg tran (((coor 1 pri)-(mini(coor 1 couttot))) 0. 0.) dini 1.6e-3 dfin 3.2e-3; *Partie de droite *---------------- pt3 = maxi (coor 1 couttot) 0. 0.; pt4 = (maxi (coor 1 couttot)) (maxi (coor 2 couttot)) 0.; ptpartd = couttot poin droit pt3 pt4 1.e-5 ; dpartd = (cont couttot) elem appuye strictement ptpartd ; pre = re 0. 0.; pd = dpartd tran (((coor 1 pre)-(maxi(coor 1 couttot))) 0. 0.) dini 1.6e-3 dfin 3.2e-3; bascout = pg et pd et couttot ; elim bascout 1.e-5 ; *Partie du haut *-------------- p5 = (mini (coor 1 bascout)) (maxi (coor 2 bascout)) 0.; p6 = (maxi (coor 1 bascout)) (maxi (coor 2 bascout)) 0.; pt_parth = bascout poin droit p5 p6 1.e-5 ; dpartd = (cont bascout) elem appuye strictement pt_parth ; ph = dpartd tran (0. ((h/2.) - (maxi(coor 2 bascout))) 0.) dini 1.6e-3 dfin t; *Structure totale EN AXI ************************************* s0 = ph et bascout ; elim s0 1.e-5 ; ************************************* pa = s0 poin proc (ri 0. 0.); pb = s0 poin proc (re 0. 0.) ; lvsup = (cont s0) elem compris pa pf ; *Definition des bords du maillage AXI lhau = cote 3 ph; lbas = (cont s0) elem compris pf pb ; pgau = s0 poin droit (ri (mini (coor 2 s0)) 0.) (ri (maxi (coor 2 s0)) 0.) 1.e-5 ; lgau = (cont s0) elem appuye strictement pgau ; n1 = 1; deg1 = 0.5 ; *Structure totale EN 3D ************************************************* v0 = s0 volu n1 rota deg1 p0 py ; elim v0 1.e-5 ; ************************************************* f1 f2 f3 = face v0; f1 = f1 coul bleu; f2 = f2 coul roug; * Creation des surfaces inferieure et superieure. surliga = lbas rota n1 deg1 p0 py; sursup = lhau rota n1 deg1 p0 py; surlev = lvsup rota n1 deg1 p0 py; *Defintion du front de fissure pfx1 = (coor 1 pf)*(cos (-1.*deg1)); pfx3 = (coor 1 pf)*(sin (-1.*deg1)); pfx = pfx1 0. pfx3; frfiss = cerc n1 pf p0 pfx; frfiss = frfiss coul cyan; elim (v0 et surliga et sursup et surlev et pfx et frfiss et f1 et f2) 1.e-5; *Trois points sur la surface f2 zm = mini (f2 coor 3) ; PC = f2 poin proc (ri 0. zm) ; PD = f2 poin proc (re 0. zm) ; PE = f2 poin proc (ri (h / 2.) zm) ; *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * FIN DU MAILLAGE *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * PARTIE CALCULS *============================================================= ************************************************************** *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * PROPRIETE MATERIAUX A 300°C E0 = 185e3; nu0 = 0.3; alfa0 = 13.08e-6; *Elevation de temperature chx = ((v0 coor 1) - ri) / (re - ri) ; cht0 = nomc 'T' (dt0 * chx); mo0 = mode v0 mecanique elastique isotrope ; Mchtref= CHAN 'CHAM' cht0 mo0 'RIGIDITE' 'CARACTERISTIQUES'; ma0 = mate mo0 YOUN E0 nu nu0 alph alfa0 talp 0. tref 0. ; rg0 = rigi mo0 ma0 ; *CONDITIONS AUX LIMITES *Blocages cl1 = bloq uy surliga ; cl2 = rela ense uy sursup ; cl3 = bloq uz f1 ; cl4 = v0 symt depl PC PD PE 1.e-5 ; cl0 = cl1 et cl2 et cl3 et cl4 ; *Traction uniaxiale (via un modele de pression) moph = MODE sursup 'CHARGEMENT' 'PRESSION' 'CONS' 'HAUT' ; *Pression sur les levres (via un modele de pression) mopl = MODE surlev 'CHARGEMENT' 'PRESSION' 'CONS' 'LEVRES' ; ***************************************** ********* SOLUTIONS ANALYTIQUES ********* ***************************************** *Fonctions d'influence i0 = 1.211 ; i1 = 0.718 ; *Contraintes imposées pour le gradient de temperature sig0 = ((E0*alfa0*dt0)/(1-nu0)) * (ri/(3*t)) * ((2*(re**2))/(ri*(re+ri)) - 1); sig1 = -1. * ((E0*alfa0*dt0)/(1-nu0)); * J analytiques JT = (1-(nu0**2)) * ((i0*(-1.*p0T)*((pi*a)**(1./2.)))**2) / E0; JP = (1-(nu0**2)) * ((i0*(-1.*p0P)*((pi*a)**(1./2.)))**2) / E0; JTH = (1-(nu0**2)) * ((((i0*sig0)+ (i1*sig1*(a/t)))*((pi*a)**(1./2.))) **2)/E0; **************************************************** * CALCUL ELASTIQUE AVEC RESO - CALCUL DE J ELASTIQUE **************************************************** * Construction des second membres maph = pres moph 'PRES' p0T ; f0T = BSIG moph maph ; mapl = pres mopl 'PRES' p0P ; f0P = BSIG mopl mapl ; sgth0 = THET mo0 ma0 cht0 ; f0TH = BSIG mo0 sgth0 ; * RESOLUTION ELASTIQUE DES 3 PROBLEMES utestT utestP utestTH = RESO (rg0 ET cl0) f0T f0P f0TH; mena ; *PROCEDURE G_THETA *cas 1 : traction seule tabJel = table ; tabJel . 'MODELE' = mo0 ET moph ; tabJel . 'CARACTERISTIQUES' = ma0 ; tabJel . 'PRESSION' = maph ; tabJel . 'BLOCAGES_MECANIQUES' = cl0 ; tabJel . 'SOLUTION_RESO' = utestT ; tabJel . 'OBJECTIF' = MOT 'J' ; tabJel . 'LEVRE_SUPERIEURE' = surlev ; tabJel . 'FRONT_FISSURE' = frfiss ; tabJel . 'COUCHE' = 5 ; g_theta tabJel ; JelT1 = tabJel.resultats.global ; *cas 2 : pression sur les levres tabJel . 'MODELE' = mo0 ET mopl ; tabJel . 'CARACTERISTIQUES' = ma0 ; tabJel . 'PRESSION' = mapl ; tabJel . 'SOLUTION_RESO' = utestP ; g_theta tabJel ; JelP1 = tabJel.resultats.global ; *cas 3 : gradient de temperature tabJel . 'MODELE' = mo0 ; tabJel . 'CARACTERISTIQUES' = ma0 ; tabJel . 'TEMPERATURES' = cht0 ; tabJel . 'SOLUTION_RESO' = utestTH ; g_theta tabJel ; JelTH1 = tabJel.resultats.global ; *Erreurs sur J entre la solution analytique et le MEF errT1 = ((JelT1-JT)/JT)*100.; errP1 = ((JelP1-JP)/JP)*100.; errTH1 = ((JelTH1-JTH)/JTH)*100.; ******************************************************* * CALCUL ELASTIQUE AVEC PASAPAS - CALCUL DE J ELASTIQUE ******************************************************* * Chargements de pression (obligatoires si modele de pression) evph = EVOL 'MANU' 'TEMP' (PROG 0. 1. 2. 3.) 'PRES' (PROG 0. 1. 0. 0.) ; chaph = CHAR 'PRES' maph evph ; evpl = EVOL 'MANU' 'TEMP' (PROG 0. 1. 2. 3.) 'PRES' (PROG 0. 0. 1. 0.) ; chapl = CHAR 'PRES' mapl evpl ; * Chargement thermique chath = CHAR 'T' cht0 (EVOL 'MANU' (PROG 0. 1. 2. 3.) (PROG 0. 0. 0. 1.)) ; *RESOLUTION AVEC PASAPAS DES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS) *AU PAS 1 : Traction seule *AU PAS 2 : Pression sur les levres *AU PAS 3 : Gradient de temperature tabT = TABL ; tabT . 'MODELE' = mo0 ET moph ET mopl ; tabT . 'CARACTERISTIQUES' = ma0 ; tabT . 'BLOCAGES_MECANIQUES' = cl0 ; tabT . 'CHARGEMENT' = chaph ET chapl ET chath ; tabT . 'TEMPS_CALCULES' = PROG 1. 2. 3. ; PASAPAS tabT ; *PROCEDURE G_THETA POUR LES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS) *ATTENTION, IL FAUT RETIRER LE CHARGERMENT MECA DE PRESSION SUR LES *LEVRES ET UTILISER LE CHARGEMENT PLEV *PROCEDURE G_THETA POUR LES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS) tabJel = TABL ; tabJel . 'SOLUTION_PASAPAS' = tabT ; tabJel . 'OBJECTIF' = MOT 'J' ; tabJel . 'LEVRE_SUPERIEURE' = surlev ; tabJel . 'FRONT_FISSURE' = frfiss ; tabJel . 'COUCHE' = 5 ; g_theta tabJel ; JelT2 = tabJel.resultats. 1 . global ; JelP2 = tabJel.resultats. 2 . global ; JelTH2 = tabJel.resultats. 3 . global ; *Erreurs sur J : solution analytique VS calcul PASAPAS + G_THETA errT2 = ((JelT2 -JT )/JT )*100.; errP2 = ((JelP2 -JP )/JP )*100.; errTH2 = ((JelTH2-JTH)/JTH)*100.; **************************************** * AFFICHAGE DES RESULTATS ET DES ERREURS **************************************** SAUT 5 'LIGNE' ; mess 'Solution Theorique : ' JT JP JTH ; mess ; mess 'Solution MEF (RESO) : ' JelT1 JelP1 JelTH1 ; mess 'Erreur en % : ' errT1 errP1 errTH1 ; mess ; mess 'Solution MEF (PASAPAS) : ' JelT2 JelP2 JelTH2 ; mess 'Erreur en % : ' errT2 errP2 errTH2 ; * Test sur les erreurs errT = MAXI 'ABS' (PROG errT1 errT2) ; si ((abs errT) > 0.3) ; erre 'Erreur sur le calcul de JelT' ; fins ; errP = MAXI 'ABS' (PROG errP1 errP2) ; si ((abs errP) > 0.3) ; erre 'Erreur sur le calcul de JelP' ; fins ; errTH = MAXI 'ABS' (PROG errTH1 errTH2) ; si ((abs errTH) > 0.4) ; erre 'Erreur sur le calcul de JelTH' ; fins ; FIN ;