* G_THETA PROCEDUR CB215821 20/07/30 21:15:10 10670
*|=====================================================================|
*| |
*| OBJECTIF : |
*| ========== |
*| |
*| 1) calculer l'integrale caracteristique de mecanique de la rupture |
*| a) J en elasto-plasticite ou en elasto-dynamique pour un |
*| materiau isotrope. |
*| b) dJ/da en elasto-plasticite, utilisable uniquement dans le cas |
*| de materiau isotrope et homogene pour les elements massifs. |
*| c) C* dans le cas de fluage secondaire stationnaire pour un |
*| materiau isotrope. |
*| d) C*H dans le cas de fluage primaire sous un chargement radial |
*| pour un materiau isotrope. |
*| |
*| 2) separer les modes K1, K2 et K3 en elasticite, utilisable |
*| uniquement dans le cas de materiau homogene et isotrope |
*| pour les elements massifs. |
*| |
*| |
*| ENTREE : |
*| ======== |
*| |
*| SUPTAB objet de type TABLE. En entree, SUPTAB sert a definir les |
*| options et les parametres du calcul. Ses indices sont des |
*| objets de type MOTS (a ecrire en toutes lettres) dont voici |
*| la liste : |
*| |
*| |
*| Arguments obligatoires dans tous les cas |
*| ---------------------------------------- |
*| |
*| SUPTAB.'OBJECTIF' = MOT pour preciser le but du calcul, vaut |
*| 1) 'J' pour calculer l'integrale J (ou G), |
*| caracteristique en elasto-plastique. |
*| 2) 'J_DYNA' pour calculer l'integrale J (ou G), |
*| caracteristique en elasto-dynamique. |
*| 3) 'C*' pour calculer l'integrale C*, |
*| caracteristique en fluage secondaire |
*| stationnaire. |
*| 4) 'C*H' pour calculer l'integrale C*(h), |
*| caracteristique en fluage primaire ou |
*| tertiaire. |
*| 5) 'DJ/DA' pour calculer l'integrale de la |
*| derivation dJ/da, caracteristique pour |
*| analyser la stabilite de propagation d'une |
*| fissure ou des fissures interagissantes. |
*| 6) 'DECOUPLAGE' pour decouper les modes mixtes, |
*| c'est a dire la separation des facteurs K1, |
*| K2 (et K3 et 3D). |
*| |
*| SUPTAB.'COUCHE' = ENTIER representant le nombre de couches |
*| d'elements autour du front de la fissure |
*| qui se deplacent pour simuler la propagtion |
*| de la fissure. Il vaut 0 si seul la pointe de |
*| la fissure se deplace, 1 si c'est la premiere |
*| couche d'elements entourant la fissure, 2 si |
*| c'est l'ensemble des premiere et deuxieme couches |
*| d'elements etc. Il convient veiller a ce que |
*| l'ensemble des elements a deplacer n'atteint pas |
*| le bord de la structure fissuree. |
*| Cet argument doit etre absent si l'on souhaite |
*| preciser soi-meme le CHAMP_THETA (cf.8.) |
*| |
*| SUPTAB.'FRONT_FISSURE' = POINT en 2D ou MAILLAGE en 3D massif |
*| representant le front de la fissure. |
*| |
*| |
*| Arguments obligatoires avec des elements standards |
*| -------------------------------------------------- |
*| |
*| SUPTAB.'LEVRE_SUPERIEURE' = Selon la convention de definition, cet |
*| objet (type MAILLAGE) representant la |
*| levre superieure de la fissure. |
*| |
*| SUPTAB.'LEVRE_INFERIEURE' = Selon la convention de definition, cet |
*| objet (type MAILLAGE) representant la |
*| la levre inferieure de la fissure. Si |
*| une seule levre est modelisee, un des |
*| des deux mots ici (LEVRE_SUPERIEURE ou |
*| LEVRE_INFERIEURE) sera suffisant pour |
*| decrire la fissure. |
*| |
*| |
*| Arguments obligatoires avec des elements enrichis (XFEM) |
*| -------------------------------------------------------- |
*| |
*| SUPTAB.'PSI' = 1ere level set (CHPOINT) decrivant la fissure dans |
*| le cas ou l'on utilise des elements XFEM . |
*| SUPTAB.'PHI' = 2eme level set. |
*| |
*| |
*| |
*| Solution obligatoire issus de la procedure PASAPAS |
*| -------------------------------------------------- |
*| |
*| SUPTAB.'SOLUTION_PASAPAS' = TABLE sortant de la procedure PASAPAS. |
*| |
*| |
*| Solution obligatoire issus de l'operateur RESO |
*| ---------------------------------------------- |
*| |
*| SUPTAB.'SOLUTION_RESO' = CHPOINT de deplacement issus de RESO. |
*| SUPTAB.'CARACTERISTIQUES' = Champ de caractristiques matrielles |
*| et eventuellement geometriques |
*| si necessaire. |
*| SUPTAB.'MODELE' = Objet modele (type MMODEL) englobant toute la |
*| structure. |
*| SUPTAB.'TEMPERATURES' = CHPOINT de temperature creant une contrainte|
*| thermique non nulle si elle existe. |
*| SUPTAB.'CHARGEMENTS_MECANIQUES' = CHPOINT representant l'ensemble |
*| des forces exterieures |
*| (surfaciques, volumiques ou |
*| ponctuelles ....) appliquees sur |
*| le systeme si elles existent, SAUF|
*| la pression sur les levres de la |
*| fissure |
*| SUPTAB.'BLOCAGES_MECANIQUES' = RIGIDITE representant le blocages |
*| mecanique du probleme, a fournir |
*| uniquement dans le cas de calcul |
*| de la derivation dJ/da. |
*| |
*| |
*| Arguments optionnels |
*| -------------------- |
*| |
*| |
*| 1 : Materiaux composites (2D massif ou 3D coque seulement) |
*| |
*| SUPTAB.'MODELES_COMPOSITES' = TABLE indicee par des entiers (1 2... |
*| M, M = nombre de Materiaux composites)|
*| pour donner les modeles des materiaux |
*| ayant des discontinutes de proprietes |
*| materielles. |
*| |
*| 2 : Pour un front de fissure tridimensionnel massif |
*| |
*| SUPTAB.'NOEUDS_AVANCES' = MAILLAGE de type POI1 pour donner les |
*| points du front pour lesquels le calcul |
*| sera effectue. Si cet argument est |
*| obsent, le calcul sera fait pour tous |
*| les noeuds sur le front de la fissure. |
*| |
*| 3 : Calcul des termes croises de la matrice dJi/daj |
*| (i non egal a j) dans le cas des fisures interagissantes. |
*| |
*| SUPTAB.'FISSURE_2' = Objet de type MAILLAGE representant une autre |
*| fissure (levres superieure + inferieure si |
*| toutes les deux levres sont presentes). |
*| SUPTAB.'FRONT_FISSURE_2' = POINT ou MAILLAGE reprsentant le front |
*| de la fissure 2 decrite ci-dessus. |
*| |
*| |
*| 4 : Cas d'une fissure circulaire dans une geometrie plane |
*| |
*| SUPTAB.'POINT_CENTRE' = centre de la fissure circulaire |
*| |
*| 5 : Cas ou l'extension de la fissure correspond a une simple |
*| translation dans un tuyauterie droite (3D). Dans ce cas |
*| on effectue dans la procedure CH_THETA une transformation |
*| de tuyau en plaque en passant au systeme de coordonnees |
*| cylindriques. Il est alors necessaire de fournir : |
*| |
*| SUPTAB.'POINT_1' = centre du systeme de coordonnees |
*| SUPTAB.'POINT_2' = POINT tel que l'axe defini par POINT_1 |
*| vers POINT_2 soit l'axe Z poisitif |
*| SUPTAB.'POINT_3' = POINT tel que le plan defini par les 3 points |
*| POINT_1 POINT_2 POINT_3 donne l'angle theta nul |
*| |
*| 6 : Cas ou l'extension de la fissure ne correspond |
*| pas a une simple translation (3D) |
*| |
*| 6.1 Fissure dans un tuyauterie droite (3D, Rotation) |
*| |
*| SUPTAB.'POINT_1' = Objet de type POINT |
*| SUPTAB.'POINT_2' = Objet de type POINT qui, avec le point POINT_1, |
*| constitue l'axe perpendiculaire a la section |
*| fissuree. |
*| |
*| 6.2 Fissure dans un coude (3D, rotation + transformation) |
*| Outre les deux points SUPTAB.'POINT_1' et SUPTAB.'POINT_2' |
*| definis en haut on donne encore : |
*| |
*| SUPTAB.'CHPOINT_TRANSFORMATION' = Objet de type CHPOINT utilise |
*| pour transformer une coude en un |
*| tuyauterie droite. |
*| SUPTAB.'OPERATEUR' = Objet de type MOT valant 'PLUS' ou 'MOIN' pour |
*| indiquer l'operateur PLUS ou MOIN a utiliser |
*| si l'on veut transformer la coude en un |
*| tuyauterie droite. |
*| |
*| 7 : Rotation rigidifiante imposee dans le calcul par PASAPAS |
*| |
*| SUPTAB.'ROTATION_RIGIDIFIANTE' = table indicee par entiers 0,1,2... |
*| donnant les champs de deplacements |
*| due a une rotation rigidifiante de |
*| la piece autour d'un point. Cette |
*| rotation rigidifiante est imposee |
*| dans le calcul par PASAPAS en tant |
*| d'un calcul en grand deplacement. |
*| |
*| 8 : Cas ou on souhaite donner soi-meme le champ THETA |
*| |
*| SUPTAB.'CHAMP_THETA' = Objet de type CHPOINT caracterisant la |
*| propagation de la fissure. Dans ce cas, |
*| ne pas fournir l'indice 'COUCHE' de SUPTAB, |
*| mais fournir 'CHAMP_THETA' a chaque appel. |
*| |
*| 9 : Cas ou on souhaite calculer une integrale dans l epaisseur |
*| d une structure en coque (rapport DMT/96-317) |
*| |
*| On utilise pour cela la technique de multicouche, qui |
*| consiste, avant d'appeler la proceduer G_THETA, a : |
*| 1) Etablir un modele multicouches (cf MODE CONS) sur un ou |
*| des element(s) proche(s) de la fissure sachant qu'il faut |
*| au moins une couche en peau inferieure, une couche en |
*| peau superieure, une couche en ligne moyenne {ces couches |
*| doivent avoir une epaisseur inferieure a 1e-4*(epaisseur |
*| totale de la coque) et donc 2 couches intermediaires. |
*| 2) Penser a donner un excentrement et un nom constituant |
*| different a ces couches. |
*| 3) Assembler le modele multicouches avec le modele du reste |
*| de la structure. |
*| 4) Effectuer le calcul des contraintes et des deplacements |
*| avec le modele total et le materiau qui en decoule. |
*| Le calcul de l'integrale avec la procedure G_THETA sera |
*| realise sur un seul element en multicouche et pour toutes les|
*| couches dans cet element qui ont une epaisseur inferieure a |
*| 1e-4*(epaisseur totale de la coque). Un tel element doit |
*| etre designe par l'argument suivant : |
*| |
*| SUPTAB.'ELEMENT_MULTICOUCHE' = Objet MAILLAGE comportant UN SEUL |
*| element modelise en multicouche. Il |
*| doit etre a l'interieur de la zone |
*| THETA, c'est a dire dans la zone |
*| definie par le nombre SUPTAB.'COUCHE'.
*| Il ne doit pas etre trop loin, ni |
*| trop proche de la pointe de la |
*| fissure. Theoriquement, l'integrale |
*| a calculer est independant du choix |
*| de l'element pres de la fissure, ce |
*| qui est numeriquement verifiable en |
*| la determinant sur des elemens en |
*| multicouche differents. NOTA : Cette |
*| technique necessite un maillage tres |
*| fin dans la zone de la pointe de la |
*| fissure. |
*| |
*| |
*| SORTIE : |
*| ======== |
*| |
*| Les resultats du calcul correspondant a un champ THETA specifie par |
*| l'objet SUPTAB.'COUCHE' (ou SUPTAB.'CHAMP_THETA' dans le cas ou on |
*| souhaite donner soi-meme un champ de type Theta) sont sauves de la |
*| maniere suivante : |
*| |
*| |
*| Dans tous les cas de calcul |
*| --------------------------- |
*| |
*| SUPTAB.'RESULTATS' = Objet contenant la valeur numerique du calcul. |
*| Son type est variable selon qu'on est en 2D ou |
*| 3D et selon la solution du probleme traite : |
*| |
*| 1) valeur de l'integrale de contour dans le cas|
*| d'une solution provenant de l'operateur RESO|
*| 2D => FLOTTANT |
*| 3D massif => TABLE indicee par |
*| .(points au front de fissure) |
*| .'GLOBAL' pour une estimation globale |
*| 3D coque => TABLE indicee par mots |
*| .'SUPERI' en peau superieure |
*| .'INFERI' en peau inferieure |
*| .'MEDIAN' au plan median et |
*| .'GLOBAL' pour une estimation globale |
*| |
*| 2) valeur de l'integrale de contour a un |
*| certain pas du calcul dans le cas d'une |
*| solution provenant de la procedure PASAPAS |
*| 2D => TABLE indicee par |
*| .(numero du pas de calcul) |
*| 3D massif => TABLE indicees par |
*| .(numero du pas de calcul).(points au |
*| front de fissure) |
*| 3D coque => TABLE indicees |
*| .(numero du pas de calcul).'SUPERI' |
*| .(numero du pas de calcul).'INFERI' |
*| .(numero du pas de calcul).'MEDIAN' et |
*| .(numero du pas de calcul).'GLOBAL' |
*| |
*| 3) valeur des F.I.C. (facteurs d'intensite des |
*| contraintes) dans le cas de decouplage des |
*| modes avec une solution provenant de |
*| l'operateur RESO |
*| 2D => TABLE indicee par mots |
*| .'I' pour KI |
*| .'II' pour KII |
*| 3D massif => TABLE indicees par |
*| .'I' .(points au front de fissure) |
*| pour KI |
*| .'II' .(points au front de fissure) |
*| pour KII |
*| .'III'.(points au front de fissure) |
*| pour KIII et |
*| .'GLOBAL'.(points au front de fissure) |
*| |
*| 4) valeur des F.I.C. (facteurs d'intensite des |
*| contraintes) a un certain pas du calcul |
*| dans le cas de decouplage des modes avec |
*| une solution provenant de la procedure |
*| PASAPAS |
*| 2D => TABLE indicees |
*| .'I' .(numero du pas de calcul) pour KI |
*| .'II'.(numero du pas de calcul) pour KII |
*| 3D massif => TABLE indicees par |
*| .'I' .(numero du pas de calcul).(point |
*| au front de fissure) pour KI |
*| .'II' .(numero du pas de calcul).(points |
*| au front de fissure) pour KII |
*| .'III'.(numero du pas de calcul).(points |
*| au front de fissure) pour KIII |
*| |
*| |
*| Dans le cas de calcul effectue pas a pas |
*| ---------------------------------------- |
*| |
*| SUPTAB.'EVOLUTION_RESULTATS' = Objet contenant l'evolution des |
*| resultats en fonction du temps. |
*| Son type est variable selon la |
*| configuration du probleme traite : |
*| |
*| 1) Evolution de l'integrale de contour |
*| 2D => EVOLUTION |
*| 3D massif => TABLE indicee par |
*| .(points au front de fissure) |
*| .'GLOBAL' evolution pour une |
*| estimation globale |
*| 3D coque => TABLE indicee par mots |
*| .'SUPERI' en peau superieure |
*| .'INFERI' en peau inferieure |
*| .'MEDIAN' au plan median et |
*| .'GLOBAL' evolution pour une |
*| estimation globale |
*| |
*| 2) Evolution des F.I.C. (facteurs |
*| d'intensite de contrainte) |
*| 2D => TABLE indicee par |
*| .'I' pour KI |
*| .'II' pour KII |
*| 3D massif => TABLE indicee par |
*| .'I'. (points au front de fissure)
*| .'II'. (points au front de fissure)
*| .'III'.(points au front de fissure)
*| .'GLOBAL' evolution pour une |
*| estimation globale |
*| |
*| |
*| Dans le cas des elements de coque |
*| --------------------------------- |
*| |
*| SUPTAB.'EPAISSEUR_RESULTATS' = representant l'evolution de la valeur|
*| des integrales dans l'epaisseur de la|
*| coque. Son type est variable selon la|
*| solution du probleme traite : |
*| 1) EVOLUTION dans le cas d'une solution |
*| provenant de l'operateur RESO |
*| 2) TABLE indicee par .(numero du pas de |
*| calcul) dans le cas d'une solution |
*| provenant de la procedure PASAPAS |
*| |
*| |
*| Dans le cas de calcul elasto-plastique |
*| -------------------------------------- |
*| |
*| SUPTAB.'CRITERE_DECHARGE' = En cas de calcul elasto-plastique |
*| isotrope ou cinematique, eventuellement thermique, on |
*| calcul un critere de decharge des contraintes defini par |
*| (si, F = courbe de traction ): crit = F(EPSeq)/ SIGeq. |
*| crit = 1. si non-decharge et crit > 1. si decharge. |
*| SUPTAB.'CRITERE_DECHARGE' est une table indicee par les |
*| temps de calcul. |
*| |
*| |
*| Dans le cas du contact frottement sur les levres |
*| -------------------------------------- |
*| |
*| A ce jour, cela est traite pour le cas xfem (these de B.Trolle). |
*| |
*|=====================================================================|
fltrac = faux ;
* fltrac = VRAI ;
flmess = VRAI ;
TRA_PRES = FAUX;
si(flmess);
'DEBUT DE LA PROCEDURE
G_THETA' '
--------------------'
;finsi;
**************************************************
************* INFORMATIONS GENERALES *************
**************************************************
SAUT 1 'LIGNE'
; VALPI
= PI
;
**************************************************
*** QUELQUES MOTS POUR SIMPLIFIER L'ECRITURE ***
**************************************************
MU1
= MOT 'UX'
; MU2
= MOT 'UY'
; MU3
= MOT 'UZ'
;MF1
= MOT 'FX'
; MF2
= MOT 'FY'
; MF3
= MOT 'FZ'
; MU1
= MOT 'UR'
; MU2
= MOT 'UZ'
; MU3
= MOT 'UT'
; MF1
= MOT 'FR'
; MF2
= MOT 'FZ'
;FINSI;
*listmots
* cas axis : faut-il ajouter mu3 ?
SINO;
MU123
= mots MU1 MU2 MU3
; MF123
= mots MF1 MF2 MF3
; MV123
= mots 'VX' 'VY' 'VZ'
;FINS;
**************************************************
***** DONNEES OBLIGATOIRES DANS TOUS LES CAS *****
**************************************************
SI (NON (EXIS SUPTAB 'OBJECTIF'
)); MESS 'ERREUR :IL FAUT SPECIFIER L INTEGRALE'
; SINON;
IINTE = 0;
SI (EGA SUPTAB.'OBJECTIF' 'J');
IINTE = 1;
FINSI;
SI (EGA SUPTAB.'OBJECTIF' 'C*');
IINTE = 2;
FINSI;
SI (EGA SUPTAB.'OBJECTIF' 'C*H');
IINTE = 3;
FINSI;
SI (EGA SUPTAB.'OBJECTIF' 'DJ/DA');
IINTE = 4;
FINSI;
SI (EGA SUPTAB.'OBJECTIF' 'J_DYNA');
IINTE = 5;
FINSI;
SI (EGA SUPTAB.'OBJECTIF' 'DECOUPLAGE');
IINTE = 99;
FINSI;
SI (EGA IINTE 0);
MESS 'ERREUR : ON NE CONNAIT PAS L INTEGRALE SPECIFIEE'
; FINSI;
FINSI;
SI (NON (EXIS SUPTAB 'FRONT_FISSURE'
)); MESS 'ERREUR : ON VEUT LE FRONT DE LA FISSURE'
; FINSI;
MESHFR1 = SUPTAB . 'FRONT_FISSURE';
SI(EGA (TYPE MESHFR1
) 'POINT'
); MESHFR1
= MANU 'POI1' MESHFR1
;FINSI;
MESHFR1
= MESHFR1
COUL 'OLIV'
;
**************************************************
****** TERMES CROISES DE LA MATRICE dJi/daj ******
**************************************************
SI (EGA IINTE 4);
SI ((EXIS SUPTAB 'FISSURE_2'
) OU (EXIS SUPTAB 'FRONT_FISSURE_2'
)); SI (NON (EXIS SUPTAB 'FISSURE_2'
)); MESS ' LES TERMES CROISES DE LA MATRICE'
; FINSI;
SI (NON (EXIS SUPTAB 'FRONT_FISSURE_2'
)); MESS 'ERREUR : ON VEUT AUSSI LE FROND DE LA FISSURE
2 POUR'
; FINSI;
SINON;
SI (EGA SUPTAB.'COUCHE' 0);
MESS 'ERREUR : LE NOMBRE DE COUCHES DOIT ETRE SUPERIEUR A'
; MESS '
0 POUR LE CALCUL DU TERME PRINCIPAL DJi
/DAi'
; FINSI;
FINSI;
FINSI;
**************************************************
****** DONNEES EN CAS DE CALCUL NONLINEAIRE ******
**************************************************
*
*initialisation des valeurs par defaut************
IGDEP = FAUX;
IGDER = FAUX;
IPERSO1 = FAUX;
IMOPRES = FAUX;
IPFISS = FAUX;
SI (EXIS SUPTAB 'SOLUTION_PASAPAS'
); IPAP = VRAI;
* recup du modele mecanique et du materiau associe depuis WTABLE ****
SI (EXIS SUPTAB.'SOLUTION_PASAPAS' 'WTABLE'
); WTAB=SUPTAB.'SOLUTION_PASAPAS'.'WTABLE';
OBJMOD=WTAB.'MOD_MEC';
OBJMAT=WTAB.'MAT_MEC';
SI (EGA IINTE 5);
SI (NON WTAB.'DYNAMIQUE');
MESS 'ERREUR : IL FAUT UNE SOLUTION ELASTO
-DYNAMIQUE'
; FINSI;
FINSI;
* recup du model et du materiau depuis SOLUTION_PASAPAS ****
* rem BP: on ne devrait jamais passer par ici ...
SINON;
MESS 'Absence de WTABLE ! l execution continue ...'
; * on reduit le modele et le materiau au seul comportement mecanique
OBJMOD
= EXTR (SUPTAB.'SOLUTION_PASAPAS'.'MODELE'
) OBJMAT
= REDU (SUPTAB.'SOLUTION_PASAPAS'.'CARACTERISTIQUES'
) OBJMOD
; SI (EGA IINTE 5);
SI (NON SUPTAB.'SOLUTION_PASAPAS'.'DYNAMIQUE');
MESS 'ERREUR : IL FAUT UNE SOLUTION ELASTO
-DYNAMIQUE'
; FINSI;
FINSI;
WTAB= SUPTAB.'SOLUTION_PASAPAS';
FINSI;
SUPTAB . 'MODELE' = OBJMOD;
* recuperation des modeles de pression
SI (EXIS (SUPTAB.'SOLUTION_PASAPAS'.'MODELE'
) '
FORM' 'CHARGEMENT'
); MODCHA
= EXTR (SUPTAB.'SOLUTION_PASAPAS'.'MODELE'
) IMOPRES = VRAI;
FINSI;
FINSI;
SI WTAB.'GRANDS_DEPLACEMENTS';
IGDEP = VRAI;
FINSI;
SI (EXIS SUPTAB 'ROTATION_RIGIDIFIANTE'
); IGDER = VRAI;
FINSI;
* cas particulier de perso1 ou l on ne calcule que ****
* le dernier pas de temps (contenu dans la table estim)
SI (EXIS SUPTAB 'PERSO1'
); IPERSO1 = SUPTAB . 'PERSO1';
si(flmess
); mess 'utilisation de PERSO1 en cours de dvlpt'
; finsi; SI (EXIS SUPTAB.'SOLUTION_PASAPAS' 'ESTIMATION'
); ESTIM = SUPTAB . 'SOLUTION_PASAPAS' . 'ESTIMATION';
SINON;
MESS 'ERREUR : il faut une ESTIMATION
dans la SOLUTION_PASAPAS'
; FINSI;
SI (NON (EXIS SUPTAB 'MAILLAGE_REDUIT'
)); MESS 'Attention! utilisation de PERSO1 sans MAILLAGE_REDUIT'
; MESS 'uniquement valable
dans le cas de fissure stationnaire'
; FINSI;
FINSI;
FINSI;
**************************************************
******** DONNEES EN CAS DE CALCUL LINEAIRE *******
**************************************************
SI (EXIS SUPTAB 'SOLUTION_RESO'
); SI (EGA (TYPE SUPTAB.'SOLUTION_RESO'
) 'CHPOINT '
); IPAP = FAUX;
SI ((EGA IINTE 2) OU (EGA IINTE 3));
MESS 'ERREUR : C
* OU C
*H EST UNE INTEGRALE'
; MESS ' CARACTERISTIQUE EN FLUAGE'
; FINSI;
SI (EGA IINTE 5);
MESS 'ERREUR : IL FAUT UNE SOLUTION DE LA PROCEDURE
PASAPAS'
; FINSI;
SI (NON (EXIS SUPTAB 'CARACTERISTIQUES'
)); MESS 'ERREUR : IL FAUT DONNER LE CHAMP CARACTERISTIQUE'
; FINSI;
SI (NON (EXIS SUPTAB 'MODELE'
)); MESS 'ERREUR : IL FAUT DONNER LE MODELE DE CALCUL'
; FINSI;
OBJMOD
= EXTR (SUPTAB.'MODELE'
) '
FORM' 'MECANIQUE'
; OBJMAT
= REDU (SUPTAB.'CARACTERISTIQUES'
) OBJMOD
; SI (EXIS (SUPTAB.'MODELE'
) '
FORM' 'CHARGEMENT'
); MODCHA
= EXTR (SUPTAB.'MODELE'
) '
FORM' 'CHARGEMENT'
; IMOPRES = VRAI;
SI (NON (EXIS SUPTAB 'PRESSION'
)) ; MESS 'ERREUR : IL MANQUE LA DONNEE DU CHAMP DE PRESSION '
; SINON ;
MATPRE
= REDU (SUPTAB.'PRESSION'
) MODPRE
; FINSI;
FINSI;
FINSI;
SI ((NON (EXIS SUPTAB 'TEMPERATURES'
)) ET (NON (EXIS SUPTAB 'CHARGEMENTS_MECANIQUES'
)) ET (NON IMOPRES));
MESS 'ERREUR : IL FAUT LES CHARGEMENTS APPLIQUES :'
; MESS ' MECANIQUES, THERMIQUES
OU LES DEUX'
; FINSI;
SI ((EGA IINTE 4) ET
(NON (EXIS SUPTAB 'BLOCAGES_MECANIQUES'
))); MESS 'ERREUR : IL FAUT DONNER LE BLOCAGES MECANIQUES'
; FINSI;
FINSI;
FINSI;
MESS 'ERREUR : IL FAUT UNE SOLUTION PROVENANT DE
PASAPAS'
; MESS '
OU DE
RESO POUR DETERMINER L INTEGRALE'
; FINSI;
**************************************************
***** DONNEES EN CAS DE CHARGEMENT THERMIQUE *****
**************************************************
SI IPAP;
CHAR1 = SUPTAB.'SOLUTION_PASAPAS'.'CHARGEMENT';
* ITHER = EXIS CHAR1 'T ';
*bp, 2014-11-13 : ajout distinction cas ITHER et ITHER1
* ITHER=V ITHER1=F <=> pb couples thermo-mecaniques
ITHER1
= (EXIS CHAR1 'T '
) ou (EXIS CHAR1 'TIMP'
); ITHER
= ITHER1
ou (EXIS CHAR1 'Q '
) ou (EXIS CHAR1 'TECO'
) ou (EXIS CHAR1 'TERA'
); SINON;
ITHER
= EXIS SUPTAB 'TEMPERATURES'
;FINSI;
**************************************************************
***** DONNEES EN CAS DE CHARGEMENT DEFORMATIONS IMPOSEES *****
**************************************************************
SI IPAP;
IDEFI
= EXIS CHAR1 'DEFI'
;SINON;
IDEFI
= EXIS SUPTAB 'DEFORMATIONS_IMPOSEES'
;FINSI;
**************************************************
***** DONNEES EN CAS DE CONTACT (ajout BP BT) ****
**************************************************
IFROT=faux;
SI IPAP;
* todo : pas developpe pour l'instant
SINON;
SI (exis SUPTAB 'MODELE_FISSURE'
); IFROT = vrai;
OBJCON = SUPTAB . 'MODELE_FISSURE';
MAICON
= extr OBJCON 'MAILLAGE'
; FINS;
FINSI;
**************************************************
******* TYPE DES ELEMENTS : COQUE OU MASSIF ******
**************************************************
* IPLAN = (EGA &ELEM 'TRI3') OU (EGA &ELEM 'QUA4') OU
* (EGA &ELEM 'TRI6') OU (EGA &ELEM 'QUA8');
*bp: pas tres robuste => on remplace par :
MAILLAGE
= EXTR OBJMOD 'MAIL'
;IPLAN
= (EXIS LELEM1 'TRI3'
) OU (EXIS LELEM1 'QUA4'
) OU (EXIS LELEM1 'TRI6'
) OU (EXIS LELEM1 'QUA8'
); ICOQU
= (&DIME EGA 3) ET IPLAN
;
**************************************************
**** MODELE MULTICOUCHES DANS LE CAS DE COQUE ****
**************************************************
SI ICOQU;
M_DETA
= EXTR OBJMOD 'ZONE'
; SI (EXIS SUPTAB 'ELEMENT_MULTICOUCHE'
); ELMULT = SUPTAB.'ELEMENT_MULTICOUCHE';
SI (NEG (TYPE ELMULT
) 'MAILLAGE'
); MESS 'ERREUR : L ELEMENT EN MULTICOUCHE DOIT'
; FINSI;
SI (NEG (NBEL ELMULT
) 1); MESS 'ERREUR : ON VEUT UN SEUL ELEMENT EN MULTICOUCHE'
; FINSI;
SINON;
MESS 'ERREUR : IL FAUT DESIGNER UN ELEMENT EN MULTICOUCHE'
; FINSI;
M_ELMU
= EXTR (REDU OBJMOD ELMULT
) 'ZONE'
; SI ((DIME M_ELMU
) '
<'
10); MESS 'ERREUR : IL FAUT AU MOINS
3 COUCHES
(peau inf, peau'
; MESS ' sup, ligne moyenne
) D EPAISSEUR INFERIEURE'
; MESS ' A 1E
-4
*(EPAISSEUR TOTALE
) + 2 COUCHES'
; MESS ' INTERMEDIAIRES POUR L ELEMENT DESIGNE EN'
; MESS ' MULTICOUCHE PROCHE DE LA FISSURE.'
; FINSI;
PEX1
= PROG; LMO1
= LECT; MODCOU
= TABLE
; EPAITO
= 0.
; REPETER NBJ5
((DIME M_ELMU
)/2); I1 = (2 * &NBJ5) - 1;
MODCOU.&NBJ5 = M_ELMU.I1;
EP1
= EXTR (REDU MODCOU.
&NBJ5 OBJMAT
) 'EPAI'
1 1 1; EPAITO = EPAITO + EP1;
PEX1
= PEX1
ET (PROG EX1
); LMO1
= LMO1
ET (LECT &NBJ5
); FIN NBJ5;
NSUPE = 0; NMOYE = 0; NINFE = 0;
REPETER NBJ6
(DIME MODCOU
); SI (EGA EX1 (EPAITO/2.) 1.E-4); NSUPE = LM1; FINSI;
SI (EGA EX1 0. 1.E-10); NMOYE = LM1; FINSI;
SI (EGA EX1 (EPAITO/(-2.)) 1.E-4); NINFE = LM1; FINSI;
FIN NBJ6;
SI (EGA NSUPE 0);
MESS 'ERREUR : IL FAUT UNE COUCHE EN PEAU SUPERIEURE'
; MESS ' D EPAISSEUR INFERIEURE A'
; MESS ' 1E
-4
*(EPAISSEUR TOTALE
) '
; FINSI;
SI (EGA NMOYE 0);
MESS 'ERREUR : IL FAUT UNE COUCHE AU PLAN MEDIAN'
; MESS ' AYANT UN EXCENTREMENT NUL'
; FINSI;
SI (EGA NINFE 0);
MESS 'ERREUR : IL FAUT UNE COUCHE EN PEAU INFERIEURE'
; MESS ' D EPAISSEUR INFERIEURE A'
; MESS ' 1E
-4
*(EPAISSEUR TOTALE
) '
; FINSI;
SUPTAB.'EPAISSEUR' = EPAITO;
M_SUPE = MODCOU.NSUPE;
M_MOYE = MODCOU.NMOYE;
M_INFE = MODCOU.NINFE;
FINSI;
**************************************************
** MAILLAGE UTILISE DANS LA RESOLUTION PAR E.F. **
**************************************************
* MAILLAGE = EXTR OBJMOD 'MAIL' ;
* -> bp:fait + haut
SI ICOQU;
TMULT = TABLE;
REPETER NBJ8
((DIME M_DETA
)/2); M1 = M_DETA.(2*&NBJ8);
ITER NBJ8;
FINSI;
TMULT.M1 = VRAI;
REPETER NBJ9
(((DIME M2
)/2) - 1); MAILLAGE
= DIFF MAILLAGE M1
; FIN NBJ9;
FINSI;
FIN NBJ8;
FINSI;
SUPTAB.'MAILLAGE' = MAILLAGE;
**************************************************
*************** TYPES D ELEMENTS ****************
**************************************************
******* ELEMENTS LINEAIRES OU NONLINEAIRES *******
ILIN = EGA NBNO1 2; IQUA = EGA NBNO1 3;
******* ELEMENTS XFEM OU STANDARD ****************
IXFEM
= EXIS OBJMOD '
ELEM' 'XQ4R' 'XC8R'
;
**************************************************
************* DEFINITON DE LA FISSURE ************
**************************************************
******* LEVELSET PSI ET PHI POUR XFEM ************
SI(IXFEM);
SI ((EXIS SUPTAB 'PSI'
) et (EXIS SUPTAB 'PHI'
)); PSI0 = SUPTAB . 'PSI';
PHI0 = SUPTAB . 'PHI';
SINO;
MESS 'ERREUR : ON VEUT PSI
et PHI LEVELSET DE LA FISSURE'
; FINSI;
******* LEVRE_SUPERIEURE ET INFERIEURE POUR STD ***
SINO;
SI (NON (EXIS SUPTAB 'LEVRE_SUPERIEURE'
)); SI (EGA IINTE 99);
MESS 'ERREUR : ON VEUT LA LEVRE SUPERIEURE DE LA FISSURE'
; FINSI;
SI (NON (EXIS SUPTAB 'LEVRE_INFERIEURE'
)); MESS 'ERREUR : IL FAUT DONNER LA FISSURE'
; MESS '
(LEVRE_SUPERIEURE
ou LEVRE_INFERIEURE
ou les
2)'
; SINON;
SUPTAB.'FISSURE' = SUPTAB.'LEVRE_INFERIEURE';
FINSI;
SINON;
SI (NON (EXIS SUPTAB 'LEVRE_INFERIEURE'
)); SI (EGA IINTE 99);
MESS 'ERREUR : ON VEUT LA LEVRE INFERIEURE DE LA FISSURE'
; FINSI;
SUPTAB.'FISSURE' = SUPTAB.'LEVRE_SUPERIEURE';
SINON;
SUPTAB.'FISSURE' = (SUPTAB.'LEVRE_SUPERIEURE') ET
(SUPTAB.'LEVRE_INFERIEURE');
FINSI;
FINSI;
SI IFROT;
MESS 'ERREUR : CONTACT via MODELE_FISSURE avec
XFEM seulement'
; FINS;
FINSI;
****************************************************
******* DETERMINATION DES CHAMPS THETA ET PI *******
******* ET DE LA ZONE DE TRAVAIL ELTETA *******
****************************************************
si(flmess
); mess 'DETERMINATION DES CHAMPS THETA'
; fins;
*** nombre de COUCHE donne => on calcule tout le reste
* CHAMP_THETA + DIRTETA
SI (EXIS SUPTAB 'COUCHE'
); SI(IXFEM);
SUPTAB.'CHAMP_THETA' UTILTETA
= CH_THETX SUPTAB
; SUPTAB.'UTILTET1' = UTILTETA;
SINO;
SUPTAB.'CHAMP_THETA' UTILTETA
= CH_THETA SUPTAB
; SUPTAB.'UTILTET1' = UTILTETA;
SI (EGA IINTE 4);
SI (NON (EXIS SUPTAB 'FRONT_FISSURE_2'
)); SUPTAB.'COUCHE' = (SUPTAB.'COUCHE') - 1;
SUPTAB.'CHAMP_PI' UTILPI
= CH_THETA SUPTAB
; SUPTAB.'COUCHE' = (SUPTAB.'COUCHE') + 1;
SINON;
P1 = SUPTAB.'FRONT_FISSURE';
SUPTAB.'FRONT_FISSURE' = SUPTAB.'FRONT_FISSURE_2';
SUPTAB.'FISSURE' = SUPTAB.'FISSURE_2';
SUPTAB.'CHAMP_PI' UTILPI
= CH_THETA SUPTAB
; SUPTAB.'FRONT_FISSURE' = P1;
FINSI;
FINSI;
FINSI;
* ELTETA = ...
si(exis SUPTAB 'MAILLAGE_REDUIT'
); * ELTETA = MAILLAGE fourni par l utilisateur (attention pas de test de compati
ELTETA = SUPTAB . 'MAILLAGE_REDUIT';
sino;
* ELTETA = MAILLAGE OU TETA N EST PAS NUL + 1 couche
SI (EGA (TYPE ( SUPTAB.'CHAMP_THETA'
)) 'CHPOINT '
); uu
= EXTR SUPTAB.'CHAMP_THETA' 'MAILLAGE'
; SINON;
uu
= EXTR SUPTAB.'CHAMP_THETA'.'GLOBAL' 'MAILLAGE'
; FINSI;
ELTETA
= ELEM MAILLAGE '
APPU' 'LARG' UU
; fins;
*** CHAMP_THETA donne, on calcule DIRTETA sur le front de fissure
*(pour faire simple, on appelle CH_THETA pour cela, mais ce n'est pas economique
SINON;
SI(EXIS SUPTAB 'CHAMP_THETA'
); MESS 'CHAMP_THETA FOURNI PAR L UTILISATEUR'
; q7 = SUPTAB.'CHAMP_THETA';
SI (NEG (TYPE q7
) 'CHPOINT '
); q7= q7 . 'GLOBAL' ;
FINSI;
* ELTETA = ...
si(exis SUPTAB 'MAILLAGE_REDUIT'
); * ELTETA = MAILLAGE fourni par l utilisateur (attention pas de test de compa
ELTETA = SUPTAB . 'MAILLAGE_REDUIT';
sino;
* ELTETA = MAILLAGE OU TETA N EST PAS NUL + 1 couche
ELTETA
= ELEM MAILLAGE '
APPU' 'LARG' UU
; fins;
* UTILTETA = ...
SI(EXIS SUPTAB 'UTILTETA'
); MESS 'UTILTETA FOURNI PAR L UTILISATEUR'
; UTILTETA = SUPTAB . 'UTILTETA';
VECTEUR1 = UTILTETA . 'DIRECTION1';
VECTEUR2 = UTILTETA . 'DIRECTION2';
SI(EGA &DIME 3); VECTEUR3
= UTILTETA . 'DIRECTION3'
; FINSI; SINON;
* DIRECTIONS dans TABUTIL
NV1
= PSCA VECTEUR1 VECTEUR1 MU123 MU123
; VECTEUR1 = (VECTEUR1 / (NV1**0.5))
CHAN 'ATTRIBUT' 'NATURE' 'DIFFUS'
; UTILTETA . 'DIRECTION1' = VECTEUR1;
VECTEUR2
= (-1.
*(EXCO VECTEUR1 'UY' 'UX'
)) et (EXCO VECTEUR1 'UX' 'UY'
); UTILTETA . 'DIRECTION2' = VECTEUR2;
SINO;
* CHT CHN CHB = FRENET SUPTAB.'FRONT_FISSURE';
* VECTEUR1 = -1.*CHN;
* VECTEUR2 = -1.*CHB;
* VECTEUR3 = -1.*CHT;
MESS ' bp: !!!
option non testee, mais on est joueur !!!'
; modfro
= MODE MESHFR1 MECANIQUE ELASTIQUE 'POUT'
; NV3
= PSCA VECTEUR3 VECTEUR3 MU123 MU123
; VECTEUR3 = (VECTEUR3 / (NV3**0.5))
CHAN 'ATTRIBUT' 'NATURE' 'DIFFUS'
; VECTEUR2
= (PVEC VECTEUR3 MU123 VECTEUR1 MU123 MU123
) CHAN 'ATTRIBUT' 'NATURE' 'DIFFUS'
; UTILTETA . 'DIRECTION2' = VECTEUR2;
UTILTETA . 'DIRECTION3' = VECTEUR3;
FINS;
FINSI;
SINO;
*** ni COUCHE ni CHAMP_THETA donne, ERREUR !
MESS 'ERREUR : ON VEUT LE NOMBRE DE COUCHEs D ELEMENTS'
; MESS ' AUTOUR DE LA FISSURE QUI SE DEPLACE'
; MESS '
ou LE CHAMP_THETA'
; MESS ' POUR SIMULER LA PROPAGATION DE LA FISSURE'
; FINSI;
FINSI;
*
si(fltrac);
q7 = SUPTAB.'CHAMP_THETA';
si(neg (type q7
) 'CHPOINT'
); q7
= q7 . 'GLOBAL'
; fins; trac vq7
(MAILLAGE
et MESHFR1
) '
TITR' 'CHAMP_THETA'
; finsi;
**************************************************
************** DIRECTIONS UTILES *****************
**************************************************
*** DIRECTION DE PROPAGATION DE LA FISSURE = DIRTETA
* SI ((EGA &DIME 2) OU ICOQU);
* DIRTETA = UTILTETA . 'DIRECTION';
DIRTETA = UTILTETA . 'DIRECTION1';
* FINSI;
* SI ((EGA &DIME 3) ET (NON ICOQU));
* IND1 = INDE (UTILTETA.'DIRECTION');
* DIRTETA = 0. 0. 0.;
* REPETER BC1 ((DIME IND1) - 1);
* DIRTETA = DIRTETA 'PLUS' (UTILTETA.'DIRECTION'.(IND1.&BC1));
* FIN BC1;
* FINSI;
* DIRTETA = DIRTETA / (NORM DIRTETA);
si(non ICOQU);
DIRNORM = UTILTETA . 'DIRECTION2';
fins;
*** DIRECTION DE CISAILLEMENT SI SEPARATION DE MODES =DIRCISA
*SI ((EGA IINTE 99) ET (EGA &DIME 3) ET (NON ICOQU));
SI ((EGA &DIME 3) ET (NON ICOQU
)); * SI(IXFEM);
DIRCISA = UTILTETA . 'DIRECTION3';
* SINON;
* F1 = PRES 'MASS' OBJMOD SUPTAB.'LEVRE_SUPERIEURE' 1.;
* N1 = NBNO SUPTAB.'FRONT_FISSURE';
* P1 = POIN SUPTAB.'FRONT_FISSURE' ((N1 + 1)/2);
* V1 = EXTR F1 MF1 P1;
* V2 = EXTR F1 MF2 P1;
* V3 = EXTR F1 MF3 P1;
* DIRCISA = PVEC DIRTETA (V1 V2 V3);
* FINSI;
* DIRCISA = DIRCISA / (NORM DIRCISA);
FINSI;
* si(fltrac);
* *si(vrai);
* dx1 = coor (&DIME + 1) MESHFR1;
* * dx1 = maxi (prog 1. ((maxi (resu dx1)) / (nbno MESHFR1)));
* dx1 = (maxi (resu dx1)) / (nbno MESHFR1);
* vdir7 = (VECT dx1 DIRTETA 'DEPL' 'BLEU')
* et (VECT dx1 DIRNORM 'DEPL' 'ROUG');
* si((EGA &DIME 3) ET (NON ICOQU));
* vdir7 = vdir7 et (VECT dx1 DIRCISA 'DEPL' 'VERT');
* trac vdir7 (MESHFR1 et (aret MAILLAGE)) TITR 'DIRECTIONS LOCALES';
* sino;
* trac vdir7 (MESHFR1 et (cont MAILLAGE)) TITR 'DIRECTIONS LOCALES';
* fins;
* fins;
**************************************************
************** ON COMPLETE ELTETA **************
**************************************************
*** ajout eventuel de ELPIa ELTETA ******
SI ((EXIS SUPTAB 'FRONT_FISSURE_2'
) ET (EGA IINTE
4)); ELPI = SUPTAB.'FRONT_FISSURE_2';
REPETER MAIL2 ((SUPTAB.'COUCHE') + 1);
ELPI
= MAILLAGE
ELEM '
APPU' 'LARG' ELPI
; FIN MAIL2 ;
ELTETA = ELTETA ET ELPI;
FINSI;
*** AJOUT DU NOEUD SUPPORT EN DEF.PL.GENERALISEES
SI (EGA &MODE 'PLANGENE'
); ELTETA
= ELTETA
ET (VALE '
MODE' 'PLANGENE'
);FINSI;
ELPOI1
= CHAN ELTETA 'POI1'
;
*** L ELEMENT SUPPORTANT LE MODELE MULTICOUCHE
*** DOIT ETRE DANS LA ZONE THETA
SI ICOQU;
N2
= NBNO (ELTETA
ET (EXTR M_MOYE 'MAIL'
)); SI (NEG N1 N2);
MESS 'ERREUR : L ELEMENT EN MULTICOUCHE DESIGNE POUR
CALCULER'
; MESS ' L INTEGRALE SE TROUVE EN DEHORS DE LA ZONE'
; MESS ' DEFINIE PAR LE NOMBRE DE COUCHES DONNE.'
; FINSI;
FINSI;
**************************************************
*********** TESTER SI REPRISE DE CALCUL **********
**************************************************
*** REPRISE DE CALCUL ? **************************
IREPRI = FAUX;
SI (IPAP ET (NON IPERSO1));
N1
= DIME (SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'
); SI ((EXIS SUPTAB 'IABC'
) ET *bp (EXIS SUPTAB 'COU1') ET
*bp (EXIS SUPTAB 'CHAMP_THETA') ET
(EXIS SUPTAB 'ELTET1'
) ET (EXIS SUPTAB 'RESULTATS'
) ET (EXIS SUPTAB 'EVOLUTION_RESULTATS'
)); IREPRI = '>' (N1 - 1) SUPTAB.'IABC';
mess 'on tente une reprise...'
; FINSI;
FINSI;
*** TESTS DE COMPATIBILITE SI REPRISE DE CALCUL ***
SI IREPRI;
* on verifie que l objectif reste le meme
SI (NEG SUPTAB.'OBJ1' SUPTAB.'OBJECTIF');
MESS 'ERREUR : REPRISE IMPOSSIBLE CAR L OBJECTIF DU'
; MESS ' CALCUL ACTUEL N EST PAS LE MEME QUE'
; MESS ' CELUI DU CALCUL PRECEDENT'
; FINSI;
* on doit avoir le meme nombre de couche (on suppose la fissure fixe)
SI ((EXIS SUPTAB 'COUCHE'
) et (EXIS SUPTAB 'COU1'
)); SI (NEG SUPTAB.'COU1' SUPTAB.'COUCHE');
MESS 'ERREUR : REPRISE IMPOSSIBLE CAR LE NOMBRE DE'
; MESS ' COUCHE ACTUEL N EST PAS LE MEME QUE'
; MESS ' CELUI UTILISE POUR LE CALCUL PRECEDENT'
; FINSI;
FINSI;
* reste a verifier la compatibilite des support de champ teta via elteta
* ELTETA doit etre inclus dans ELTET1
ELTET1 = SUPTAB.'ELTET1';
MESS 'ERREUR : REPRISE IMPOSSIBLE CAR LE SUPPORT DU '
; MESS ' CHAMP_THETA FOURNI N EST PAS INCLUS
DANS'
; MESS ' CELUI UTILISE POUR LE CALCUL PRECEDENT'
; fins;
MESS 'REPRISE DU CALCUL AUTORISE !'
; FINSI;
*** TESTS DE COMPATIBILITE SI UTILISATION DE PERSO1 ***
SI (IPERSO1
et (EXIS SUPTAB 'ELTET1'
)); ELTET1 = SUPTAB.'ELTET1';
MESS 'ERREUR : REPRISE IMPOSSIBLE CAR LE SUPPORT DU '
; MESS ' CHAMP_THETA FOURNI N EST PAS INCLUS
DANS'
; MESS ' CELUI UTILISE POUR LE CALCUL PRECEDENT'
; fins;
si(flmess
); MESS 'POURSUITE DU CALCUL via PERSO1 AUTORISE !'
; fins; FINS;
**************************************************
** MODELES ET MATERIAUX DANS LA ZONE DE TRAVAIL **
**************************************************
*** VERIFICATION DES DONNEES D ENTREE POUR MODELES_COMPOSITES
SI (EXIS SUPTAB 'MODELES_COMPOSITES'
); SI ((&DIME EGA 3) ET (NON ICOQU
)); MESS 'ERREUR : ON NE PEUT ENCORE TRAITER LES PROBLEMES'
; MESS ' DE MATERIAUX COMPOSITES EN 3D MASSIF'
; FINSI;
N1
= DIME SUPTAB.'MODELES_COMPOSITES'
; SI ('<' N1 2);
MESS 'ERREUR : IL FAUT AU MOINS DEUX MODELES POUR'
; MESS ' DETERMINER LA LIGNE COMMUNE D INTERFACE'
; FINSI;
REPETER BIN4 N1;
T1
= TYPE SUPTAB.'MODELES_COMPOSITES'.
&BIN4
; SI (NEG T1 'MMODEL ');
MESS ' TABLE MODELES_COMPOSITES EST INCORRECTE'
; FINSI;
SI (EGA &BIN4 1);
M2
= EXTR SUPTAB.'MODELES_COMPOSITES'.
&BIN4 'MAIL'
; SINON;
M2 = M2 ET
(EXTR SUPTAB.'MODELES_COMPOSITES'.
&BIN4 'MAIL'
); FINSI;
FIN BIN4;
MESS 'ERREUR : TOUS LES MODELES DE MATERIAUX'
; MESS ' COMPOSITES NE SONT PAS DONNES'
; FINSI;
FINSI;
*** CREATION DE OBJMOD ET TABMOD *****************
SI (EXIS SUPTAB 'MODELES_COMPOSITES'
); * CAS DE MODELES COMPOSITES (AVEC DISCONTINUITE) : ON A DU TRAVAIL
REPETER BIN1
(DIME SUPTAB.'MODELES_COMPOSITES'
); M1 = SUPTAB.'MODELES_COMPOSITES'.&BIN1;
N3
= NBNO (ELTETA
ET M2
);* si on a des noeuds en commun, ...
SI (NEG (N1 + N2) N3);
* ... on les recupere
E1
= (DIFF ELPOI1 M2
) DIFF (ELPOI1
ET M2
);* si tous les noeuds en commun sont sur le contour
* => pas d elements a recuperer => on passe au modele suivant
SI (EGA N1 N2);
ITER BIN1;
FINSI;
* sinon, on recupere les elements concernes et le modele reduit
SI (EGA N1 1);
OBJMOD = TABMOD.N1;
SINON;
OBJMOD = OBJMOD ET TABMOD.N1;
FINSI;
FINSI;
FIN BIN1;
SINON;
* CAS DE MODELES SANS DISCONTINUITE : ON A MOINS DE TRAVAIL
OBJMOD
= REDU OBJMOD ELTETA
; TABMOD.1 = OBJMOD;
FINSI;
* list OBJMOD;
OBJMAT
= REDU OBJMAT OBJMOD
;
*** CHAMP EPAISSEUR DANS LA ZONE DE TRAVAIL ******
SI ICOQU;
EPAICH
= (CHAN (EXCO OBJMAT 'EPAI' 'SCAL'
) 'STRESSES' OBJMOD
) CHAN '
TYPE' 'SCALAIRE'
;FINSI;
**************************************************
** CALCUL DE C* PAR DEUX TYPES DE MODELE FLUAGE **
**************************************************
*** ITYPEF = 1 MODELE FLUAGE POUR LEQUEL ON A UNE EXPRESSION
*** EXPLICITE DE L'INTEGRATION DE LA VITESSE DE
*** DEFORMATION DE FLUAGE SUR LE TEMPS
*** ITYPEF = 2 MODELE FLUAGE POUR LEQUEL ON N'OBTIENT PAS
*** FACILEMENT CETTE EXPRESSION EXPLICITE
*** ITYPEF = 99 SI EN ELASTO OU THERMO-ELASTO-PLASTICITE
*** VERIF DES DONNEES ****************************
SI ((EGA IINTE 2) OU (EGA IINTE 3));
MESS 'ERREUR : LA FORMULATION DU PROBLEME NE PERMET'
; FINSI;
FINSI;
SI (EGA IINTE 3);
SI ((EXIS OBJMOD '
MATE' 'FLUAGE' 'BLACKBURN'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'RCCMR_316'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'RCCMR_304'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'POLYNOMIAL'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'LEMAITRE'
)); MESS 'ERREUR : IL FAUT UN MODELE DE FLUAGE NORTON'
; FINSI;
FINSI;
*** DETERMINATION DE ITYPEF *********************
ITYPEF = 99;
SI ((EGA IINTE 2) OU (EGA IINTE 3));
SI ((EXIS OBJMOD '
MATE' 'FLUAGE' 'NORTON'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'POLYNOMIAL'
)); ITYPEF = 1;
FINSI;
SI ((EXIS OBJMOD '
MATE' 'FLUAGE' 'BLACKBURN'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'RCCMR_316'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'RCCMR_304'
) OU (EXIS OBJMOD '
MATE' 'FLUAGE' 'LEMAITRE'
)); ITYPEF = 2;
FINSI;
FINSI;
**************************************************
******* INTERFACES DANS LA ZONE DE TRAVAIL *******
**************************************************
*** CREATION DES INTERFACES INTER-MODELE *********
IPARAL
= VRAI
; LINTER
= TABL;SI ((EXIS SUPTAB 'MODELES_COMPOSITES'
) ET ('
>'
(DIME TABMOD
) 1)); * on boucle sur les modeles qui appartiennent a ELTETA
REPETER BIN2
((DIME TABMOD
) - 1); M1
= EXTR (TABMOD.
&BIN2
) 'MAIL'
; IIN3 = &BIN2;
NIN3
= (DIME TABMOD
) - &BIN2
; REPETER BIN3 NIN3;
IIN3 = IIN3 + 1;
M2
= EXTR (TABMOD .
IIN3) 'MAIL'
;* On itere si (M1 inclut dans M2) ou (M2 inclut dans M1)
ITER BIN3;
FINSI;
* On itere si M1 et M2 n ont pas de noeuds communs
ITER BIN3;
FINSI;
* On recupere l interface M1-M2
* LO1=vrai <=> il existe des noeuds communs a M1 et M2
LO1
= NEG (N1
+ N2
) (NBNO (M1
ET M2
));* LO2=vrai <=> il n'y a pas 1 noeud commun a M1 et M2
LO2
= NEG ('
ABS'
((N1
+ N2
) - (NBNO (M1
ET M2
)))) 1;* LO4 = M1 et M2 forment bien une interface et ne se chevauchent pas
* SI (LO1 ET LO2 ET LO3);
SI (LO1 ET LO2 ET LO4);
* on ajoute l interface car on a >1 noeuds en commun a M1 et M2
LINTER.LE1 = L1;
* IPARAL=vrai <=> toutes les interfaces sont // a la fissure
* rem: si IPARAL=faux, alors il faut ajouter des termes d interfaces au c
*petite modif car DIRTETA doit etre un chpoint desormais (a verifier)...
xDIRTETA
= extr PDIRTETA Presu 'UX'
; yDIRTETA
= extr PDIRTETA Presu 'UY'
; PDIRTETA = xDIRTETA yDIRTETA;
sino;
zDIRTETA
= extr PDIRTETA Presu 'UZ'
; PDIRTETA = xDIRTETA yDIRTETA zDIRTETA;
fins;
PDIRTETA
= PDIRTETA
/ (norm PDIRTETA
); LO1 = ((EGA P1 PDIRTETA 1.E-6) OU
(EGA P1 (-1.*PDIRTETA) 1.E-6));
IPARAL = IPARAL ET LO1;
FINSI;
FIN BIN3;
FIN BIN2;
FINSI;
**** dans le cas decouplage seulement :
**** TEST SI FRONT_FISSURE EST DANS UNE INTERFACE ****
IDANS = FAUX;
SI ((NEG (dime LINTER
) 0) ET (EGA IINTE
99)); REPETER BIN4
(DIME IND1
); LE1 = IND1.&BIN4;
M1
= CHAN LINTER.
LE1 'POI1'
; N2
= NBNO (M1
ET SUPTAB.'FRONT_FISSURE'
); SI (EGA N1 N2);
IDANS = VRAI;
QUIT BIN4;
SINON;
ITER BIN4;
FINSI;
FIN BIN4;
FINSI;
*** SI OUI (IDANS),ON DETERMINE LES MODELES SUP ET INF *******
MODINF = 0; MODSUP = 0;
SI IDANS;
* on redefinit : IPARAL=vrai <=> l'interface a laquelle appartient la fissure e
IPARAL=FAUX;
LSUP = SUPTAB.'LEVRE_SUPERIEURE';
LINF = SUPTAB.'LEVRE_INFERIEURE';
N1SUP
= NBNO (M1
ET LSUP
); N2INF
= NBNO (M2
ET LINF
); SI ( ((N1 + NLSUP - N1SUP) > 1) ET ((N2 + NLINF - N2INF) > 1));
* LSUP et Mod1 ont plus d'1 point commun ET idem pour LINF et Mod2
MODSUP
= TABMOD.
(EXTR LE1
1); MODINF
= TABMOD.
(EXTR LE1
2); SINON;
N1INF
= NBNO (M1
ET LINF
); N2SUP
= NBNO (M2
ET LSUP
); SI ( ((N1 + NLINF - N1INF) > 1) ET ((N2 + NLSUP - N2SUP) > 1));
MODSUP
= TABMOD.
(EXTR LE1
2); MODINF
= TABMOD.
(EXTR LE1
1); SINON;
MESS 'ERREUR : INCOMPATIBILITE ENTRE LE MODELES_COMPOSITES'
; MESS '
ET LES LEVRE_SUPERIEURE
ET _INFERIEURE'
; FINSI;
FINSI;
* LA FISSURE EST BIEN DANS LE PROLONGEMENT DE L' INTERFACE
IPARAL=VRAI;
FINSI;
* REM: il faudrait egalement verifier que MODSUP et MODINF suffisent a decrire E
**************************************************
*** MODPLA : table indicee par entier pour stocker les modeles
*** mecaniques de chaque objet MMODEL.
**************************************************
*** Elle est vide si le modele est elastique ou elastoplastique
*** avec une courbe de traction independante de la temperature.
*** Dans le cas contraire la table vaut :
*** 1 si le modele est plastique isotrope. Alors une
*** nouvelle courbe de traction EPSE-SIGMA est faite.
*** 2 si le modele est plastique cinematique
*** 3 si le modele est plastique parfaite
YOUVARI = FAUX; NUVARI = FAUX;
ALFVARI = FAUX; MODPLA = TABLE; TABTRA = TABLE;
REPETER BCMOD1 NBOBJ;
*
*** YOUVARI **************************************
SI (EGA TYPYO 'EVOLUTIO');
YOUVARI = VRAI;
SINON;
SI (TEST1 '>' 1.E-10);
YOUVARI = VRAI;
FINSI;
FINSI;
*
*** NUVARI **************************************
SI (EGA TYPNU 'EVOLUTIO');
NUVARI = VRAI;
SINON;
SI (TEST1 '>' 1.E-10);
NUVARI = VRAI;
FINSI;
FINSI;
*
*** ALFVARI **************************************
SI ITHER;
SI (EGA TYPAL 'EVOLUTIO');
ALFVARI = VRAI;
SINON;
SI (TEST1 '>' 1.E-10);
ALFVARI = VRAI;
FINSI;
FINSI;
FINSI;
*
*** courbe de TRACtion ***************************
SI (EGA TYPTR 'NUAGE ');
MODPLA.&BCMOD1 = 1;
REPETER BNUA1
(DIME TRA1 'UPLE'
); SI (EGA &BNUA1 1);
SINON;
NUA1
= EXTR TRA1 '
SUPE' NOMFLO1
(T1
+ 1.
E-10
); FINSI;
VYOU1
= (EXTR 2 PSIG1
) / (EXTR 2 PEPS1
); REPETER BSIG1
((DIME PSIG1
) - 1); VA1
= (EXTR (&BSIG1
+ 1) PEPS1
) - ((EXTR (&BSIG1
+ 1) PSIG1
) / VYOU1
); PEPS2
= PEPS2
ET (PROG VA1
); FIN BSIG1;
SI (&BNUA1 EGA 1);
SINON;
FINSI;
FIN BNUA1;
TABTRA.&BCMOD1 = TRA2;
*** On enleve la courbe de traction si elle depend de
*** la temperature (operation trop couteuse pour VARI)
REPETER BCOM1
(DIME LCOMP1
); SI (EGA &BCOM1 1);
SINON;
FINSI;
FINSI;
FIN BCOM1;
FINSI;
*** SIGY (et pas TRAC) *************************
SINON;
SI ((EGA TYPH 'EVOLUTIO') OU
(EGA TYPSI 'EVOLUTIO'));
MODPLA.&BCMOD1 = 2;
FINSI;
SINON;
SI (EGA TYPSI 'EVOLUTIO');
MODPLA.&BCMOD1 = 3;
FINSI;
FINSI;
FINSI;
FINSI;
FIN BCMOD1;
MATVARI
= YOUVARI
OU NUVARI
OU ALFVARI
OU ((DIME MODPLA
) '
>'
0);
**************************************************
************ CAS IMPOSSIBLE A TRAITER ************
**************************************************
SI (EXIS SUPTAB 'MODELES_COMPOSITES'
); * SI ((EGA IINTE 99) ET (NON IPARAL));
SI ((EGA IINTE 99) ET (IDANS ET (NON IPARAL)));
MESS 'ERREUR : ON NE PEUT ENCORE DECOUPLER LES MODES'
; MESS '
DANS LE CAS DES MATERIAUX COMPOSITES'
; MESS '
SI LA FISSURE N APPARTIENT PAS A L INTERFACE'
; FINSI;
* SI ((EGA IINTE 4) ET (NON IPARAL));
SI ((EGA IINTE
4) ET (NEG (dime LINTER
) 0)); MESS '
DANS LE CAS DES MATERIAUX COMPOSITES'
; FINSI;
SINON;
SI ((EGA IINTE 99) ET MATVARI);
MESS 'ERREUR : ON NE PEUT DECOUPLER LES MODES'
; MESS '
DANS LE CAS DES CARACTERISTIQUES'
; MESS ' MATERIELS VARIABLES
DANS L ESPACE'
; FINSI;
SI ((EGA IINTE 4) ET MATVARI);
MESS ' LE CAS DES CARACTERISTIQUES MATERIELS'
; FINSI;
FINSI;
*
SI ((EGA IINTE 99) ET ICOQU);
MESS 'ERREUR : ON NE PEUT ENCORE DECOUPER LES MODES'
; MESS '
DANS LE CAS DES ELEMENTS DE COQUE'
; FINSI;
SI ((EGA IINTE 4) ET ICOQU);
MESS '
DANS LE CAS DES ELEMENTS DE COQUE'
; FINSI;
SI ((EGA IINTE 3) ET ITHER);
MESS ' C
*(H
) DANS LE CAS DE CHARGEMENT THERMIQUE'
; FINSI;
SI ((EGA IINTE 2) ET ITHER);
MESS ' C
* DANS LE CAS DE CHARGEMENT THERMIQUE'
; FINSI;
*
* SI (IGDEP ET (NEG IINTE 1));
*BP : on propose de calculer malgre tout...
SI (IGDEP ET ((NEG IINTE 1) et (NEG IINTE 5) et (NEG IINTE 99)));
MESS ' SPECIFIEE EN GRANDS
-DEPLACEMENTS'
; FINSI;
**************************************************
*** TITRES A AFFICHER SELON LE PROBLEME TRAITER **
**************************************************
SI IPAP;
TXMECANI
= MOT ' Mecanique'
; TXTERMI
= MOT ' Thermique'
; TXPRESS
= MOT ' Volumique'
;SINON;
TXMECANI
= MOT ' Mecanique'
; TXTERMI
= MOT ' Thermique'
; TXPRESS
= MOT ' Volumique'
;FINSI;
MC1
= MOT ' EN VISCO
-THERMO
-PLASTIQUE'
;MC2
= MOT ' EN VISCO
-PLASTICITE'
;MC3
= MOT ' EN THERMO
-PLASTICITE'
;MC4
= MOT ' EN ELASTO
-PLASTICITE'
;si(exis SUPTAB 'COUCHE'
); MC10
= CHAI '
(Theta ' SUPTAB.'COUCHE' '
)'
;sino;
MC10
= CHAI '
(Theta utilisateur
)'
;fins;
SI (EGA IINTE 1);
CHA1
= CHAI 'INTEGRALE J EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' INTEGRALE J' MC3
; SINON;
TX1
= CHAI ' INTEGRALE J' MC4
; FINSI;
FINSI;
SI (EGA IINTE 2);
CHA1
= CHAI 'INTEGRALE C
* EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' INTEGRALE C
*' MC1
; SINON;
TX1
= CHAI ' INTEGRALE C
*' MC2
; FINSI;
FINSI;
SI (EGA IINTE 3);
CHA1
= CHAI 'INTEGRALE C
*H EN FONCTION DU TEMPS'
;
SI ITHER;
TX1
= CHAI ' INTEGRALE C
*H' MC1
; SINON;
TX1
= CHAI ' INTEGRALE C
*H' MC2
; FINSI;
FINSI;
SI (EGA IINTE 4);
SI (EXIS SUPTAB 'FISSURE_2'
); CHA1
= CHAI 'INTEGRALE CROISEE DJi
/DAj EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' INTEGRALE CROISEE DJi
/DAj' MC3
; SINON;
TX1
= CHAI ' INTEGRALE CROISEE DJi
/DAj' MC4
; FINSI;
SINON;
CHA1
= CHAI 'INTEGRALE DJ
/DA EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' INTEGRALE DJ
/DA' MC3
; SINON;
TX1
= CHAI ' INTEGRALE DJ
/DA' MC4
; FINSI;
FINSI;
FINSI;
SI (EGA IINTE 5);
CHA1
= CHAI 'INTEGRALE J DYNAMIQUE EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' INTEGRALE J EN THERMO
-ELASTO
-DYNAMIQUE'
; SINON;
TX1
= CHAI ' INTEGRALE J EN ELASTO
-DYNAMIQUE'
; FINSI;
FINSI;
SI (EGA IINTE 99);
CHA1
=CHAI 'F.
I.
C.
Ki EN FONCTION DU TEMPS'
; SI ITHER;
TX1
= CHAI ' SEPARATION DES F.
I.
C.' MC3
; SINON;
TX1
= CHAI ' SEPARATION DES F.
I.
C.' MC4
; FINSI;
FINSI;
*
TX2
= CHAI ' Contribution due au chargement' MC10
;TX3
= CHAI ' ________________________________________'
;si(exis SUPTAB 'COUCHE'
); SI ('>' SUPTAB.'COUCHE' 9);
TX3
= CHAI ' _________________________________________'
; FINSI;
SI ('>' SUPTAB.'COUCHE' 99);
TX3
= CHAI ' __________________________________________'
; FINSI;
sino;
TX3
= CHAI ' _____________________________________________'
;fins;
***************************************************
**************** AFFICHAGE DU TITRE **************
***************************************************
SI IPAP;
SI (EGA IINTE 99);
MESS '
Mode No.
Pas ' TXMECANI TXTERMI TXPRESS MOTCO
; SINON;
MESS ' No.
Pas ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
FINSI;
SI ((EGA &DIME 3) ET (NON ICOQU
)); SI (EGA IINTE 99);
MESS '
Mode Noeud No.
Pas ' TXMECANI TXTERMI TXPRESS MOTCO
; SINON;
MESS ' Noeud No.
Pas ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
FINSI;
SI ICOQU;
MESS ' Plan No.
Pas ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
SINON;
SI (EGA IINTE 99);
MESS '
Mode ' TXMECANI TXTERMI TXPRESS MOTCO
; SINON;
MESS ' ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
FINSI;
SI ((EGA &DIME 3) ET (NON ICOQU
)); SI (EGA IINTE 99);
MESS '
Mode Noeud ' TXMECANI TXTERMI TXPRESS MOTCO
; SINON;
MESS ' Noeud ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
FINSI;
SI ICOQU;
MESS ' Plan ' TXMECANI TXTERMI TXPRESS MOTCO
; FINSI;
FINSI;
*******************************************************
* CONDITIONS AUX LIMITES POUR LE DECOUPLAGE DES MODES *
* utile dans le cas de l utilisation d une METHODE MECANIQUE
* pour la creation des champs auxiliaires
*******************************************************
SI (EGA IINTE 99);
PM = SUPTAB.'FRONT_FISSURE';
SI ((EGA &DIME 3) ET (NON IDANS
)); * SI (EGA &DIME 2);
* X1 Y1 = COOR MAILLAGE; X0 Y0 = COOR PM;
* DIS1 = (((X1 - X0)**2) + ((Y1 - Y0)**2))**0.5;
* SINON;
* X1 Y1 Z1 = COOR MAILLAGE;
SI (NON ICOQU
); PM
= POIN PM 'INIT'
; FINSI; * X0 Y0 Z0 = COOR PM;
* DIS1 = (((X1 - X0)**2) + ((Y1 - Y0)**2) + ((Z1 - Z0)**2))**0.5;
* FINSI;
* PLOIN1 = POIN 1 (POIN 'MAXI' DIS1);
* SI (EGA &DIME 2);
* X1 Y1 = COOR (MAILLAGE DIFF (SUPTAB.'LEVRE_INFERIEURE'
* ET SUPTAB.'LEVRE_SUPERIEURE'));
* X0 Y0 = COOR PLOIN1;
* DIS1 = (((X1 - X0)**2) + ((Y1 - Y0)**2))**0.5;
* SINON;
* X1 Y1 Z1 = COOR (MAILLAGE DIFF (SUPTAB.'LEVRE_INFERIEURE'
* ET SUPTAB.'LEVRE_SUPERIEURE'));
* X0 Y0 Z0 = COOR PLOIN1;
* DIS1 = (((X1 - X0)**2) + ((Y1 - Y0)**2) + ((Z1 - Z0)**2))**0.5;
* FINSI;
* PLOIN2 = POIN 1 (POIN 'MAXI' DIS1);
* BLOQ0 = (BLOQ 'DEPL' 'ROTA' PLOIN2) ET (BLOQ MU2 PLOIN1);
*BP : on fait + simple pour les CL utiles au calcul des champs aux.
* + tard on fera probablement du tout analytique comme pour les xfem
* car la methode mecanique n est valable que pour fissure plane et
* front rectiligne a cause de la direction du chargement difficile a
* definir sinon
si(NON ICOQU);
sino;
fins;
FINSI;
FINSI;
******************************************
* FISSURE DANS LE REPERE GLOBAL ET LOCAL *
******************************************
* ELEMENT FINI STANDARD ******************
SI ((EGA IINTE 99) et (non IXFEM));
*** CAS 2D *******************************
* Inclinaison de la fissure par rapport a l'axe global
SEG1
= ORDO (SUPTAB.'LEVRE_SUPERIEURE'
ELEM '
APPU' 'LARG' PM
) ; P_SUP
= POIN SEG1 'INIT'
; SI (EGA P_SUP PM);
P_SUP
= POIN SEG1 'FINA'
; FINSI;
SEG1
= ORDO (SUPTAB.'LEVRE_INFERIEURE'
ELEM '
APPU' 'LARG' PM
) ; P_INF
= POIN SEG1 'INIT'
; SI (EGA P_INF PM);
P_INF
= POIN SEG1 'FINA'
; FINSI;
ALPHA1 = ATG (YG0 - ((YP1 + YP2)/2.)) (XG0 - ((XP1 + XP2)/2.));
* Coordonnees dans le repere Global et Local
XL1 = ((XG1 - XG0)*(COS ALPHA1)) + ((YG1 - YG0)*(SIN ALPHA1));
YL1 = ((YG1 - YG0)*(COS ALPHA1)) - ((XG1 - XG0)*(SIN ALPHA1));
SI ('
<'
('
MESU'
('
DROI'
1 P_SUP P_INF
)) 1.
E-10
); L1
= SUPTAB.'LEVRE_SUPERIEURE'
ELEM '
APPU' ELTETA
; L2
= SUPTAB.'LEVRE_INFERIEURE'
ELEM '
APPU' ELTETA
; YL1 = YL1 + C1 + C2;
FINSI;
* Coordonnees cylindriques RAY1 TETA1 (1.E-30 pour eviter erreur atg 0 0
TETA1 = ATG YL1 (XL1 + 1.E-30);
RAY1 = (((XL1*XL1) + (YL1*YL1))**0.5) + 1.E-10;
* On inverse afin d'avoir YL1 > 0 pour modsup et <0 pour modinf
SI (('<' VA1 0.) ET ('>' VA2 0.));
PPPP = P_SUP; P_SUP = P_INF; P_INF = PPPP;
MMDD = MODSUP; MODSUP = MODINF; MODINF = MMDD;
FINSI;
SI ((EGA XP1 XP2 1.E-10) ET (EGA YP1 YP2 1.E-10));
TETA_S
= REDU TETA1 SUPTAB.'LEVRE_SUPERIEURE'
; TETA_F
= REDU TETA1 SUPTAB.'LEVRE_INFERIEURE'
; TETA1 = TETA1 - TETA_S - TETA_F;
* On se debrouille pour avoir exactement +/-180 sur levre sup/inf
SI (('>' VA1 0.) ET ('<' VA2 0.));
TETA1 = TETA1 + ((TETA_S*0.) + 180.) + ((TETA_F*0.) - 180.);
SINON;
TETA1 = TETA1 + ((TETA_F*0.) + 180.) + ((TETA_S*0.) - 180.);
FINSI;
FINSI;
* valeur en radian
TETA1rad = TETA1*VALPI/180.;
* cas d une interface: on construit 3 zones :
* PM1 = points du modele 1 / interface
* PM2 = points du modele 2 / interface
* L1 = elem de l interface
SI IDANS;
FINSI;
* si(fltrac);
* trac RAY1 ELTETA 'TITR' 'RAY1';
* trac TETA1 ELTETA 'TITR' 'TETA1' (prog -180. PAS 20. 180.);
* fins;
*** CAS 3D *******************************
SINO;
mess 'on ne fait pas ici le passage local global en 3D ...?'
; mess '
et on ne calule pas RAY1 TETA1
non plus?'
; FINSI;
FINSI;
* XFEM************************************
SI ((EGA IINTE 99) et (IXFEM));
*** CAS 2D *******************************
* On recupere les level set
LV7
= (PSI1
NOMC 'UX'
) ET (PHI1
NOMC 'UY'
);*trac GLV7 objmod 'TITR' 'GLV7';
* ce repere est il direct? (si oui/non, SDIR1=+/-1)
GLV7po
= CHAN '
CHPO' OBJMOD GLV7 'MOYE'
; XDIR1
= ((EXCO GLV7po 'UX,X' 'SCAL'
) * (EXCO GLV7po 'UY,Y'
)) - ((EXCO GLV7po 'UX,Y' 'SCAL'
) * (EXCO GLV7po 'UY,X'
)); * Angle ALPHA1 de passage local -> global
NGPSI1
= ( ((EXCO GLV7 'UX,X' 'SCAL'
)**2) + ((EXCO GLV7 'UX,Y' 'SCAL'
)**2) )**(-0.5) ; NGPHI1
= ( ((EXCO GLV7 'UY,X' 'SCAL'
)**2) + ((EXCO GLV7 'UY,Y' 'SCAL'
)**2) )**(-0.5) ; COS1A
= 0.5 * ( ( (EXCO GLV7 'UX,X' 'SCAL'
) * NGPSI1
) + (SDIR1
* ( (EXCO GLV7 'UY,Y' 'SCAL'
) * NGPHI1
)) ); SIN1A
= 0.5 * ( ( (EXCO GLV7 'UX,Y' 'SCAL'
) * NGPSI1
) - (SDIR1
* ( (EXCO GLV7 'UY,X' 'SCAL'
) * NGPHI1
)) ); ALPHA1
= (MASQ SIN1A 'EGSUPE'
0.
) * (ACOS COS1A
);*trac ALPHA1 objmod 'TITR' 'ALPHA1';
COSA2 = COS1A ** 2 ;
SINA2 = SIN1A ** 2 ;
SINCOSA = SIN1A * COS1A ;
* repere local de la fissure
XL1
= (NOMC PSI1 'SCAL'
) CHAN 'CHAM' OBJMOD 'STRESSES'
; YL1
= (SDIR1
* (NOMC PHI1 'SCAL'
)) CHAN 'CHAM' OBJMOD 'STRESSES'
;* Coordonnees cylindriques RAY1 TETA1
TETA1
= CHAN (ATG YL1
(XL1
+ 1.
D-30)) '
TYPE' 'SCALAIRE'
; RAY1 = ((XL1**2) + (YL1**2))**0.5;
*** CAS 3D *******************************
SINO;
mess 'DECOUPLAGE 3D
XFEM en cours de dvpt...'
; * On recupere les level set
LV7
= (PSI1
NOMC 'UX'
) ET (PHI1
NOMC 'UY'
) et (MANU '
CHPO' ELTETA
1 'UZ '
0. 'NATURE' 'DIFFUS'
); SDIR1 = 1.;
* creation de la matrice de rotation
V1 = UTILTETA . 'V1';
V2 = UTILTETA . 'V2';
V3 = UTILTETA . 'V3';
ROT1
= (EXCO V1
(mots 'UX' 'UY' 'UZ'
) (mots 'UX,X' 'UY,X' 'UZ,X'
)) ET (EXCO V2
(mots 'UX' 'UY' 'UZ'
) (mots 'UX,Y' 'UY,Y' 'UZ,Y'
)) ET (EXCO V3
(mots 'UX' 'UY' 'UZ'
) (mots 'UX,Z' 'UY,Z' 'UZ,Z'
)); (mots UX,X UX,Y UX,Z UY,X UY,Y UY,Z UZ,X UZ,Y UZ,Z
); ROT1
= CHAN 'CHAM' ROT1 OBJMOD 'STRESSES' 'GRADIENT'
;* SUPTAB . 'ROT1' = ROT1;
* repere local de la fissure
XL1
= (NOMC PSI1 'SCAL'
) CHAN 'CHAM' OBJMOD 'STRESSES'
; YL1
= (NOMC PHI1 'SCAL'
) CHAN 'CHAM' OBJMOD 'STRESSES'
;* Coordonnees cylindriques RAY1 TETA1
TETA1
= CHAN (ATG YL1
(XL1
+ 1.
D-30)) '
TYPE' 'SCALAIRE'
; RAY1 = ((XL1**2) + (YL1**2))**0.5;
FINSI;
* si(fltrac);
* trac ROT1 OBJMOD 'TITR' 'ROT1';
* trac RAY1 objmod 'TITR' 'RAY1';
* trac TETA1 objmod 'TITR' 'TETA1' (prog -180. PAS 20. 180.);
* finsi;
FINSI;
***************************************************
****** FORCE, DEPLACEMENT ET GRADIENT NULS ******
***************************************************
ZER1
= MANU 'CHML' OBJMOD 'SCAL'
0. '
TYPE' 'SCALAIRE' 'STRESSES'
;SI (EGA &MODE 'PLANGENE'
); 'FY'
0. '
TITR' 'FORCES ' 'NATURE' 'DIFFUS'
; 3 'FZ' 0. 'MX' 0. 'MY' 0.
'
TITR' 'FORCES ' 'NATURE' 'DIFFUS'
); 'UY'
0. '
TITR' 'DEPLACEM' 'NATURE' 'DIFFUS'
; 3 'UZ' 0. 'RX' 0. 'RY' 0.
'
TITR' 'DEPLACEM' 'NATURE' 'DIFFUS'
);FINSI;
CMD000
= CHAN 'NOEUD' OBJMOD
('
ZERO' OBJMOD 'DEPLACEM'
);CMD001
= CHAN 'STRESSES' OBJMOD
('
ZERO' OBJMOD 'DEPLACEM'
);GRA000
= '
ZERO' OBJMOD 'GRADIENT'
;VAR000
= '
ZERO' OBJMOD 'VARINTER'
;
**************************************************
* NOMBRE DE BOUCLE POUR LE CALCUL DES INTEGRALES *
**************************************************
SI IPAP;
NBG = -1;
NBDEP
= DIME (SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'
); SI IREPRI ;
NBG = SUPTAB.'IABC';
NBDEP = NBDEP - 1 - NBG;
FINSI;
SI IPERSO1;
NBG = (WTAB . 'PAS') - 1;
NBDEP = 1;
FINSI;
SINON;
NBG = -1;
NBDEP = 1;
FINSI;
***************************************************
** SOLUTION DU PAS PRECEDENT SI REPRISE DE CALCUL * (ou si perso1)
***************************************************
SI (IREPRI ou (IPERSO1 et (NBG >eg 0))) ;
SIG1 = (SUPTAB.'SOLUTION_PASAPAS'.'CONTRAINTES'.NBG)
SI (EXIS (SUPTAB.'SOLUTION_PASAPAS'
) 'VARIABLES_INTERNES'
); VAR1 = (SUPTAB.'SOLUTION_PASAPAS'.'VARIABLES_INTERNES'.NBG)
SINON;
VAR1 = VAR000;
FINSI;
MAT1 = SUPTAB.'MAT1';
WELAS
= 0.5*('
ENER' OBJMOD SIG1
('
ELAS' OBJMOD SIG1 MAT1
));* WPLAS = SUPTAB.'END1' - WELAS;
WPLAS
= (redu OBJMOD SUPTAB.'END1'
) - WELAS
; SI (EGA ITYPEF 2);
VDI1 = SUPTAB.'VDI1';
FINSI;
SI (((DIME MODPLA
) '
>'
0) ET ITHER
); WVMIS = SUPTAB.'ENV1';
FINSI;
si(flmess
); mess 'RECUP DU PAS PRECEDENT OK'
; fins; FINSI;
***************************************************
******** TABLE INFORMATION COMPLEMENTAIRE *********
***************************************************
INFTAB.'MOTTI' = MOTTI;
INFTAB.'MODCOU' = MODCOU;
INFTAB.'TABMOD' = TABMOD;
INFTAB.'LINTER' = LINTER;
INFTAB.'ICOQU' = ICOQU;
INFTAB.'IGDEP' = IGDEP;
INFTAB.'IGDER' = IGDER;
INFTAB.'IREPRI' = IREPRI;
INFTAB.'IPAP' = IPAP;
INFTAB.'ILIN' = ILIN;
INFTAB.'ITHER' = ITHER;
INFTAB.'ITHER1' = ITHER1;
INFTAB.'IPARAL' = IPARAL;
INFTAB.'MATVARI' = MATVARI;
INFTAB.'YOUVARI' = YOUVARI;
INFTAB.'ALFVARI' = ALFVARI;
INFTAB.'IINTE' = IINTE;
INFTAB.'IQUA' = IQUA;
INFTAB.'ELTETA' = ELTETA;
INFTAB.'OBJMOD' = OBJMOD;
INFTAB.'MODPLA' = MODPLA;
INFTAB.'FOR000' = FOR000;
INFTAB.'DEP000' = DEP000;
INFTAB.'CMD000' = CMD000;
INFTAB.'CMD001' = CMD001;
INFTAB.'GRA000' = GRA000;
INFTAB.'IXFEM' = IXFEM;
INFTAB.'ITYPEF' = ITYPEF;
INFTAB.'IPERSO1' = IPERSO1;
* ajout sm
INFTAB.'IDEFI' = IDEFI;
*ajout BP BT pour le contact frottant
INFTAB . 'IFROT' = IFROT ;
si (IFROT);
INFTAB . 'OBJCON' = OBJCON ;
fins;
*fin ajout BP BT
si(IPERSO1); INFTAB . 'ESTIMATION' = ESTIM; fins;
***********************************************
***********************************************
********* BOUCLE SUR LE PAS DE CALCUL *********
***********************************************
***********************************************
REPETER BOUCEXT NBDEP ;
IABC = NBG + &BOUCEXT;
***************************************************
** DEPLACEMENTS,CONTRAINTES ... A L INSTANT INST **
***************************************************
*** SOLUTION_PASAPAS ******************************
SI IPAP;
*** Cas PERSO1
SI IPERSO1;
INST = ESTIM . 'TEMPS' ;
DEPINT
= (ESTIM . 'DEPLACEMENTS'
) REDU ELTETA
; SIGF
= (ESTIM . 'CONTRAINTES'
) REDU OBJMOD
; SI (IGDEP ET (NON IGDER));
SIGF
= '
CAPI' SIGF DEPINT OBJMOD
; FINSI;
SI IGDER;
SI (NON (EXIS (SUPTAB.'ROTATION_RIGIDIFIANTE'
) IABC
)); MESS 'ERREUR : Le deplacement du a une rotation'
; MESS ' rigidifiante au pas ' IABC ' n est pas donne'
; FINSI;
DEPINT = DEPINT -
(REDU SUPTAB.'ROTATION_RIGIDIFIANTE'.
IABC ELTETA
) ; FINSI;
SI (EXIS ESTIM 'VARIABLES_INTERNES'
); VARF
= (ESTIM . 'VARIABLES_INTERNES'
) REDU OBJMOD
; SINON;
VARF = VAR000;
FINSI;
SI (EGA IINTE 2);
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;
DEPINT
= (REDU ELTETA DEPINT
) / DELTAT
; VITDFI
= (REDU ELTETA VITDFI
) / DELTAT
; FINSI;
SI (EGA IINTE 5);
VITF
= (ESTIM .'VITESSES'
) REDU ELTETA
; ACCF
= (ESTIM .'ACCELERATIONS'
) REDU ELTETA
; FINSI;
*** Cas ou on appelle g_theta apres pasapas
SINO;
INST = SUPTAB.'SOLUTION_PASAPAS'.'TEMPS'.IABC ;
DEPINT = (SUPTAB.'SOLUTION_PASAPAS'.'DEPLACEMENTS'.IABC)
SIGF = (SUPTAB.'SOLUTION_PASAPAS'.'CONTRAINTES'.IABC)
SI (IGDEP ET (NON IGDER));
SIGF
= '
CAPI' SIGF DEPINT OBJMOD
; FINSI;
SI IGDER;
SI (NON (EXIS (SUPTAB.'ROTATION_RIGIDIFIANTE'
) IABC
)); MESS 'ERREUR : Le deplacement du a une rotation'
; MESS ' rigidifiante au pas ' IABC ' n est pas donne'
; FINSI;
DEPINT = DEPINT -
(REDU SUPTAB.'ROTATION_RIGIDIFIANTE'.
IABC ELTETA
) ; FINSI;
SI (EXIS (SUPTAB.'SOLUTION_PASAPAS'
) 'VARIABLES_INTERNES'
); VARF = (SUPTAB.'SOLUTION_PASAPAS'.'VARIABLES_INTERNES'.IABC)
SINON;
VARF = VAR000;
FINSI;
SI (EGA IINTE 2);
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;
DEPINT
= (REDU ELTETA DEPINT
) / DELTAT
; VITDFI
= (REDU ELTETA VITDFI
) / DELTAT
; FINSI;
SI (EGA IINTE 5);
VITF = (SUPTAB.'SOLUTION_PASAPAS'.'VITESSES'.IABC)
ACCF = (SUPTAB.'SOLUTION_PASAPAS'.'ACCELERATIONS'.IABC)
FINSI;
FINSI;
*** SOLUTION_RESO **********************************
SINON;
DEPINT
= REDU (SUPTAB.'SOLUTION_RESO'
) ELTETA
; SIGF
= SIGM 'LINE' DEPINT OBJMOD OBJMAT
; FINSI;
*** ON CHANGE LE DEPLACEMENT DEPINT EN MCHAML AU NOEUD ******
*bp: utilite?
DEPINT
= CHAN 'CHAM' DEPINT OBJMOD 'NOEUD' 'DEPLACEMENTS'
;
****************************************************
* MODIFICATION DES CHAMPS THETA ET PI SI GRANDE ROT
****************************************************
SI IPAP;
SI IGDER;
FORM SUPTAB.'ROTATION_RIGIDIFIANTE'.
IABC; SUPTAB.'CHAMP_THETA' UTILTETA
= CH_THETA SUPTAB
; SI (EGA IINTE 4);
SI (NON (EXIS SUPTAB 'FRONT_FISSURE_2'
)); SUPTAB.'COUCHE' = (SUPTAB.'COUCHE') - 1;
SUPTAB.'COUCHE' = (SUPTAB.'COUCHE') + 1;
SINON;
P1 = SUPTAB.'FRONT_FISSURE';
SUPTAB.'FRONT_FISSURE' = SUPTAB.'FRONT_FISSURE_2';
SUPTAB.'FISSURE' = SUPTAB.'FISSURE_2';
SUPTAB.'FRONT_FISSURE' = P1;
FINSI;
FINSI;
FINSI;
FINSI;
***************************************************
********** TEMPERATURES A L INSTANT INST **********
***************************************************
SI ITHER ;
SI IPAP;
SI (ITHER1);
TEPINT
= '
TIRE' CHAR1 INST 'T'
; SINON;
* bp, 2014-11-13 : ajout distinction cas ITHER et ITHER1
TEPINT
= REDU ELTETA SUPTAB.'SOLUTION_PASAPAS'.'TEMPERATURES'.
IABC ; FINSI;
SINON;
TEPINT = SUPTAB.'TEMPERATURES';
FINSI;
FINSI;
***************************************************
********** DEF IMPOSEE A L INSTANT INST **********
***************************************************
SI IDEFI;
SI IPAP;
DEFINT
= '
TIRE' CHAR1 INST 'DEFI'
; SINON;
DEFINT = SUPTAB . 'DEFORMATIONS_IMPOSEES';
FINSI;
FINSI;
***************************************************
********** CONTACT FROTTANT **********
***************************************************
SI IFROT;
SI IPAP;
*...todo
SINON;
* a priori DEPLACEMENT_FISSURE pas tres utile ...
SI (EXIS SUPTAB 'DEPLACEMENT_FISSURE'
); WDEP = SUPTAB . 'DEPLACEMENT_FISSURE';
SINO;
WDEP
= REDU (SUPTAB.'SOLUTION_RESO'
) MAICON
; FINS;
MESS 'ERREUR : IL FAUT DEPLACEMENT_FISSURE
si MODELE_FISSURE'
; FINS;
SI (EXIS SUPTAB 'PRESSION_FISSURE'
); SIGCON = SUPTAB . 'PRESSION_FISSURE';
SINO;
SIGCON
= REDU SIGF OBJCON
; FINS;
* peut etre faire un test sur SIGCON ...
FINSI;
FINS;
*****************************************************
* CONTRAINTE RECALCULEE SI ITHER = VRAI et NON IPAP *
*****************************************************
SI (ITHER ET (NON IPAP));
SIGF
= SIGF
- ('
THET' OBJMOD OBJMAT TEPINT
); FINSI;
***************************************************
************ MATERIAU A L INSTANT INST ************
***************************************************
SI (MATVARI ET ITHER ET IPAP);
SINON;
MAT1 = OBJMAT;
*bp: ne devrait on pas ecrire ci dessous..? inutile?
* MAT1 = REDU OBJMOD OBJMAT
FINSI;
***************************************************
********* RIGIDITE TOTALE A L INSTANT INST ********
***************************************************
*bp: utilite?
SI ( (EGA IINTE 4) OU ((EGA IINTE 99)
et (non IXFEM
) et (EGA &DIME 3) ET (NON IDANS
)) ); SI IPAP;
SI (MATVARI ET ITHER);
(SUPTAB.'SOLUTION_PASAPAS'.'CARACTERISTIQUES');
SINON;
M1 = SUPTAB.'SOLUTION_PASAPAS'.'CARACTERISTIQUES';
FINSI;
RIGTOT
= RIGI M1
(SUPTAB.'SOLUTION_PASAPAS'.'MODELE'
); SINON;
RIGTOT
= RIGI (SUPTAB.'CARACTERISTIQUES'
) (SUPTAB.'MODELE'
); FINSI;
FINSI;
********************************************************************
****** CHARGEMENT MECANIQUE (FORCES NODALES) A L INSTANT INST ******
********************************************************************
PREINT = FOR000;
*** SOLUTION_PASAPAS ******************************
SI IPAP;
PREINT
= PREINT
+ ((TIRE CHAR1 INST 'MECA'
) REDU ELTETA
); FINSI;
SI IMOPRES;
* Dans le cas ou il y a un modele de pression
MAILPTOT
= EXTR MODPRE 'MAIL'
;* --> on isole la partie du MMODEL de pression qui est appliquee
* sur la fissure
MAILPF
= INTE MAILPTOT
(SUPTAB.'FISSURE'
); IPFISS = VRAI;
MODPF
= REDU MODPRE MAILPF
; MATPF
= REDU MATPRE MODPF
; FINSI;
* --> on calcule les forces nodales equivalentes aux pressions
* appliquees hors de la fissure
MAILPEXT
= DIFF MAILPTOT MAILPF
; SI ((NBEL MAILPEXT
) > 0); MODPEXT
= REDU MODPRE MAILPEXT
; MATPEXT
= REDU MATPRE MODPEXT
; PREINT
= PREINT
+ (BSIG MODPEXT MATPEXT
); FINSI;
FINSI;
*** SOLUTION_RESO **********************************
SINON ;
SI (EXIS SUPTAB 'CHARGEMENTS_MECANIQUES'
); PREINT
= PREINT
+ (SUPTAB.'CHARGEMENTS_MECANIQUES'
REDU ELTETA
); FINSI;
SI IMOPRES;
* Dans le cas ou il y a un modele de pression
MAILPTOT
= EXTR MODPRE 'MAIL'
;* --> on isole la partie du MMODEL de pression qui est appliquee sur
* la fissure
MAILPF
= INTE MAILPTOT
(SUPTAB.'FISSURE'
); IPFISS = VRAI;
MODPF
= REDU MODPRE MAILPF
; MATPF
= REDU MATPRE MODPF
; FINSI;
* --> on calcule les forces nodales equivalentes aux pressions
* appliquees hors de la fissure
MAILPEXT
= DIFF MAILPTOT MAILPF
; SI ((NBEL MAILPEXT
) > 0); MODPEXT
= REDU MODPRE MAILPEXT
; MATPEXT
= REDU MATPRE MODPEXT
; PREINT
= PREINT
+ (BSIG MODPEXT MATPEXT
); FINSI;
FINSI;
FINSI;
**********************************************************************
******** POUR LES PRESSIONS APPLIQUEES SUR LA FISSURE
* on calcule le MCHAML vectoriel de [pression * normale]
**********************************************************************
SI IPFISS;
* champ de normale unitaire au maillage de la fissure
NF2
= CHAN 'CHAM' NF1 MODPF 'STRESSES'
; XX
= PSCA NF2 NF2 LCOMP LCOMP
; NF = NF2 / (XX**0.5);
* champ de pression normale au maillage de la fissure
PF
= CHAN 'STRESSES' MODPF MATPF
; PNF = PF * NF LCP LCOMP LCOMP;
FINSI;
****************************************************
*** SOLUTIONS AUXILAIRES SI DECOUPLAGE DES MODES ***
****************************************************
* NBMIXT = nbre d integrale a calculer (=1 si J, =2 si K1 K2, =3
* si K1 K2 K3)
NBMIXT = 1;
IM = 0;
SI (EGA IINTE 99);
NBMIXT = 2;
NBMIXT = 3;
FINSI;
**** CONSTANTES MATERIAUX **************************
CHAM1
= (EXCO 'YOUN' MAT1
) ET (EXCO 'NU ' MAT1
);* -CAS D UN MATERIAU HOMOGENE
SI (NON IDANS);
* on construit les champ auxilaires d'une solution d'un materiau
* homogene => on prend les valeurs de E et nu en pointe de fissure
si(IXFEM);
fins;
SI (EGA (TYPE PM
) 'MAILLAGE'
); SINON;
VYO_1
= EXTR CHPO1 'YOUN' PM
; VNU_1
= EXTR CHPO1 'NU ' PM
; FINSI;
* Constante de Kolosov
SI (EGA &MODE 'PLANCONT'
); KAP_1 = (3. - VNU_1) / (1. + VNU_1);
SINON;
KAP_1 = (3. - (4. * VNU_1));
FINSI;
* Module de cisaillement
MU_1 = VYO_1 / (2.*(1. + VNU_1));
* Constante C_MATE (= 1 / E^etoile)
SI (EGA &MODE 'PLANCONT'
); C_MATE = 1. / VYO_1;
SINON;
C_MATE = (1. - (VNU_1*VNU_1)) / VYO_1;
FINSI;
* -CAS D UN BI-MATERIAU
SINON;
VYO_1
= EXTR CHPO1 'YOUN' PM
; VNU_1
= EXTR CHPO1 'NU ' PM
; VYO_2
= EXTR CHPO1 'YOUN' PM
; VNU_2
= EXTR CHPO1 'NU ' PM
; si(flmess);
fins;
* Constante de Kolosov
SI (EGA &MODE 'PLANCONT'
); KAP_1 = (3. - VNU_1) / (1. + VNU_1);
KAP_2 = (3. - VNU_2) / (1. + VNU_2);
SINON;
KAP_1 = (3. - (4. * VNU_1));
KAP_2 = (3. - (4. * VNU_2));
FINSI;
* Module de cisaillement
MU_1 = VYO_1 / (2.*(1. + VNU_1));
MU_2 = VYO_2 / (2.*(1. + VNU_2));
* Constante bi-metallique EPS1
VA1 = (KAP_1/MU_1) + (1./MU_2);
VA2 = (KAP_2/MU_2) + (1./MU_1);
EPS1 = (1./(2.*VALPI)) * (LOG (VA1/VA2));
* Constante C_MATE (= 1 / E^etoile)
COSH1 = VALPI * EPS1;
COSH1 = ((EXP COSH1) + (EXP (COSH1*(-1.)))) / 2.;
VA1 = (MU_1 + (KAP_1*MU_2)) * (MU_2 + (KAP_2*MU_1));
VA2 = MU_1 * MU_2 * ((MU_1*(1. + KAP_2)) + (MU_2*(1. + KAP_1)));
C_MATE = (COSH1*COSH1*VA1) / (4.*VA2);
mess 'EPS1
=' EPS1 ' C_MATE
=' C_MATE
; FINSI;
* rappel: G = C_MATE * (K1^2 + K2^2) + 1/MU * K3^2
* M = 2*C_MATE * (K1*K1^aux + K2*K2^aux) + 2/MU * K3*K3^aux
* on evite le cas 3d + EF standard
SI (NON ((EGA &DIME 3) ET (NON IXFEM
)));
**** CHAMPS POUR SIMPLIFIER L ECRITURE ************
* (CHPOINT pour EF std // CHAMELEM si XFEM)
* -CAS D UN MATERIAU HOMOGENE
SI (NON IDANS);
SIN1T = sin TETA1; COS1T = cos TETA1;
SIN05T = sin (0.5*TETA1); COS05T = cos (0.5*TETA1);
SI IXFEM;
RM05 = (2.*VALPI*RAY1) ** -0.5 ;
COE_GU = RM05 / (4. * MU_1);
SIN15T = sin (1.5*TETA1); COS15T = cos (1.5*TETA1);
* KAP_1 = KAP_1 * UN1;
SINO;
COE_1 = ((RAY1/(2.*VALPI))**0.5) / (2.*MU_1);
FINS;
* -CAS D UN BI-MATERIAU
SINON;
EPSLGR = EPS1 * (LOG RAY1);
VA1 = COS (EPSLGR*180./VALPI);
VA2 = SIN (EPSLGR*180./VALPI);
BTA1 = ((0.5*VA1) + (EPS1*VA2)) / (0.25 + (EPS1*EPS1));
BTAPM1 = ((0.5*VA2) - (EPS1*VA1)) / (0.25 + (EPS1*EPS1));
DTA_1 = EXP (0. - ((VALPI - TETA1rad)*EPS1));
DTA_2 = EXP ((VALPI + TETA1rad)*EPS1);
GAM_1 = (KAP_1*DTA_1) - (DTA_1**(-1.));
GAM_2 = (KAP_2*DTA_2) - (DTA_2**(-1.));
GAMPM_1 = (KAP_1*DTA_1) + (DTA_1**(-1.));
GAMPM_2 = (KAP_2*DTA_2) + (DTA_2**(-1.));
COS05T = COS (TETA1/2.); SIN05T = SIN (TETA1/2.);
D_1 = (BTA1*GAM_1*COS05T) + (BTAPM1*GAMPM_1*SIN05T);
D_2 = (BTA1*GAM_2*COS05T) + (BTAPM1*GAMPM_2*SIN05T);
DPM_1 = (BTAPM1*GAM_1*COS05T) - (BTA1*GAMPM_1*SIN05T);
DPM_2 = (BTAPM1*GAM_2*COS05T) - (BTA1*GAMPM_2*SIN05T);
GTAR1 = EPSLGR + (0.5*TETA1rad);
CVA_1 = (SIN TETA1) * (SIN (GTAR1*180./VALPI));
CVA_2 = (SIN TETA1) * (COS (GTAR1*180./VALPI));
COE_1 = ((RAY1/(2.*VALPI))**0.5) / (4.*MU_1);
COE_2 = ((RAY1/(2.*VALPI))**0.5) / (4.*MU_2);
FINSI;
fins;
FINSI;
* BOUCLE SUR LES INTEGRALES A CALCULER ==============================*
REPETER BOUCMIX NBMIXT;
IM = IM + 1;
****************************************************
**** CHAMPS AUXILIAIRES SI DECOUPLAGE **************
* Il existe plusieurs manieres de creer les champs aux :
* 1. en utilisant l expression analytique de grad(U) et sigma
* 2. en utilisant l expression analytique de U et en calculant
* sigma, Fint=Bsigma, grad(U) ...
* 3. en appliquant une pression/cisaillement sur les faces de la fissure
* et en resolvant le pb associe
* On utilise la 1ere lorsqu'on peut, la 2eme sinon
SI (EGA IINTE 99);
**** MOTMIX et MOTMIA **************************
SI (IM
EGA 1); MOTMIX
= MOT 'I'
; MOTMIA
= MOT ' I'
; FINSI; SI (IM
EGA 2); MOTMIX
= MOT 'II'
; MOTMIA
= MOT ' II'
; FINSI; SI (IM
EGA 3); MOTMIX
= MOT 'III'
;MOTMIA
= MOT 'III'
; FINSI;
**** METHODE ANALYTIQUE PURE *******************
* -CAS D UN MATERIAU HOMOGENE 2D et 3D --------------------
* SI (IXFEM et (EGA &DIME 2) ET (NON IDANS));
SI (IXFEM ET (NON IDANS));
si(flmess);
mess 'MATERIAU HOMOGENE
XFEM: METHODE ANALYTIQUE'
; fins;
* DERIVEES REPERE CYLINDRIQUE / COORDONNEES GLOBALES
* R,X = DRDX R,Y = DRDY
* T,X = (1/R)*DTDX T,Y = (1/R)*DTDY
DRDX
= ( COS1T
*(EXCO GLV7 'UX,X' 'SCAL'
)) + (SDIR1
*SIN1T
*(EXCO GLV7 'UY,X' 'SCAL'
)) ; DRDY
= ( COS1T
*(EXCO GLV7 'UX,Y' 'SCAL'
)) + (SDIR1
*SIN1T
*(EXCO GLV7 'UY,Y' 'SCAL'
)) ; DTDX
= ( -1.
*SIN1T
*(EXCO GLV7 'UX,X' 'SCAL'
)) + (SDIR1
*COS1T
*(EXCO GLV7 'UY,X' 'SCAL'
)); DTDY
= ( -1.
*SIN1T
*(EXCO GLV7 'UX,Y' 'SCAL'
)) + (SDIR1
*COS1T
*(EXCO GLV7 'UY,Y' 'SCAL'
)); DRDZ
= ( COS1T
*(EXCO GLV7 'UX,Z' 'SCAL'
)) + (SDIR1
*SIN1T
*(EXCO GLV7 'UY,Z' 'SCAL'
)) ; DTDZ
= ( -1.
*SIN1T
*(EXCO GLV7 'UX,Z' 'SCAL'
)) + (SDIR1
*COS1T
*(EXCO GLV7 'UY,Z' 'SCAL'
)); fins;
* -debut du cas contact frottant IFROT (btrolle 19/02/2013)
* ajout des solutions analytiques du saut sur les levres de la fissure
* projection de psi1 sur la fissure et calcul de psi,x
'SI' IFROT;
* PSI en CHAML aux PG sur la fissure utilise pour calcul du gradient
* on selectionne la partie de psi2 non nulle (= dans le champ theta)
* contour du domaine theta
con1 = 'CONTOUR' ELTETA;
* maillage support pour l'integration
mai1 = 'INCLUSION' MAICON con1 'STRI';
* OBJCON2 = MODE mai1 'MECANIQUE' 'ZCO2';
OBJCON2
= REDU OBJCON mai1
; 'SINON';
* maillage support pour l'integration
mai1
= 'INCLUSION'
('EXTRAIRE' OBJCON 'MAIL'
) ELTETA '
VOLU'
'STRI';
* OBJCON2 = 'MODELISER' mai1 'MECANIQUE' 'ZCO3';
OBJCON2
= REDU OBJCON mai1
; 'FINSI';
PSI1e
= CHAN 'CHAM' PSI1 OBJMOD 'NOEUDS'
; PSI2
= PROI OBJCON2 PSI1e 'STRESSES'
;* PSI en CHPOINT sur la fissure utilise pour calcul du repere local
PSI3
= PROI MAICON PSI1e
; TesPsi = PSI2 'MASQUE' 'INFERIEUR'(-1E-15) ;
* PSI2B = 'CHANGER' 'CHAM' (TesPsi * (-1.*PSI2) ('MOTS' 'PSI') ('MOTS'
* ('MOTS' 'PSI')) OBJCON 'STRESSES';
PSI2B
= (TesPsi
* (-1.
*PSI2
) ('
MOTS' 'PSI'
) ('
MOTS' 'PSI'
)* 'MESSAGE' ' PSI2B = '; 'LISTE' PSI2B;
* terme sqrt(r/2pI) et sa derivee,r
RM05B = ((PSI2B/(2.*VALPI)) **0.5);
RM05BR =0.5*((2.*VALPI*PSI2B)**-0.5);
* Change le nom pour pouvoir calculer grad avec ZCO
et (MANU '
CHPO'
('EXTRAIRE' OBJCON2 'MAIL'
) 1 'AY '
0.
'NATURE' 'DIFFUS')
et (MANU '
CHPO'
('EXTRAIRE' OBJCON2 'MAIL'
) 1 'AZ '
0.
'NATURE' 'DIFFUS');
GLV72
= 'CHANGER'
(GRAD LV72 OBJCON2
) '
TYPE' 'SCALAIRE'
;* 'MESSAGE' ' GLV72 =';'LISTE' GLV72 ;
* Angle ALPHA1 de passage local -> global
NGPSI2
= ( ((EXCO GLV72 'AX,X' 'SCAL'
)**2) + ((EXCO GLV72 'AX,Y' 'SCAL'
)**2) )**(0.5) ; COS1AB
= ( (EXCO GLV72 'AX,Y' 'SCAL'
) / NGPSI2
) ; SIN1AB
= ( (EXCO GLV72 'AX,X' 'SCAL'
) / NGPSI2
) ; 'SINON';
* Angle ALPHA1 de passage local -> global
NGPSI2
= ( ((EXCO GLV72 'AX,X' 'SCAL'
)**2) + ((EXCO GLV72 'AX,Y' 'SCAL'
)**2) + ((EXCO GLV72 'AX,Z' 'SCAL'
)**2) )**(0.5) ; * CHAM premiere tangente = grad de psi
V12X
= ('
NOMC' 'UX'
((EXCO GLV72 'AX,X' 'SCAL'
) '
/' NGPSI2
)); V12Y
= ('
NOMC' 'UY'
((EXCO GLV72 'AX,Y' 'SCAL'
) '
/' NGPSI2
)); V12Z
= ('
NOMC' 'UZ'
((EXCO GLV72 'AX,Z' 'SCAL'
) '
/' NGPSI2
)); V12 = V12X 'ET' V12Y 'ET' V12Z;
* CHAM des normales
V22Y
