Download ccreep.eso
Back to the list
1 : C CCREEP SOURCE FANDEUR 14/01/22 21:15:01 7891
2 : SUBROUTINE CCREEP ( wrk52, wrk53, wrk54,
3 : & IFOURB, IB, IGAU, NBPGAU,
4 : & wcreep, iecou, xecou )
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8(A-H,O-Z)
7 : C-----------------------------------------------------------------------
8 : C
9 : C DESCRIPTION FONCTIONNELLE :
10 : C -------------------------
11 : C
12 : C Integration des lois 'VISCO_EXTERNE' par un schema Runge-Kutta 1.2
13 : C
14 : C Entrees/sorties : segments wrk52, wrk53, wrk54 de l'objet DECHE
15 : C Entrees : IFOURB = -3 EN DEFORMATIONS PLANES GENERALISEES
16 : C -2 EN CONTRAINTES PLANES
17 : C -1 EN DEFORMATIONS PLANES
18 : C 0 EN AXISYMETRIE
19 : C 1 EN SERIE DE FOURIER
20 : C 2 EN TRIDIM
21 : C IB = NUMERO DE L'ELEMENT COURANT
22 : C IGAU = NUMERO DU POINT COURANT
23 : C NBPGAU = NBRE DE POINTS DE GAUSS
24 : C Variables de travail : segments wcreep, iecou, xecou
25 : C
26 : C-----------------------------------------------------------------------
27 : -INC CCOPTIO
28 : -INC DECHE
29 : C
30 : C Segment de travail pour les lois 'VISCO_EXTERNE'
31 : C
32 : SEGMENT WCREEP
33 : C Entrees/sorties constantes de la routine CREEP
34 : REAL*8 SERD
35 : CHARACTER*16 CMNAMC
36 : INTEGER LEXIMP, NSTTVC, LAYERC, KSPTC
37 : C Entrees/sorties de la routine CREEP pouvant varier
38 : REAL*8 STV(NSTV), STV1(NSTV), STVP1(NSTV),
39 : & STVP2(NSTV), STV12(NSTV), STVP3(NSTV),
40 : & STVP4(NSTV), STV13(NSTV), STVF(NSTV),
41 : & TMP12, TMP, TMP32,
42 : & DTMP12, DTMP,
43 : & PRD12(NPRD), PRD(NPRD), PRD32(NPRD),
44 : & DPRD12(NPRD), DPRD(NPRD)
45 : INTEGER KSTEPC
46 : C Autres indicateurs et variables de travail
47 : LOGICAL LTMP, LPRD, LSTV
48 : INTEGER IVIEX, NPAREC
49 : REAL*8 dTMPdt, dPRDdt(NPRD)
50 : ENDSEGMENT
51 : C
52 : SEGMENT IECOU
53 : INTEGER NYOG ,NYNU ,NYALFA,NYSMAX,NYN ,NYM ,NYKK ,
54 : & NYALF1,NYBET1,NYR ,NYA ,NYRHO ,NSIGY ,NNKX ,NYKX ,icow16,
55 : & icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
56 : & icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
57 : & icow32,icow33,NSTRS1,MFR1, NBGMAT,NELMAT,MSOUPA,
58 : & icow39,icow40,icow41,icow42,icow43,icow44
59 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
60 : . icow51,icow52,icow53,icow54,icow55,icow56
61 : . icow57,icow58
62 : ENDSEGMENT
63 : C
64 : SEGMENT XECOU
65 : REAL*8 xcow1, xcow2,xcow3,DTT, DT, TREF, xcow7
66 : ENDSEGMENT
67 : C
68 : C Variables locales
69 : C
70 : C Tableaux de travail
71 : C
72 : REAL*8 SIG(8),SIG1(8),SIG12(8),SIG13(8),
73 : & DSPT(8),XX(8),
74 : & EVP1(8),EVP2(8),EVP3(8),EVP4(8),
75 : & EPSV(8),EPSV1(8),EPSV12(8),EPSV13(8),
76 : & CRIGI(12),
77 : & DECRA(5),DESWA(5),TIME12(2),TIME(2),TIME32(2)
78 : C
79 : C Variables scalaires
80 : C
81 : LOGICAL DTLIBR
82 : C
83 : C------------------- Debut du code executable --------------------------
84 : C
85 : C=======================================================================
86 : C 1 - INITIALISATIONS
87 : C Parametres de pilotage des iterations internes
88 : C=======================================================================
89 : C
90 : PRELOC=1.d-8
91 : KERRE = 0
92 : C* -> Les tests de "Restriction" sont maintenant effectues lors de la
93 : C* ceation du modele (cf MODELI).
94 : C*C
95 : C*C Restriction pour l'instant a 'VISCO_EXTERNE' 'GENERAL'
96 : C*C
97 : C* IF ( IVIEX.NE.1 ) THEN
98 : C* KERRE = 958
99 : C* RETURN
100 : C* ENDIF
101 : C*C
102 : C*C Restriction aux elements massifs avec option de calcul 3D
103 : C*C
104 : C* IF ( MFR1.NE.1.OR.IFOURB.NE.2 ) THEN
105 : C* KERRE = 950
106 : C* RETURN
107 : C* ENDIF
108 : C*C
109 : C* Normalement ce test sur DT < 0 devrait etre fait avant (COMVAL ?)
110 : IF (DT.LT.0.0) THEN
111 : KERRE = 414
112 : RETURN
113 : ENDIF
114 : IF (DT.EQ.0.0) DT = 1.e-20
115 : C
116 : DTLIBR =.TRUE.
117 : C
118 : DTLEFT = DT
119 : BORNE = 2.0
120 : RMAX = 1.3
121 : RMIN = 0.7
122 : DIV = 7.0
123 : FAC = 3.0
124 : C
125 : XMAX = XMAT(1)*1.D-3
126 : C
127 : C=======================================================================
128 : C 2 - Prediction elastique de l'increment de contraintes
129 : C=======================================================================
130 : C
130 : CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,N2EL,
131 : & N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST,NBPGAU,MELE,
132 : & NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,
133 : & D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
134 : C
135 : IF (IRTD.NE.1) THEN
136 : KERRE = 714
137 : RETURN
138 : ENDIF
139 : C
140 : C=======================================================================
141 : C 3 - Pre-traitements avant iterations en sous-increments
142 : C=======================================================================
143 : C
144 : C.....Deformations inelastiques au debut du pas :
145 : C Passage gamma -> epsilon pour les termes extradiagonaux
146 : C
147 : C ATTENTION : adaptations si extension a d'autres formulations EF
148 : C ou a d'autres options de calcul
149 : DO 30 I=4,NSTRS1
150 : EPIN0(I)=0.5D0*EPIN0(I)
151 : 30 CONTINUE
152 : C
153 : C.....Calcul des derivees temporelles de la temperature et des
154 : C parametres externes
155 : C Ces parametres de chargement sont supposes varier lineairement
156 : C au cours du pas de temps
157 : C
158 : IF ( LTMP ) THEN
159 : dTMPdt = (TUREF(1)-TURE0(1))/DT
160 : ENDIF
161 : IF ( LPRD ) THEN
162 : DO 35 I=1,NPAREC
163 : dPRDdt(I) = (PAREXF(I)-PAREX0(I))/DT
164 : 35 CONTINUE
165 : ENDIF
166 : C
167 : C=======================================================================
168 : C 4 - INITIALISATIONS AVANT ITERATIONS EN SOUS-INCREMENTS
169 : C=======================================================================
170 : C
171 : ITERO = 0
172 : 6543 CONTINUE
173 : C
174 : ITERO = 1 + ITERO
175 : IF ( ITERO.NE.1) THEN
176 : DTLIBR = .TRUE.
177 : preloc = preloc * 7.D0
178 : IF (ITERO.GT.3) THEN
179 : KERRE = 268
180 : RETURN
181 : ENDIF
182 : ENDIF
183 : C
184 : DTLEFT = DT
185 : TAU = DTLEFT
186 : TAU12 = TAU*0.5D0
187 : C
188 : TIME12(1) = 0.5D0*DT
189 : TIME12(2) = temp0+TIME12(1)
190 : TIME(1)= DT
191 : TIME(2)= tempf
192 : TIME32(1) = 1.5D0*DT
193 : TIME32(2) = temp0+TIME32(1)
194 : C
195 : IF ( LTMP ) THEN
196 : DTMP = TUREF(1)-TURE0(1)
197 : DTMP12 = 0.5D0*DTMP
198 : TMP12 = TURE0(1)+DTMP12
199 : TMP = TUREF(1)
200 : TMP32 = TUREF(1)+DTMP12
201 : ENDIF
202 : IF ( LPRD ) THEN
203 : DO 36 I=1,NPAREC
204 : DPRD(I) = PAREXF(I)-PAREX0(I)
205 : DPRD12(I) = 0.5D0*DPRD(I)
206 : PRD12(I) = PAREX0(I)+DPRD12(I)
207 : PRD(I) = PAREXF(I)
208 : PRD32(I) = PAREXF(I)+DPRD12(I)
209 : 36 CONTINUE
210 : ENDIF
211 : C
212 : ASIG = SQRT(PROCON(SIG0,SIG0,NSTRS1))
213 : ERRABS = preloc*ASIG
214 : IF (XMAX.GT.ASIG) ERRABS = preloc*XMAX
215 : C
216 : DO 40 I=1,NSTRS1
217 : SIG(I) = SIG0(I)
218 : EPSV(I) = EPIN0(I)
219 : DSPT(I) = DSIGT(I)/DT
220 : 40 CONTINUE
221 : C
222 : EC0 = VAR0(1)
223 : ESW0 = VAR0(2)
224 : P = VAR0(3)
225 : QTLD = VAR0(4)
226 : IF ( LSTV ) THEN
227 : DO 50 I=1,NSTTVC
228 : STV(I)=VAR0(4+I)
229 : 50 CONTINUE
230 : ENDIF
231 : C
232 : C=======================================================================
233 : C 5 - ITERATIONS EN SSINCREMENTS /FIN SI DTLEFT = 0
234 : C=======================================================================
235 : C
236 : NSSINC = 0
237 : NITERA = 0
238 : 60 CONTINUE
239 : C
240 : NSSINC = NSSINC + 1
241 : IF (NSSINC.GT.MSOUPA) THEN
242 : DTLIBR=.FALSE.
243 : GOTO 6543
244 : ENDIF
245 : C
246 : C Evaluation initiale des vitesses de deformation inelastique
247 : C sur la base du dernier sous-increment converge
248 : LEND = 0
248 : CALL CREEP(DECRA,DESWA,STV,serd,EC0,ESW0,P,QTLD,
249 : & TMP,DTMP,PRD,DPRD,TIME,TAU,
250 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
251 : & layerc,ksptc,KSTEPC,NSSINC)
252 : IF (KSTEPC.NE.1) RETURN
252 : CALL INCREE(TAU,SIG,P,QTLD,DECRA,DESWA,
253 : & EVP1,EC0P1,ESW0P1,
254 : & NSTRS1,IVIEX,MFR1,IFOURB)
255 : C
256 : C ---------------------------------------------------------------
257 : C DEBUT DES ITERATIONS SUR TAU OPTIMAL /FIN SI RA PETIT
258 : C
259 : NITERA = 0
260 : 70 CONTINUE
261 : C
262 : NITERA = NITERA + 1
263 : C
264 : C ____________________________________________________________
265 : C Premiere evaluation de l'etat a t+TAU12
265 : CALL AVCREE(TAU12,SIG,EPSV,EC0,ESW0,
266 : & SIG1,EPSV1,EC01,ESW01,P1,QTLD1,
267 : & DSPT,EVP1,EC0P1,ESW0P1,XMAT,
268 : & NSTRS1,IVIEX,MFR1,IFOURB)
269 : C MAJ des variables internes supplementaires le cas echeant
270 : IF ( LSTV ) THEN
271 : LEND = 1
271 : CALL CREEP(DECRA,DESWA,STV1,serd,EC0,ESW0,P1,QTLD1,
272 : & TMP12,DTMP12,PRD12,DPRD12,TIME12,TAU12,
273 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
274 : & layerc,ksptc,KSTEPC,NSSINC)
275 : IF (KSTEPC.NE.1) RETURN
276 : ENDIF
277 : C
278 : C ____________________________________________________________
279 : C Reevaluation des vitesses de deformation inelastique sur la
280 : C base de l'etat a t+TAU12 calcule precedemment, puis moyenne
281 : LEND = 0
281 : CALL CREEP(DECRA,DESWA,STV1,serd,EC01,ESW01,P1,QTLD1,
282 : & TMP32,DTMP,PRD32,DPRD,TIME32,TAU,
283 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
284 : & layerc,ksptc,KSTEPC,NSSINC)
285 : IF (KSTEPC.NE.1) RETURN
285 : CALL INCREE(TAU,SIG1,P1,QTLD1,DECRA,DESWA,
286 : & EVP2,EC0P2,ESW0P2,
287 : & NSTRS1,IVIEX,MFR1,IFOURB)
288 : C
289 : DO 701 I=1,NSTRS1
290 : EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
291 : 701 CONTINUE
292 : EC0P2 = 0.5D0*(EC0P1+EC0P2)
293 : ESW0P2 = 0.5D0*(ESW0P1+ESW0P2)
294 : C
295 : C ____________________________________________________________
296 : C Reevaluation de l'etat a t+TAU12
296 : CALL AVCREE(TAU12,SIG,EPSV,EC0,ESW0,
297 : & SIG12,EPSV12,EC012,ESW012,P12,QTLD12,
298 : & DSPT,EVP2,EC0P2,ESW0P2,XMAT,
299 : & NSTRS1,IVIEX,MFR1,IFOURB)
300 : C MAJ des variables internes supplementaires le cas echeant
301 : IF ( LSTV ) THEN
302 : LEND = 1
302 : CALL CREEP(DECRA,DESWA,STV12,serd,EC0,ESW0,P12,QTLD12,
303 : & TMP12,DTMP12,PRD12,DPRD12,TIME12,TAU12,
304 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
305 : & layerc,ksptc,KSTEPC,NSSINC)
306 : IF (KSTEPC.NE.1) RETURN
307 : ENDIF
308 : C
309 : C ____________________________________________________________
310 : C Reevaluation des vitesses de deformation inelastique sur la
311 : C base de l'etat a t+TAU12 calcule precedemment
312 : LEND = 0
312 : CALL CREEP(DECRA,DESWA,STV12,serd,EC012,ESW012,P12,QTLD12,
313 : & TMP,DTMP12,PRD,DPRD12,TIME,TAU12,
314 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
315 : & layerc,ksptc,KSTEPC,NSSINC)
316 : IF (KSTEPC.NE.1) RETURN
316 : CALL INCREE(TAU12,SIG12,P12,QTLD12,DECRA,DESWA,
317 : & EVP3,EC0P3,ESW0P3,
318 : & NSTRS1,IVIEX,MFR1,IFOURB)
319 : C
320 : C ____________________________________________________________
321 : C Premiere evaluation de l'etat a t+TAU
321 : CALL AVCREE(TAU12,SIG12,EPSV12,EC012,ESW012,
322 : & SIG13,EPSV13,EC013,ESW013,P13,QTLD13,
323 : & DSPT,EVP3,EC0P3,ESW0P3,XMAT,
324 : & NSTRS1,IVIEX,MFR1,IFOURB)
325 : C MAJ des variables internes supplementaires le cas echeant
326 : IF ( LSTV ) THEN
327 : LEND = 1
327 : CALL CREEP(DECRA,DESWA,STV13,serd,EC012,ESW012,P13,QTLD13,
328 : & TMP,DTMP12,PRD,DPRD12,TIME,TAU12,
329 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
330 : & layerc,ksptc,KSTEPC,NSSINC)
331 : IF (KSTEPC.NE.1) RETURN
332 : ENDIF
333 : C
334 : C ____________________________________________________________
335 : C Reevaluation des vitesses de deformation inelastique sur la
336 : C base de l'etat a t+TAU calcule precedemment, puis moyenne
337 : LEND = 0
337 : CALL CREEP(DECRA,DESWA,STV13,serd,EC013,ESW013,P13,QTLD13,
338 : & TMP32,DTMP12,PRD32,DPRD12,TIME32,TAU12,
339 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
340 : & layerc,ksptc,KSTEPC,NSSINC)
341 : IF (KSTEPC.NE.1) RETURN
341 : CALL INCREE(TAU12,SIG13,P13,QTLD13,DECRA,DESWA,
342 : & EVP4,EC0P4,ESW0P4,
343 : & NSTRS1,IVIEX,MFR1,IFOURB)
344 : C
345 : DO 711 I=1,NSTRS1
346 : EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
347 : 711 CONTINUE
348 : EC0P4 = 0.5D0*(EC0P3+EC0P4)
349 : ESW0P4 = 0.5D0*(ESW0P3+ESW0P4)
350 : C
351 : C ____________________________________________________________
352 : C Reevaluation de l'etat a t+TAU
352 : CALL AVCREE(TAU12,SIG12,EPSV12,EC012,ESW012,
353 : & SIGF,EPINF,EC0F,ESW0F,PF,QTLDF,
354 : & DSPT,EVP4,EC0P4,ESW0P4,XMAT,
355 : & NSTRS1,IVIEX,MFR1,IFOURB)
356 : C MAJ des variables internes supplementaires le cas echeant
357 : IF ( LSTV ) THEN
358 : LEND = 1
358 : CALL CREEP(DECRA,DESWA,STVF,serd,EC012,ESW012,PF,QTLDF,
359 : & TMP,DTMP12,PRD,DPRD12,TIME,TAU12,
360 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
361 : & layerc,ksptc,KSTEPC,NSSINC)
362 : IF (KSTEPC.NE.1) RETURN
363 : ENDIF
364 : C
365 : C ____________________________________________________________
366 : C Reevaluation des vitesses de deformation inelastique sur la
367 : C base de l'etat a t+TAU calcule precedemment, puis moyenne
368 : LEND = 0
368 : CALL CREEP(DECRA,DESWA,STVF,serd,EC0F,ESW0F,PF,QTLDF,
369 : & TMP32,DTMP12,PRD32,DPRD12,TIME32,TAU12,
370 : & cmnamc,leximp,LEND,coorga,nsttvc,ib,igau,
371 : & layerc,ksptc,KSTEPC,NSSINC)
372 : IF (KSTEPC.NE.1) RETURN
372 : CALL INCREE(TAU12,SIGF,PF,QTLDF,DECRA,DESWA,
373 : & EVP4,EC0P4,ESW0P4,
374 : & NSTRS1,IVIEX,MFR1,IFOURB)
375 : C
376 : DO 721 I=1,NSTRS1
377 : EVP2(I) = (EVP1(I)+EVP4(I))/6.D0+EVP3(I)*2.D0/3.D0
378 : 721 CONTINUE
379 : EC0P2 = (EC0P1+EC0P4)/6.D0+EC0P3*2.D0/3.D0
380 : ESW0P2 = (ESW0P1+ESW0P4)/6.D0+ESW0P3*2.D0/3.D0
381 : C
382 : C ____________________________________________________________
383 : C Derniere evaluation de l'etat a t+TAU pour test convergence
383 : CALL AVCREE(TAU,SIG,EPSV,EC0,ESW0,
384 : & SIG1,EPSV1,EC01,ESW01,P1,QTLD1,
385 : & DSPT,EVP2,EC0P2,ESW0P2,XMAT,
386 : & NSTRS1,IVIEX,MFR1,IFOURB)
387 : C Pas de MAJ des variables internes supplementaires (inutile)
388 : C
389 : C ____________________________________________________________
390 : C CALCUL DU RAPPORT : ERREUR CALCULEE / ERREUR ADMISE
391 : C
392 : DO 730 I=1,NSTRS1
393 : XX(I) = SIGF(I)-SIG1(I)
394 : 730 CONTINUE
395 : RA = SQRT(PROCON(XX,XX,NSTRS1))/(ERRABS)
396 : SQRA = SQRT(RA)
397 : C
398 : C ____________________________________________________________
399 : C TEST DE FIN D'ITERATIONS / MISE A JOUR DE TAU
400 : C DIV =7 BORNE = 2
401 : C SI SQRA>7 TAU = TAU/7 ET NOUVEL ESSAI
402 : C SI 2RMAX on vise RA = 1
504 : C Fin des boucles en ss increments si tau = DTLEFT
505 : C
506 : IF ( TAU.LT.DTLEFT ) THEN
507 : DTLEFT = DTLEFT - TAU
508 : IF ( FAC*SQRA.LT.1.D0) THEN
509 : TAU=TAU*FAC
510 : ELSE IF ( (SQRA.LT.RMIN).OR.(SQRA.GT.RMAX) ) THEN
511 : TAU=TAU/SQRA
512 : ENDIF
513 : IF (TAU.GT.DTLEFT) TAU = DTLEFT
514 : TAU12 = TAU*0.5D0
515 : TIME12(1) = TIME(1) + TAU12
516 : TIME12(2) = TIME(2) + TAU12
517 : TIME32(1) = TIME(1) + 1.5D0*TAU
518 : TIME32(2) = TIME(2) + 1.5D0*TAU
519 : TIME(1) = TIME(1) + TAU
520 : TIME(2) = TIME(2) + TAU
521 : IF ( LTMP ) THEN
522 : DTMP = dTMPdt * TAU
523 : DTMP12 = dTMPdt * TAU12
524 : TMP12 = TMP + DTMP12
525 : TMP32 = TMP + 1.5D0*DTMP
526 : TMP = TMP + DTMP
527 : ENDIF
528 : IF ( LPRD ) THEN
529 : DO 100 I=1,NPAREC
530 : DPRD(I) = dPRDdt(I) * TAU
531 : DPRD12(I) = dPRDdt(I) * TAU12
532 : PRD12(I) = PRD(I) + DPRD12(I)
533 : PRD32(I) = PRD(I) + 1.5D0*DPRD(I)
534 : PRD(I) = PRD(I) + DPRD(I)
535 : 100 CONTINUE
536 : ENDIF
537 : GOTO 60
538 : ENDIF
539 : C
540 : IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
541 : WRITE ( IOIMP,* ) ' PROBLEME TAU > DTLEFT '
542 : KERRE = 223
543 : RETURN
544 : ENDIF
545 : C
546 : C=======================================================================
547 : C 6 - Post-traitements pour stockage des resultats
548 : C=======================================================================
549 : C
550 : C.....Deformations inelastiques debut/fin de pas :
551 : C Passage epsilon -> gamma pour les termes extradiagonaux
552 : C ATTENTION : adaptations si extension a d'autres formulations EF
553 : C ou a d'autres options de calcul
554 : DO 110 I=1,NSTRS1
555 : A=1.D0
556 : IF (I.GT.3) A=2.D0
557 : EPIN0(I)=EPIN0(I)*A
558 : EPINF(I)=EPINF(I)*A
559 : 110 CONTINUE
560 : C
561 : C.....Variables internes a la fin du pas
562 : C
563 : VARF(1) = EC0F
564 : VARF(2) = ESW0F
565 : VARF(3) = PF
566 : VARF(4) = QTLDF
567 : IF ( LSTV ) THEN
568 : DO 120 I=1,NSTTVC
569 : VARF(4+I)=STVF(I)
570 : 120 CONTINUE
571 : ENDIF
572 : C
573 : RETURN
574 : END
575 :
576 :
577 :
578 :
579 :
580 :
581 :
582 :