* fichier : harttrac3d.dgibi ************************************************************************ * Section : Mecanique Plastique ************************************************************************ *======================================================================* * MODELE HYPERELASTIQUE HART-SMITH QUASI-INCOMPRESSIBLE * * EN GRANDES TRANSFORMATIONS * * * * TEST DE VALIDATION DU MODELE : TRACTION (3D) SIMPLE SELON AXE Z * * COMPARAISON AVEC LA SOLUTION ANALYTIQUE INCOMPRESSIBLE * * * * Contribution de Laurent Gornet - Ecole Centrale de Nantes (2010) * *======================================================================* * Pour plus d'informations, voir la presentation de L. Gornet lors * * du Club Cast3m 2009, disponible sur le site Web de Cast3m. * *======================================================================* * Exemple d'utilisation d'un modele UMAT en grandes transformations * * * * Note : Actuellement en grandes deformations dans PASAPAS, le modele * * ne peut contenir que des modeles de type UMAT. On ne peut * * pas melanger les derivees objectives et les modeles de C3m. * *======================================================================* 'OPTION' 'DIME' 3 'MODE' 'TRID' 'ECHO' 0 ; *opti trac PSC; * Mettre VRAI si l'on souhaite divers traces. GRAPH = VRAI ; GRAPH=FAUX; * Mettre 1 ou 2 selon le type de calcul voulu ICALCUL = 2; * 'SI' ('>' ICALCUL 2) ; 'MESS' 'ICALCUL doit valoir 1 ou 2' ; 'ERREUR' 5 ; 'FINSI' ; title = 'CHAINE' 'HART-SMITH- ' 'TRACTION UNIAXIALE Z' ; *======================================================================* * Geometrie - Maillage * *======================================================================* * Longueur (direction x) du pave : Lg_x = 1. ; * Largeur (direction y) du pave : Lg_y = 1. ; * Hauteur (direction z) du pave : Lg_z = 1. ; * Nombre d'elements selon les directions x, y et z : 'SI' ('EGA' ICALCUL 1) ; Nel_x = 2 ; Nel_y = 2 ; Nel_z = 2 ; * 'TET4' 'PRI6' 'PYR5' 'TET4' 'CUB8' 'OPTION' 'ELEM' 'CUB8' ; 'FINSI' ; 'SI' ('EGA' ICALCUL 2) ; Nel_x = 3 ; Nel_y = 3 ; Nel_z = 3 ; *CU20 'OPTION' 'ELEM' 'CUB8'; 'FINSI' ; *'CU20' P1 = 0. 0. 0. ; P2 = Lg_x 0. 0. ; P3 = Lg_x Lg_y 0. ; P4 = 0. Lg_y 0. ; * L1 = 'DROITE' Nel_x P1 P2 ; L2 = 'DROITE' Nel_y P2 P3 ; L3 = 'DROITE' Nel_x P3 P4 ; L4 = 'DROITE' Nel_y P4 P1 ; * SBas = 'DALLER' L1 L2 L3 L4 ; Pave = 'VOLUME' SBas 'TRANS' Nel_z (0. 0. Lg_z) ; Shau = 'FACE' 2 Pave 'COULEUR' 'BLEU' ; 'SI' GRAPH ; 'TRACER' Pave 'TITRE' ('CHAINE' title ' - MAILLAGE') 'CACH' ; 'FINSI' ; *======================================================================* * Modele - Materiau - Caracteristiques (en MPa) * *======================================================================* * Ne pas oublier de definir les parametres lies a l'elasticite. * Meme si ce n'est pas utilise dans le modele, cela est utile pour * l'operateur de convergence mecanique de PASAPAS-INCREME. * LCMAT = MOTS 'YOUN' 'NU ' 'G' 'K1' 'K2' 'D '; MO = MODE Pave 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'NON_LINEAIRE' 'UTILISATEUR' 'NUME_LOI' 34 'C_MATERIAU' LCMAT ; * * Pour calculer le module d'Young, on utilise les * Coefficients du modele de Mooney-Rivlin (en MPa) : * C1 = 0.183 ; C2 = 0.0034 ; * * On fixe le coefficient de Poisson XNU a une valeur proche de 0.5 * du fait de l'incompressibilite inherente au modele. * Le module de Young YOU est alors connu, car, pour ce modele, le * module de cisaillement MU vaut : MU = YOU/(2*(1+XNU)) = 2.(C1+C2) * Il s'agit de la valeur initiale et de la borne inferieure dans le cas * de la traction. En fonction du niveau de deformation atteinte en * traction, il faut augmenter cette valeur afin de pouvoir faire * converger les calculs (module tangent en fin de calculs). * Prendre des valeurs superieures n'entraine pas de modification des * resultats, cela modifie seulement le nombre d'iterations mecaniques. *XNU = 0.497 ; YOUini = 3.*(2.*(C1+C2)) ; * 'SI' ('EGA' ICALCUL 1) ; XNU = 0.497 ; CoeD = 1.E-4 ; YOU = 10000. * YOUini ; 'FINSI' ; 'SI' ('EGA' ICALCUL 2) ; XNU = 0.497; CoeD = 1.E-4 ; * YOU = 50000. * YOUini ; YOU = 10000. * YOUini ; 'FINSI' ; * *Parametres du modèle HS : essais Treloar/Kawabata MPa G = 0.175 ; K1 = 2.86E-4 ; K2 = 0.311 ; * * * MA = MATE MO 'YOUN' YOU 'NU ' XNU 'G' G 'K1' K1 'K2' K2 'D ' CoeD ; *==================================================================* * Conditions aux limites * *==================================================================* BL1 = 'BLOQUER' 'UZ ' Shau ; BL2 = 'BLOQUER' 'UZ ' Sbas ; BL3 = 'BLOQUER' 'UX ' L4 ; BL4 = 'BLOQUER' 'UY ' L1 ; BLTOT = BL1 'ET' BL2 'ET' BL3 'ET' BL4 ; * * Definition des instants du chargement : t_deb = 0. ; t_fin = 10. ; L_tps = 'PROG' t_deb t_fin ; * Deplacement suivant Z : L_UZ = 'PROG' 0. (6. * Lg_z) ; FF_z = 'DEPIMP' BL1 1. ; EV_z = 'EVOL' 'MANU' 'TEMPS' L_tps 'LAMZ' L_UZ ; CHARTOT = 'CHARGEMENT' 'DIMP' FF_z EV_z ; *======================================================================* * Initialisation de la table pour appel a PASAPAS * *======================================================================* TAB1 = 'TABLE' ; TAB1.'MODELE' = MO ; TAB1.'CARACTERISTIQUES' = MA ; TAB1.'BLOCAGES_MECANIQUES' = BLTOT ; TAB1.'CHARGEMENT' = CHARTOT ; *TAB1.'PAS_AJUSTE'= VRAI; *TAB1.'PRECISION' = 1.E-7 ; *TAB1.'FTOL' = 1.E-5 ; *TAB1.'MTOL' = 1.E-5 ; TAB1.'CONVERGENCE_FORCEE' = FAUX ; TAB1.'GRANDS_DEPLACEMENTS' = VRAI ; TAB1.'HYPOTHESE_DEFORMATIONS' = MOT 'UTILISATEUR' ; TAB1.'TEMPS_CALCULES' = 'PROG' t_deb 'PAS' 0.1 t_fin ; TAB1.'TEMPS_SAUVES' = 'PROG' t_deb 'PAS' 0.5 t_fin ; TAB1.'REAC_GRANDS'=500.; * L_abs = TAB1.'TEMPS_SAUVES' ; n_abs = 'DIMENSION' L_abs ; * PASAPAS TAB1 ; * * Quelques traces de controle apres calculs 'SI' GRAPH ; Defo_0 = 'DEFORMEE' Pave (TAB1.'DEPLACEMENTS'.(n_abs-1)) 0. ; Defo_1 = 'DEFORMEE' Pave (TAB1.'DEPLACEMENTS'.(n_abs-1)) 1. 'VERT' ; 'TRACER' (Defo_0 'ET' Defo_1) 'TITRE' ('CHAINE' title ' - DEFORMEES INITIALE ET FINALE') ; 'TRACER' MO (TAB1.'CONTRAINTES'.(n_abs-1)) 'TITRE' ('CHAINE' title ' - CONTRAINTES EN FIN DE CALCUL') ; 'FINSI' ; * *======================================================================* * Construction de la solution analytique INCOMPRESSIBLE * *======================================================================* * Definitions : * - Allongement selon direction y : Lamz = 1 + (UY/Lg_y) * - Densite d'energie de deformation hyperelastique : W(I1,I2) * - I1, I2 : trois invariants du tenseur de Cauchy-Green droit * Dans le cas du modele de HART-SMITH : Incompressible * dW/dI1 = G * exp (k1 * (I1 - 3)**2) et dW/dI2 = G * k2 * 1 / I2 * Les contraintes de Cauchy sont calculables analytiquement : * - SCzz = 2.(Lamy**2 - 1./Lamy).(dW/dI1 + 1./Lamy.dW/dI2) * L_Un = 'PROG' n_abs '*' 1. ; Lamz = L_Un + (('IPOL' L_abs L_tps L_UZ) / Lg_z) ; * SCxx_th = 0. * L_Un ; SCyy_th = 0. * L_Un ; L_z1 = Lamz * Lamz ; L_z2 = L_Un / Lamz ; L_tr = L_Un * 3.; I1 = L_z1 + (2. * L_z2); I2 = (2. * Lamz) + ( L_Un / L_z1 ); ************************************************************************ dWI1= G * (exp (k1 *((I1 - L_tr)**2))); dWI2= G * k2 * L_Un / I2; ************************************************************************ SCzz_th =(L_z1 - L_z2) * ((2.*dWI1*L_Un) + (2.*dWI2*L_z2)) ; SCxy_th = 0. * L_Un ; SCxz_th = 0. * L_Un ; SCyz_th = 0. * L_Un ; * *======================================================================* * Comparaison des resultats avec la solution analytique * *======================================================================* * La comparaison s'effectue entre les valeurs moyennes des contraintes * calculees et les solutions analytiques correspondantes. * On ne cherche pas a verifier l'uniformite du champ de contraintes. * (Faire le calcul en mettant GRAPH a VRAI et voir les isovaleurs !) * TabD = TAB1.'DEPLACEMENTS' ; TabS = TAB1.'CONTRAINTES' ; Confini = 'FORM' ; ChmUn = 'MANU' 'CHML' MO 'SCAL' 1. ; * SCxx = 'PROG' 0. ; SCyy = 'PROG' 0. ; SCzz = 'PROG' 0. ; SCxy = 'PROG' 0. ; SCxz = 'PROG' 0. ; SCyz = 'PROG' 0. ; 'REPETER' Boucle (n_abs - 1) ; 'FORM' (TabD.&Boucle) ; VolSU = 'INTG' MO ChmUn ; * mess ' volsu ' volsu; SCxx = SCxx 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMXX')/VolSU)) ; SCyy = SCyy 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMYY')/VolSU)) ; SCzz = SCzz 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMZZ')/VolSU)) ; SCxy = SCxy 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMXY')/VolSU)) ; SCxz = SCxz 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMXZ')/VolSU)) ; SCyz = SCyz 'ET' ('PROG' (('INTG' MO (TabS. &Boucle) 'SMYZ')/VolSU)) ; 'FORM' Confini ; 'FIN' Boucle ; * * LG lamb L_abs = Lamz; * 'SI' GRAPH ; tlege = 'TABLE' ; tlege. 1 = 'MARQ CROI' ; tlege.'TITRE' = 'TABLE' ; tlege.'TITRE'. 1 = 'Numerique' ; tlege.'TITRE'. 2 = 'Analytique' ; Evxx = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCXX' SCxx ; Evxx_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCXX' SCxx_th ; 'DESSIN' (Evxx 'ET' Evxx_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY XX (Pa)') ; Evyy = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCYY' SCyy ; Evyy_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCYY' SCyy_th ; 'DESSIN' (Evyy 'ET' Evyy_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY YY (Pa)') ; Evzz = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCZZ' SCzz ; Evzz_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCZZ' SCzz_th ; 'DESSIN' (Evzz 'ET' Evzz_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY ZZ (Pa)') ; Evxy = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCXY' SCxy ; Evxy_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCXY' SCxy_th ; 'DESSIN' (Evxy 'ET' Evxy_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY XY (Pa)'); Evxz = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCXZ' SCxz ; Evxz_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCXZ' SCxz_th ; 'DESSIN' (Evxz 'ET' Evxz_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY XY (Pa)'); Evyz = 'EVOL' 'ROUG' 'MANU' 'LAMB' L_abs 'SCYZ' SCyz ; Evyz_th = 'EVOL' 'BLEU' 'MANU' 'LAMB' L_abs 'SCYZ' SCyz_th ; 'DESSIN' (Evyz 'ET' Evyz_th) 'LEGE' tlege 'TITRE' ('CHAINE' title ' - CONTRAINTE DE CAUCHY XY (Pa)'); 'FINSI' ; * * Tests de bon fonctionnement : Smaxref = 'MAXIMUM' ('ABS' SCzz_th) ; r_xx = 'MAXIMUM' ('ABS' (SCxx - SCxx_th)) / Smaxref ; r_yy = 'MAXIMUM' ('ABS' (SCyy - SCyy_th)) / Smaxref ; r_zz = 'MAXIMUM' ('ABS' (SCzz - SCzz_th)) / Smaxref ; uu= ('ABS' (SCzz - SCzz_th)) / Smaxref; list sczz; list SCzz_th; list uu; r_xy = 'MAXIMUM' ('ABS' (SCxy - SCxy_th)) / Smaxref ; r_xz = 'MAXIMUM' ('ABS' (SCxz - SCxz_th)) / Smaxref ; r_yz = 'MAXIMUM' ('ABS' (SCyz - SCyz_th)) / Smaxref ; * MESS ' RESULTATS : ' title ; MESS ' ------------------------------------------------- '; 'SAUTER' 1 'LIGNE' ; 'MESS' ' Tests de bon fonctionnement :' ; 'MESS' ' -------------------------------' ; 'MESS' ' Comparaison effectuee sur les contraintes de Cauchy' ; 'MESS' ' Ecart relatif maximal entre la valeur moyenne calculee' ; 'MESS' ' et la solution analytique' ; 'MESS' ' Composante XX : ' r_xx ; 'MESS' ' Composante YY : ' r_yy ; 'MESS' ' Composante ZZ : ' r_zz ; 'MESS' ' Composante XY : ' r_xy ; 'MESS' ' Composante XZ : ' r_xz ; 'MESS' ' Composante YZ : ' r_yz ; 'SAUTER' 1 'LIGNE' ; * Ecart relatif maximal tolere 'SI' ('EGA' ICALCUL 1) ; Sigref = 1.E-3 ; 'FINSI' ; 'SI' ('EGA' ICALCUL 2) ; Sigref = 1.E-2 ; 'FINSI' ; * L_z = 'PROG' r_xx r_yy r_zz r_xy r_xz r_yz ; 'SI' ('>EG' ('MAXIMUM' L_z) Sigref) ; 'MESS' ' ---------------------' ; 'MESS' ' ECHEC DU CAS-TEST !' ; 'MESS' ' ---------------------' ; 'ERREUR' 5 ; 'SINON' ; 'MESS' ' ----------------------' ; 'MESS' ' SUCCES DU CAS-TEST !' ; 'MESS' ' ----------------------' ; 'FINSI' ; 'SAUTER' 1 'LIGNE' ; 'FIN' ;