5. Mode opératoire : les routines à écrire#
C’est à ce niveau que doit être fait le choix du type d’intégration. Il existe quatre possibilités:
Utiliser l’architecture de l’environnement d’intégration explicite par un schéma de Runge-Kutta d’ordre 2 [R5.03.14] (ALGO_INTE=”RUNGE_KUTA” ):
il s’agit de la méthode la plus simple. Outre la récupération des données matériaux, il suffit d’écrire une routine calculant les dérivées des variables internes
le calcul de l’opérateur tangent n’est pas disponible sous cet environnement, c’est l’opérateur de rigidité élastique qui est utilisé
Implantation «complète» du nouveau comportement dans l’environnement d’intégration implicite PLASTI [R5.03.14] ( ALGO_INTE=”NEWTON” ):
résolution du système non linéaire local par la méthode de Newton. Outre la récupération des données matériaux, il faut écrire plusieurs routines appelées dans l’algorithme de Newton local (évaluation du seuil, calcul du résidu, calcul analytique de la matrice jacobienne…)
L’opérateur cohérent est obtenu directement à partir de la jacobienne du système local, et l’opérateur tangent en prédiction est par défaut l’opérateur de rigidité élastique
PLASTI ne permet pas d’obtenir des modèles optimisés en temps calcul
Implantation «facile» du nouveau comportement dans l’environnement d’intégration implicite PLASTI avec [R5.03.14] ( ALGO_INTE=”NEWTON_PERT” ):
le système non-linéaire local peut être réécrit de telle sorte que l’évaluation du résidu ne nécessite de la part du développeur que de spécifier l’expression des dérivées des variables internes écrites dans le cadre de la méthode 2 (intégration explicite par RK2). On peut donc également réaliser une intégration implicite dans PLASTI avec les deux seules routines nécessaires à l’intégration explicite.
La jacobienne du système local est calculée par perturbation, le calcul est donc encore plus coûteux qu’avec la méthode 2. De la même manière, l’opérateur tangent cohérent est obtenu directement à partir de la jacobienne du système local
Créer une routine autonome d’intégration complète du comportement:
permet souvent d’obtenir les modèles les plus performants (par exemple, en réduisant le système à résoudre à une seule équation scalaire, non linéaire, voir par exemple [R5.03.04], [R5.03.16], [R5.03.21], …)
nécessite plus de travail « sur le papier » pour optimiser les équations
Attention:
Dans le cas 4, on choisira un numéro de routine nn et on écrira la routine lc00nn. Dans les autres cas on choisira comme point d’entrée le numéro 32: LC0032 appelle PLASTI ou NMVPRK (Runge-Kutta) suivant la valeur de ALGO_INTE choisie par l’utilisateur.
5.1. Première possibilité : introduire un nouveau comportement explicite – schéma de RUNGE - KUTTA#
Ce type d’intégration correspond à ALGO_INTE=”RUNGE_KUTTA”, c’est la façon la plus rapide d’introduire un nouveau comportement.
Il faut écrire dans un premier temps une routine XXXMATappelée par la routine d’aiguillage LCMATEafin de récupérer les paramètres matériaux et la taille du système différentiel non-linéaire à intégrer.
SUBROUTINE XXXMAT(FAMI,KPG,KSP,MOD,IMAT,NMAT,MATERD,MATERF,
MATCST,NDT,NDI,NR,NVI,VIND)
Arguments en entrée :
FAMI,KPG,KSP : famille et numéro de point de gauss / sous-point
IMAT : adresse du matériau
MOD : type de modelisation
NMAT : dimension de MATERD / MATERF
Arguments en sortie :
MATERD : coefficients matériau a t
MATERF : coefficients matériau a t+dt
MATERx(*,1) = caractéristiques élastiques
MATERx(*,2) = caractéristiques plastiques
MATCST : 'oui' si matériau a t = matériau a t+dt
'non' sinon
NDT : nb total de composantes des tenseurs
NDI : nb de composantes directes des tenseurs
NR : nb de composantes du système non-linéaire
NVI : nb de variables internes
On donne ci-dessous un exemple pour chacune des deux fonctions principales que doit remplir cette routine
AFFECTATION DES DIMENSIONS DU PROBLEME LOCAL (NDT,NDI,NR,NVI)
NVI=7
IF ( MOD .EQ.'3D' ) THEN
NDT = 6
NDI = 3
NR = NDT+2
ELSE IF ( MOD .EQ. 'D_PLAN' .OR. MOD .EQ. 'AXIS') THEN
NDT = 4
NDI = 3
NR = NDT+2
ELSE
CALL U2MESS('F',...)
ENDIF
RECUPERATION DU MATERIAU
NOMC(1) = 'E'
NOMC(2) = 'NU'
NOMC(3) = 'ALPHA'
NOMC(4) = 'SY'
NOMC(5) = 'D_SIGM_EPSI'
CALL RCVALB(FAMI,KPG,KSP,'-',IMAT,' ','ELAS',0,' ',
& 0.D0,3,NOMC(1),MATERD(1,1),ICODRE,1)
CALL RCVALB(FAMI,KPG,KSP,'-',IMAT,' ','ECRO_LINE',0,' ',
& 0.D0,2,NOMC(4),MATERD(1,2),ICODRE,1)
Il faut ensuite écrire une routine RKDXXX appelée par la routine d’aiguillage LCDVIN et donnant les dérivées temporelles des variables internes.
Exemples de routine RKDXXX: RKDCHA, RKDVEC, RKDHAY.
5.2. Deuxième possibilité : introduction « complète » d’un nouveau comportement dans PLASTI (implicite)#
Cet type d’intégration correspond à ALGO_INTE=”NEWTON”. L’environnement PLASTI permet d’intégrer de manière systématique des relations de comportement non-linéaires par une méthode de Newton locale (au niveau du point de Gauss). Connaissant les contraintes et les variables internes à l’instant \(i-1\) ainsi que l’incrément de déformation totale \(\Delta {\epsilon}_{i}^{n}\) donné par l’algorithme de Newton global, le système local d’équations à résoudre sous forme purement implicite est écrit de la façon suivante:
\(R(\Delta y)=(\begin{array}{c}g(\Delta y)\\ l(\Delta y)\\ f(\Delta y)\end{array})=0\) avec \(\Delta y=(\begin{array}{c}\Delta \sigma \\ \Delta \text{vari}\\ \Delta p\end{array})\)
La première équation représente par exemple la relation contrainte-déformation élastique (6 équations à 6 inconnues), avec \(\mathrm{A}\) l’opérateur d’élasticité (éventuellement modifié pour les lois avec endommagement), \(\Delta {\varepsilon}^{p}\) la variation de déformation plastique et \(\Delta {\varepsilon}^{\mathit{th}}\) la variation de déformation thermique:
\(g(\Delta y)=\Delta \sigma -\mathrm{A}(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{\mathit{th}}-\Delta {\varepsilon}^{p})=0\)
la seconde représente l’ensemble des lois d’évolution des différentes variables internes scalaires et/ou vectorielles (\({n}_{v}\) équations scalaires à \({n}_{v}\) inconnues)
\(l(\Delta y)=0\)
la dernière représente le critère éventuel de plasticité (1 équation)
\(f(\Delta y)=0\)
Ce système de \(6+{n}_{v}(+1)\) équations à \(6+{n}_{v}(+1)\) inconnues est résolu par une méthode de Newton:
\(\lbrace \begin{array}{c}\frac{\partial R}{\partial \Delta y}(\Delta {y}_{k}).d[\Delta {y}_{k}]=-R(\Delta {y}_{k})\\ \Delta {y}_{k+1}=\Delta {y}_{k}+d(\Delta {y}_{k})\end{array}\)
A convergence,on obtient donc les incréments de contraintes et de variables internes. L’opérateur tangent cohérent est quant à lui calculé de manière systématique à partir de la jacobienne du système local par la routine LCOPTG (voir [R5.03.14] pour le détail des équations). Il faut donc programmer a minima une routine définissant le résidu (formé des équations ci-dessus) ainsi qu’une routine construisant la matrice jacobienne. On décrit brièvement ci-dessous l’architecture générale de PLASTI, en indiquant au fur et à mesure la liste des routines à écrire.
Architecture générale de PLASTI:
CALL LCMATE(...)
→écriture nécessaired'une routine spécifique XXXMAT de récupération du matériau identique à RUNGE_KUTTA
IF ( OPT .EQ. 'RAPH_MECA' .OR. OPT .EQ. 'FULL_MECA' ) THEN
INTEGRATION ELASTIQUE SUR DT
CALL LCELAS(...)
→écriture éventuelle d'une routine spécifique XXXELA (par défaut LCELIN: élasticité linéaire)
PREDICTION ETAT ELASTIQUE A T+DT
CALL LCCNVX(...,SEUIL)
→écriture nécessaired'une routine spécifique XXXCVX d'évaluation du seuil
IF ( SEUIL .GE. 0.D0 ) THEN
CALL LCPLAS(...)
ENDIF
ENDIF
La routine LCPLAS appelle LCPLNL, qui réalise la boucle de Newton dont la structure est la suivante:
Notations:
YD=(SIGD,VIND): vecteur des inconnues (de dimension \(6+{n}_{v}\) ) à l’instant T
YF=(SIGF,VINF): vecteur des inconnues à l’instant T+DT
DY: incrément du vecteur des inconnues entre les instants Tet T+DT
DDY: incrément de vecteur des inconnues entre deux itérations de Newton successives
R:résidu
DRDY:jacobienne
On résout donc:(DY) = 0
Par une méthode de Newton DRDY(DYK) DDYK = - R(DYK)
DYK+1 = DYK + DDYK (DY0 DEBUT)
et on réactualise YF = YD + DY
CALCUL DE LA SOLUTION D ESSAI INITIALE DU SYSTEME NL EN DY
CALL LCINIT(...DY,...)
→écriture éventuelled'une routine spécifique XXXINI (par défaut DY est initialisé à 0)
ITERATIONS DE NEWTON
ITER = 01 CONTINUE
ITER = ITER + 1
INCREMENTATION DE YF = YD + DY
CALL LCSOVN(NR,YD,DY,YF)
CALCUL DES TERMES DU SYSTEME A T+DT = -R(DY)
CALL LCRESI(...,DY,R,IRET)
→écriture nécessaired'une routine spécifique XXXRES de calcul du résidu
CALCUL DU JACOBIEN DU SYSTEME A T+DT = DRDY(DY)
si ALGO_INTE='NEWTON' calcul exact de la jacobienne
CALL LCJACB(...DY,...DRDY,IRET)
→écriture nécessaire d'une routine spécifique XXXJAC
sinon si ALGO_INTE='NEWTON_PERT', calcul par perturbation (cf §:ref:`5.3 <RefNumPara__3809_1482954243>` )
CALL LCJACP(...DRDY,...)
RESOLUTION DU SYSTEME LINEAIRE DRDY(DY).DDY = -R(DY)
CALL LCEQMN(NR,DRDY,DRDY1)
CALL LCEQVN(NR, R , DDY)
CALL MGAUSS('NCWP',DRDY1,DDY,NR,NR,1,RBID,IRET )
REACTUALISATION DE DY = DY + DDY
CALL LCSOVN ( NR , DDY , DY , DY )
RECHERCHE LINEAIRE dans le cas ALGO_INTE='NEWTON_RELI'
CALL LCRELI ( …)
ESTIMATION DE LA CONVERGENCE
CALL LCCONV(DY,DDY,NR,ITMAX,TOLER,...,R,...,IRTET)
→écriture éventuelled'une routine spécifique XXXCVG du critère de convergence (critère relatif par défaut dans LCCONG)
IF ( IRTET.GT.0 ) GOTO 1
CONVERGENCE -> INCREMENTATION DE YF = YD + DY
CALL LCSOVN ( NDT+NVI , YD , DY , YF )
MISE A JOUR DE SIGF , VINF
CALL LCEQVN ( NDT , YF(1) , SIGF )
CALL LCEQVN ( NVI-1 , YF(NDT+1) , VINF )
En résumé, pour une introduction «complète» d’un nouveau comportement dans PLASTI, il faut nécessairement écrire les routines spécifiques suivantes:
XXXMAT appelée par LCMATE: récupération du matériau et de la taille du problème local
XXXCVX appelée par LCCNVX: évaluation du seuil
XXXRES appelée par LCRESI: calcul du résidu
XXXJAC appelée par LCJACB: calcul de la jacobienne
Il peut aussi être utile, selon le besoin, d’écrire les routines spécifiques suivantes:
XXXELA appelée par LCELAS: intégration élastique (si élasticité non-linéaire)
XXXINI appelée par LCINIT: initialisation (pour une initialisation autre que DY0=0)
XXXCVG appelée par LCCONV: pour modifier le critère de convergence
5.3. Troisième possibilité : utilisation des routines de l’intégration explicite dans une intégration implicite avec PLASTI#
5.3.1. Principe#
Ce dernier cas correspond à ALGO_INTE_=”NEWTON_PERT” , il s’agit de la méthode d’implantation «facile» du nouveau comportement dans l’environnement d’intégration implicite PLASTI. Il est possible d’utiliser directement les deux routines XXXMAT et RKDXXX (récupération des données matériau, et dérivées des variables internes) utilisées avec ALGO_INTE_=”RUNGE_KUTTA” pour réaliser un intégration implicite. En effet, le système d’équations différentielles résolu par RUNGE_KUTTA peut s’écrire :
\(\lbrace \begin{array}{c}\Delta \sigma =\mathrm{A}(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{\mathit{th}}-\Delta {\varepsilon}^{p}(Y))\\ \frac{dY}{dt}=F(Y,t;\sigma )\end{array}\)
où \(Y\) représente l’ensemble des variables internes du modèle. La relation entre le tenseur des contraintes et la partie élastique du tenseur des déformations est généralement linéaire, mais peut être évaluée de façon non linéaire par une expression spécifique.
Une fois programmée la routine RKDXXX permettant de calculer \(\frac{dY}{dt}=F(Y,t;\sigma )\) , il est possible de l’utiliser pour une intégration implicite, ce qui consiste à résoudre ( cf. [R5.03.14]) :
\(R(\Delta Z)=0=\left[\begin{array}{c}{R}_{1}(\Delta Z)\\ {R}_{2}(\Delta Z)\end{array}\right]\) , avec \(\Delta Z=(\begin{array}{c}\Delta \sigma \\ \Delta Y\end{array})=Z(t+\Delta t)-Z(t)\)
Le premier système d’équations représente la relation contrainte - déformation élastique
\({R}_{1}(\Delta Z)={A}^{-1}\sigma -(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{\mathit{th}}-\Delta {\varepsilon}^{p}(Y))={A}^{-1}\sigma -G(Y)=0\)
Par convention, les premières valeurs de \(Y\) représentent la variation de déformation plastique, pour faciliter le calcul de \(G(Y)\) (voir la routine LCRESA pour plus de détails)
Le deuxième exprime les lois d’évolution des différentes variables internes, soit après discrétisation temporelle par un schéma d’Euler implicite:
\({R}_{2}(\Delta Z)=\Delta Y-\Delta t.F(Y,\sigma )=0\)
Ce système est résolu par la méthode de Newton proposée dans l’environnement PLASTI et décrite au paragraphe précédent:
\(\lbrace \begin{array}{c}\frac{\partial R}{\partial \Delta Z}d(\Delta {Z}_{k})\text{= -}R(\Delta {Z}_{k})\\ \Delta {Z}_{k+1}=\Delta {Z}_{k}+d(\Delta {Z}_{k})\end{array}\)
Les quantités \(G\) et \(F\) intervenant dans le résidu sont calculées par la routine «explicite» RKDXXX à écrire, et le résidu est construit automatiquement par la routine LCRESA. Le matrice jacobienne est calculée automatiquement par perturbation (routine LCJACP). L’opérateur tangent cohérent est quant à lui calculé de manière systématique à partir de la jacobienne (routine LCOPTG, voir [R5.03.14] pour le détail des équations).
En résumé, ce procédé permet, avec les deux seules routines nécessaires à l’intégration explicite, (coefficients matériau et calcul des dérivées des variables internes) d’utiliser une intégration implicite, et de bénéficier d’une matrice tangente. Ce procédé est économique en termes de temps de développement, mais a priori moins efficace en temps CPU qu’une matrice jacobienne programmée explicitement.
On détaille ci-dessous le calcul de l’opérateur tangent par perturbation, ainsi que le critère de convergence de l’algorithme de Newton local.
Calcul par perturbation (LCJACP): différences finies d’ordre 2
Initialisation de la perturbation: \(\eta ={10}^{-7}∣∣\Delta Z∣∣\)
boucle sur les colonnes \(j\) de la matrice à remplir:
calcul de \(R(\Delta Z+\eta {I}_{j})\) avec \({I}_{j}={[00\dots 1\dots 00]}^{T}\) vecteur nul sauf à la ligne \(j\)
calcul de \(R(\Delta Z-\eta {I}_{j})\)
calcul de la colonne \(j\) : \({\left[\frac{\partial R}{\partial \Delta Z}\right]}_{\text{...},j}\mathrm{\simeq }\frac{R(\Delta Z+\eta {I}_{j})-R(\Delta Z-\eta {I}_{j})}{2\eta }\)
Critère de convergence (LCCONG):
On sépare les deux blocs \({R}_{1}\) et \({R}_{2}\) du résidu pour éviter les problèmes dus à des ordres de grandeur différents. A chaque itération \(k\) de l’algorithme de Newton local, on calcule:
\({\mathit{err}}_{1}=\frac{{∣∣{R}_{1}(\Delta {Z}_{k})∣∣}_{\infty}}{{∣∣{R}_{1}(\Delta {Z}_{0})∣∣}_{\infty}}\) , avec \({R}_{1}(\Delta {Z}_{0})=\Delta {\varepsilon}_{i}^{n}\) car on a \(\Delta {Z}_{0}=0\) à l’initialisation. (cf. § 5.2 )
\({\mathit{err}}_{2}=\frac{{∣∣{R}_{2}(\Delta {Z}_{k})∣∣}_{\infty}}{{∣∣Y(t)+\Delta {Y}_{k}∣∣}_{\infty}}\)
Le critère d’arrêt est alors le suivant:
\(max({\mathit{err}}_{1,}{\mathit{err}}_{2})<\xi\) , où \(\xi\) est donné par RESI_INTE_RELA.
5.3.2. Exemple#
Un exemple: la loi visco-élasto-plastique de Hayhurst.
Routine de lecture des coefficients matériau (appelée par la routine d’aiguillage LCMATT):
HAYMAT
Routine de calcul des dérivées des variables internes (appelée par la routine d’aiguillage LCDVIN):
RKDHAY
Les développements explicite / implicite de cette relation de comportement sont testés et comparés dans le cas-test ssnv225 [V6.04.225]. Il s’agit d’un test au point matériel de fluage en grandes déformations permettant de valider les capacités du modèle de HAYHURST à représenter le fluage primaire, secondaire et tertiaire. Voici les caractéristiques d’exécution pour ce test en version 11.2 :
Modélisation A: ALGO_INTE=”RUNGE_KUTTA”
temps CPU: 90.59 s
Nombre de pas de temps: 1660
Nombre d’itérations de Newton: 7615
Modélisation B: ALGO_INTE=”NEWTON_PERT”:
temps CPU: 27.38s
Nombre de pas de temps: 520
Nombre d’itérations de Newton: 1404
Les résultats sont quasi identiques :
5.4. Quatrième possibilité : écrire une routine lc00nn autonome#
5.4.1. lc00nn : routine relative à un point d’intégration d’un élément, spécifique à une loi de comportement#
Chercher dans le répertoire bibfor/algorith un numéro de routine lc00nn non utilisé \((50<\mathit{nn}<100)\) ) et partir de cette routine vide.
Remarque:
L’appel à lc00nn par lc0000 (qui est la routine appelant toute les routines d’intégration des comportement disponibles dans code_aster) est déjà écrit. Cependant, il faut s’assurer que les arguments (ainsi que l’ordre de ces arguments) choisis à la déclaration de lc00nn par le développeur soient les mêmes qu’à l’appel dans lc0000 .
Les arguments d’entrée d’une routine lc00nn sont a minima :
FAMI famille de points de gauss (RIGI,MASS,…)
KPG,KSP numéro du point de gauss et du sous-point
NDIM dimension de l espace (\(\mathrm{3d}=3\) , \(\mathrm{2d}=2\) , \(\mathrm{1d}=1\) )
IMATE adresse du matériau
COMPOR infos sur le comportement
compor(1) = relation de comportement (
vmis_cine_…)compor(2) = nombre de variables internes
compor(3) = type de déformation (petit,green…)
CRIT critères locaux
crit(1) = nombre d itérations maxi a convergence (ITER_INTE_MAXI)
crit(3) = valeur de la tolérance de convergence (RESI_INTE_RELA)
INSTAM instant t-
INSTAP instant t = t- + dt
EPSM déformation totale a t- (ou éventuellement le gradient de transformation suivant le type de déformation : les arguments de la routine appelante, lc0000, contiennent la dimension de EPSM, et DEPS),
DEPS incrément de déformation totale (même remarque)
SIGM contrainte a t-
VIM variables internes a t-
OPTION option de calcul
RIGI_MECA_TANG -> dsidep(t)
FULL_MECA -> dsidep(t+dt) , sig(t+dt)
RAPH_MECA -> sig(t+dt)
ANGMAS les trois angles (nautiques ou d’Euler) du mot_clef massif
TYPMOD type de modélisation : 3D, AXIS, D_PLAN, C_PLAN,…
ICOMP compteur pour le redécoupage local du pas de temps
NVI nombre de variables internes du comportement
les arguments de sortie sont, selon l’option de calcul:
VIP variables internes a l’instant actuel (options RAPH_MECA et FULL_MECA)
SIGP terme d’ordre 0 dans la linéarisation de la contrainte (option RIGI_MECA_TANG) – cf. § 4
contraintes a l’instant actuel (options RAPH_MECA et FULL_MECA)
DSIDEP opérateur tangent cohérent ou en vitesse (option FULL_MECA ou RIGI_MECA_TANG).
CODRET code retour permettant d’indiquer (s’il est non nul) une problème d’intégration
locale, donc d’effectuer un redécoupage du pas de temps
Remarques :
Le cas échéant, il est possible d’utiliser également un type dérivé en entrée ( variable BEHinteg dans LC0000/LCnnnn ). Ce type contient des arguments supplémentaires, par exemple une longueur caractéristique dans le cas des modèles non locaux…
De même, il est possible de transférer des arguments supplémentaires en sortie de la routine LCnnnn/LC0000 via la variable BEHinteg.
Attention, dans le cas RIGI_MECA_TANG, les tableaux relatifs aux contraintes et variables internes en fin de pas de temps ne sont pas alloués. Il ne faut donc pas les utiliser pour calculer la matrice tangente de prédiction.
5.4.2. Organisation de la routine à écrire#
On prend comme au § 3.3 l’exemple de la loi de Von Mises à écrouissage cinématique linéaire VMIS_CINE_LINE (num_lc=3 dans vmis_cine_line.py):
SUBROUTINE LC0003(FAMI,KPG,KSP,NDIM,IMATE,COMPOR,CRIT,INSTAM,
& INSTAP,EPSM,DEPS,SIGM,VIM,OPTION,ANGMAS,SIGP,VIP,
& TAMPON,TYPMOD,ICOMP,NVI,DSIDEP,CODRET)
Cette routine est en fait une routine d’aiguillage pour les comportements VMIS_CINE_LINE et VMIS_ECMI_* . L’intégration de VMIS_CINE_LINE est réalisée dans la routine NMCINE, appelée par lc0003 et dont le contenu est décrit ici brièvement.
On notera que dans cet exemple, la loi suit le mécanisme deform_ldc = “OLD”, et qu’à ce titre, les déformations thermiques sont calculées et retranchées à l’intérieur de la loi (cf. § 4 ).
LECTURE DES CARACTERISTIQUES ELASTIQUES DU MATERIAU (TEMPS T = +)
NOMRES(1)='E'
NOMRES(2)='NU'
CALL RCVALB(FAMI,KPG,KSP,'+',IMATE,' ','ELAS',0,' ',
& 0.D0,2,NOMRES,VALRES,ICODRE,2)
E = VALRES(1)
NU = VALRES(2)
Remarque :
RCVALB est une routine générale permettant d’interpoler les valeurs des coefficients matériaux par rapport aux variables de commandes dont ils dépendent (voir les utilitaires).
RCVARC est une routine qui permet de récupérer la valeur d’une variable de commande (température, séchage , irradiation, …) à l’instant considéré, et au point de Gauss considéré (voir les utilitaires). Exemple:
CALL RCVARC(” “,”TEMP”,”-“,FAMI,KPG,KSP,TM,IRET)
CALL RCVARC(” “,”TEMP”,”+”,FAMI,KPG,KSP,TP,IRET)
LECTURE DES CARACTERISTIQUES D’ECROUISSAGE
NOMRES(1)='D_SIGM_EPSI'
NOMRES(2)='SY'
CALL RCVALB(FAMI,KPG,KSP,'+',IMATE,' ','ECRO_LINE',0,' ',
& 0.D0,2,NOMRES,VALRES,ICODRE, 2)
DSDE=VALRES(1)
SIGY=VALRES(2)
C = 2.D0/3.D0*DSDE/(1.D0-DSDE/E)
CALCUL DES CONTRAINTES ELASTIQUES ET DU CRITERE DE VON MISES
DO 110 K=1,3
DEPSTH(K) = DEPS(K) -EPSTHE
DEPSTH(K+3) = DEPS(K+3)
110 CONTINUE
EPSMO = (DEPSTH(1)+DEPSTH(2)+DEPSTH(3))/3.D0
DO 115 K=1,NDIMSI
DEPSDV(K) = DEPSTH(K) - EPSMO * KRON(K)
115 CONTINUE
SIGMO = (SIGM(1)+SIGM(2)+SIGM(3))/3.D0
SIELEQ = 0.D0
DO 114 K=1,NDIMSI
SIGDV(K) = SIGM(K) - SIGMO*KRON(K)
SIGDV(K) = DEUXMU/DEUMUM*SIGDV(K)
SIGEL(K) = SIGDV(K) + DEUXMU * DEPSDV(K)
SIELEQ = SIELEQ + (SIGEL(K)-C/CM*VIM(K))**2
114 CONTINUE
SIGMO = TROISK/TROIKM * SIGMO
SIELEQ = SQRT(1.5D0*SIELEQ)
SEUIL = SIELEQ – SIGY
DP = 0.D0
PLASTI=VIM(7)
CALCUL DES CONTRAINTES ET DES VARIABLES INTERNES
Les expressions des contraintes et des variables internes (RAPH_MECA et FULL_MECA) sont données dans [R5.03.02]
IF ( OPTION(1:9) .EQ. 'RAPH_MECA' .OR.
& OPTION(1:9) .EQ. 'FULL_MECA' ) THEN
IF (SEUIL.LT.0.D0) THEN
VIP(7) = 0.D0
DP = 0.D0
SIELEQ = 1.D0
A1 = 0.D0
A2 = 0.D0
ELSE
VIP(7) = 1.D0
DP = SEUIL/(1.5D0*(DEUXMU+C))
A1 = (DEUXMU/(DEUXMU+C))*(SEUIL/SIELEQ)
A2 = (C /(DEUXMU+C))*(SEUIL/SIELEQ)
ENDIF
PLASTI=VIP(7)
DO 160 K = 1,NDIMSI
SIGDV(K) = SIGEL(K) - A1*(SIGEL(K)-VIM(K)*C/CM)
SIGP(K) = SIGDV(K) + (SIGMO + TROISK*EPSMO)*KRON(K)
VIP(K) = VIM(K)*C/CM + A2*(SIGEL(K)-VIM(K)*C/CM)
160 CONTINUE
ENDIF
CALCUL DE L’OPERATEUR TANGENT “DSIPSEP”: RIGI_MECA_TANG(en vitesse) ou FULL_MECA(cohérent)
IF ( OPTION(1:14) .EQ. 'RIGI_MECA_TANG' .OR.
& OPTION(1:9) .EQ. 'FULL_MECA' ) THEN
CALL MATINI(6,6,0.D0,DSIDEP)
DO 120 K=1,6
DSIDEP(K,K) = DEUXMU
120 CONTINUE
IF ( OPTION(1:14) .EQ. 'RIGI_MECA_TANG') THEN
DO 174 K = 1,NDIMSI
SIGDV(K) = SIGDV(K) - VIM(K)*C/CM
174 CONTINUE
ELSE
DO 175 K = 1,NDIMSI
SIGDV(K) = SIGDV(K) - VIP(K)
175 CONTINUE
ENDIF
SIGEPS = 0.D0
DO 170 K = 1,NDIMSI
SIGEPS = SIGEPS + SIGDV(K)*DEPSDV(K)
170 CONTINUE
A1 = 1.D0/(1.D0+1.5D0*(DEUXMU+C)*DP/SIGY)
A2 = (1.D0+1.5D0*C*DP/SIGY)*A1
IF(PLASTI.GE.0.5D0.AND.SIGEPS.GE.0.D0) THEN
COEF = -1.5D0*(DEUXMU/SIGY)**2 / (DEUXMU+C) * A1
DO 135 K=1,NDIMSI
DO 135 L=1,NDIMSI
DSIDEP(K,L) = A2 * DSIDEP(K,L) + COEF*SIGDV(K)*SIGDV(L)
135 CONTINUE
LAMBDA = LAMBDA + DEUXMU**2*A1*DP/SIGY/2.D0
ENDIF
DO 130 K=1,3
DO 131 L=1,3
DSIDEP(K,L) = DSIDEP(K,L) + LAMBDA
131 CONTINUE
130 CONTINUE
ENDIF
5.4.3. Variables de commande (external state variables)#
La dépendance des paramètres matériau aux variables de commande (séchage, température, etc.) est prise en compte par l’usage de la routine RCVALB. Il existe en sus des variables de commandes spécifiques à certaines lois de comportement qui doivent être calculées (en général au niveau de l’élément) et non transmises par le mécanisme AFFE_MATERIAU/AFFE_VARC, il s’agit de:
ELTSIZE1: calcul de la taille d’un élément (pour la loi BETON_DOUBLE_DP);
ELTSIZE2: calcul de la taille d’un élément (pour la loi ENDO_PORO_BETON);
COORGA: coordonnées des point de Gauss de l’élément (pour la loi META_LEMA_ANI);
GRADVELO: pour le gradient de la vitesse de déformation (pour la loi MONOCRISTAL);
HYGR: hygrométrie (à partir de la fonction de désorption FONC_DESORP).
Les lois de comportement utilisant ces variables utilisent le mot-clef exte_vari dans leur catalogue. L’ajout d’une autre variable de ce style nécessite un développement spécifique qui ne peut donc pas être détaillé ici.
Par contre, on peut utiliser les variables définies plus haut. Il faut:
impacter le catalogue de la loi de comportement et ajouter cette variable dans exte_vari;
récupérer la valeur de cette variable via le module calcul
Variables |
Variable (s) dans le module calcul |
HYGR |
|
ELTSIZE1 |
|
ELTSIZE2 |
ca_vext_eltsize2_(9) |
COORGA |
ca_vext_coorga_(27,3) |
GRADVELO |
ca_vext_gradvelo_(9) |
Par exemple:
use calcul_module, only : ca_vext_gradvelo_
real (kind=8):: l(3,3)
integer:: i, j
do i = 1, 3
do j = 1, 3
l(i,j)=ca_vext_gradvelo_(3*(i-1)+j)
end do
end do