Download ccerac.eso
Back to the list
1 : C CCERAC SOURCE KICH 09/09/11 21:15:12 6505
2 : SUBROUTINE CCERAC(wrk52,wrk53,wrk54,NVARI,
3 : 1 NSSINC,INV,IFOURB,IB,IGAU,NBPGAU,iecou,xecou)
4 : C CERACA SOURCE AM 00/12/13 21:15:32 4045
5 : c SUBROUTINE CERACA(WRK0,WR00,WRK1,WRK5,WRK7,WRK8,WRK9,WTRAV,
6 : c 1 INPLAS,MFR,DT,NSTRS,NVARI,NCOMAT,PRECIS,MSOUPA,JECHER,DTT,
7 : c 2 NSSINC,INV,KERRE,ICARA,IFOURB,NYOG,NYNU,NYALFA,NYSMAX,
8 : c 3 NYN,NYM,NYKK,NYALF1,NYBET1,NYR,NYA,NYKX,NNKX,NYRHO,NSIGY,T0,TF,
9 : c 5 TREF,TLIFE,ITHHER,NCOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST,
10 : c 7 NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,CRIGI,
11 : c 8 KERREU1)
12 : C
13 : C---------------------------------------------------------------------
14 : C Objet: Calculer au cours d'un pas de temps DT, l'evolution des
15 : C variables internes a l'aide d'un schema Runge-Kutta 1.2
16 : C ---------------------------------------------------------------------
17 : C MFR1 <- MFR, XCARB <- XCAR, NSTRS1 <- NSTRS,
18 : C
19 : C---------------------------------------------------------------------
20 : C Entree: INPLAS type de materiau
21 : C MFR1 indice de la formulation mecanique(seulement massif ou coque
22 : C pour les materiaux endommageables)
23 : C DEPST(NSTRS1) increment des deformations totales
24 : C SIG0(NSTRS1) contraintes initiales
25 : C EPIN0(NSTRS1) deformations viscoplastiques initiales
26 : C VAR0(NVARI) variables internes initiales
27 : C NVARI nombre de variables internes
28 : C SIGY(NSIGY) courbe de la limite elastique en fonction de T°C
29 : C XMAT(NCOMAT) materiau
30 : C XCARB(ICARA) caracteristiques geometriques
31 : C YSMAX(NYSMAX) intervient ds. le test de convergence des iter.
32 : C PRECIS precision relative sur SIGMA
33 : C MSOUPA nombre maximal de sous pas autorises
34 : C JECHER = 0 avancer
35 : C = 1 rechercher sortie avec DTT
36 : C IFOURB = -2 EN CONTR.PLANES
37 : C -1 EN DEFORM. PLANES
38 : C 0 EN AXISYMETRIE
39 : C 1 EN SERIE DE FOURIER
40 : C 2 EN TRIDIM
41 : C CMATE = NOM DU MATERIAU
42 : C VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
43 : C VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
44 : C N2EL = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
45 : C N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
46 : C IB = NUMERO DE L ELEMENT COURANT
47 : C IGAU = NUMERO DU POINT COURANT
48 : C EPAIST= EPAISSEUR
49 : C NBPGAU= NBRE DE POINTS DE GAUSS
50 : C MELE = NUMERO DE L ELEMENT FINI
51 : C NPINT = NBRE DE POINTS D INTEGRATION
52 : C SECT = SECTION
53 : C LHOOK = TAILLE DE LA MATRICE DE HOOKE
54 : C DT pas de temps
55 : C ITHHER = 0 pas de chargement thermique et materiau constant
56 : C = 1 chargement thermique et materiau constant
57 : C = 2 chargement thermique et materiau(T)
58 : C-----------------------------------------------------------------------
59 : C
60 : C-----------------------------------------------------------------------
61 : C Sortie: SIGF(NSTRS) contraintes finales
62 : C EPINF(NSTRS) deformations viscoplastiques finales
63 : C VARF(NVARI) variables internes finales
64 : C DTT sous-increment de temps optimal (si JECHER=1)
65 : C NSSINC nombre de sous-increments si JECHER=0
66 : C INV = 1 si inversion
67 : C 0 sinon
68 : C KERRE = 0 si tout OK
69 : C <> 0 si entrees incoherentes
70 : C-----------------------------------------------------------------------
71 : C
72 : IMPLICIT INTEGER(I-N)
73 : IMPLICIT REAL*8(A-H,O-Z)
74 : -INC CCOPTIO
75 : -INC DECHE
76 :
77 : SEGMENT IECOU
78 : * COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
79 : INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
80 : C INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM, NYKK,
81 : 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
82 : C . NYALF1,NYBET1,NYR, NYA, NYRHO,NSIGY, NNKX, NYKX, IND,
83 : 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
84 : C . NSOM, NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT, NCXMAT,
85 : 3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
86 : C . LTRAC, MFR, IELE, NHRM, NBNN, NBELEM,ICARA,
87 : 4 icow32,icow33,NSTRS1,MFR1,NBGMAT,NELMAT,MSOUPA,
88 : C . LW2, NDEF, NSTRSS,MFR1, NBGMAT,NELMAT,MSOUPA,
89 : 5 icow39,icow40,icow41,icow42,icow43,icow44
90 : C . NUMAT1,LENDO, NBBB, NNVARI,KERR1, MELEME
91 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
92 : . icow51,icow52,icow53,icow54,icow55,icow56
93 : . icow57,icow58
94 : ENDSEGMENT
95 :
96 : SEGMENT XECOU
97 : * COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
98 : REAL*8 xcow1, xcow2,xcow3,DTT,DT,xcow6, xcow7
99 : C REAL*8 DTOPTI,TSOM, TCAR, DTT, DT, TREFA,TEMP00
100 : ENDSEGMENT
101 :
102 :
103 :
104 : *
105 : *
106 : DIMENSION SIG(8),EPSV(8),VAR(100)
107 : DIMENSION DSPT(8),XX(8),SIG1(8),EPSV1(8)
108 : DIMENSION VAR1(100),EVP1(8),VARP1(100)
109 : DIMENSION EVP2(8),VARP2(100)
110 : DIMENSION CRIGI(12)
111 : C**** debut ajout Eloi
112 : logical dtlibr
113 : C**** fin ajout Eloi
114 : NCOMAT = NMATT
115 : C**** debut ajout Eloi
116 : dtlibr=.TRUE.
117 : C**** fin ajout Eloi
118 : C ==================================================
119 : C RECHERCHE DU NUMERO DE LA VARIABLE INTERNE EGALE A 1
120 : C SI ON A ENDOMMAGEMENT GENERALISE
121 : C
122 : IBID2 = NVARI-1
123 : C*************************************
124 : C On cherche le numero de la propriété du matériau correspondant à ENDG
125 : C
126 : IF (MFR1.EQ.1.AND.IFOMOD.EQ.2) THEN
127 : IBID1 = 20
128 : ELSE
129 : IBID1 = 15
130 : ENDIF
131 : C Test pour savoir si on a dépassé la limite de déformation de fluage
132 : C
133 : IF (VAR0(IBID2).EQ.1.D0) THEN
134 : DO 10 IJ=1,NVARI
135 : VARF(IJ) = VAR0(IJ)
136 : 10 CONTINUE
137 : DO 20 K=1,NSTRS1
138 : EPINF(K) = EPIN0(K)
139 : SIGF(K) = 0.D0
140 : 20 CONTINUE
141 : GOTO 1000
142 : ENDIF
143 : DO 30 I=1,NSTRS1
144 : TOTO = DABS(EPIN0(I))
145 : IF (TOTO.GE.XMAT(IBID1)) THEN
146 : DO 40 K=1,NSTRS1
147 : EPINF(K) = EPIN0(K)
148 : SIGF(K) = 0.D0
149 : 40 CONTINUE
150 : DO 50 IJ=1,(NVARI-1)
151 : VARF(IJ) = VAR0(IJ)
152 : 50 CONTINUE
153 : VARF(IBID2) = 1.D0
154 : GOTO 1000
155 : ENDIF
156 : 30 CONTINUE
157 : C*************************************
158 :
159 :
160 : C =========================================================
161 : KERRE = 0
162 : IF (MFR1.NE.1.AND.MFR1.NE.3.OR.IFOURB.EQ.1) THEN
163 : KERRE = 99
164 : RETURN
165 : ENDIF
166 :
167 : IF (MFR1.EQ.3) THEN
168 : THICK = XCARB(1)
169 : ALFA = XCARB(2)
170 : ENDIF
171 : BORNE = 2.0
172 : RMAX = 1.3
173 : RMIN = 0.7
174 : DIV = 7.0
175 : FAC = 3.0
176 : C
177 : C CALCUL DES INCREMENTS DE CONTRAINTES ELASTIQUES
178 : C
178 : CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,N2EL,N2PTEL,
179 : 1 MFR1,IFOURB,IB,IGAU,EPAIST,NBPGAU,MELE,NPINT,NBGMAT,
180 : 2 NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,
181 : 3 CRIGI,DSIGT,IRTD)
182 : C PRINT *,'DSIGT =',DSIGT
183 : *
184 : IF (IRTD.NE.1) THEN
185 : KERRE=69
186 : GOTO 1000
187 : ENDIF
188 : *
189 : IF (MFR1.EQ.3) THEN
190 : DO 60 I=1,NSTRS1/2
191 : SIG0( I) = SIG0( I)/THICK
192 : SIG0(NSTRS1/2+I) = SIG0(NSTRS1/2+I)*6.0D0/THICK/THICK
193 : DSIGT( I) = DSIGT( I)/THICK
194 : DSIGT(NSTRS1/2+I)= DSIGT(NSTRS1/2+I)*6.0D0/THICK/THICK
195 : 60 CONTINUE
196 : IF (IFOURB.EQ.-2) THEN
197 : SIG0 (2) = 0.D0
198 : SIG0 (4) = 0.D0
199 : DSIGT(2) = 0.D0
200 : DSIGT(4) = 0.D0
201 : ENDIF
202 : ENDIF
203 : C
204 : C------------------------------------------
205 : C CONTROLE DE LA COHERENCE DES ENTREES
206 : C Supprimée dans ce cas
207 : C------------------------------------------
208 : IF (DT.LE.0.D0) KERRE = 414
209 : MOTERR(1:8) = ' CONST '
210 : C IF (INCOMAT.NE. 15)KERRE = 146
211 : C IF (IFOURB.EQ.-2.AND.NCOMAT.NE. 9)KERRE = 146
212 : IF (IFOURB.NE.2.AND.IFOURB.NE.-2.AND.IFOURB.NE.0
213 : & .AND.IFOURB.NE.-1) THEN
214 : KERRE = 194
215 : MOTERR(1:8) = 'FLUAGE'
216 : ENDIF
217 : IF (MFR1.EQ.1.AND.IFOMOD.EQ.2) THEN
218 : IBID = 15
219 : ELSE
220 : IBID = 10
221 : ENDIF
222 : XMAX=XMAT(IBID)
223 : *
224 : * TEST SUR XMAX MILL 8/3/91
225 : *
226 : IF (XMAX.EQ.0.D0) THEN
227 : XMAX=XMAT(1)*1.D-3
228 : ENDIF
229 : C
230 : C-----------------------------
231 : IF (KERRE.NE.0) THEN
232 : GOTO 1100
233 : ENDIF
234 : C
235 : C===========================================================
236 : C A PARTIR DE MAINTENANT, LES DEFORMATIONS
237 : C DE CISAILLEMENT NE SONT PLUS
238 : C DEFINIES PAR DES GAMA.
239 : C ON DIVISE DONC LES TERMES DE CISAILLEMENT PAR 2.
240 : C SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CERACARO:
241 : C MFR1=1 (MASSIF)
242 : C MFR1=3 (COQUES MINCES)
243 : C
244 : C
245 : C Cas de la formulation massive
246 : C Les termes de cisaillement apparaissent
247 : C au delà de la troisieme composante
248 : C
249 : IF (MFR1.EQ.1) THEN
250 : DO 70 I=1,NSTRS1
251 : A=1.D0
252 : IF (I.GT.3) A=2.D0
253 : EPIN0(I)=EPIN0(I)/A
254 : 70 CONTINUE
255 : C
256 : C Cas des coques minces
257 : C Les termes de cisaillement apparaissent
258 : C pour la troisieme et la sixieme composante
259 : C uniquement dans les cas de calculs
260 : C tridimensionnels ou d'analyse de Fourier
261 : C
262 : ELSE IF (MFR1.EQ.3) THEN
263 : IF (IFOURB.EQ.2) THEN
264 : DO 80 I=1,NSTRS1
265 : A=1.D0
266 : IF (I.EQ.3) A=2.D0
267 : IF (I.EQ.6) A=2.D0
268 : EPIN0(I)=EPIN0(I)/A
269 : 80 CONTINUE
270 : ENDIF
271 : ENDIF
272 : C
273 : C===========================================================
274 : C
275 : C ----------------
276 : C INITIALISATION
277 : C ----------------
278 : C**** debut ajout Eloi
279 : ITERO = 0
280 : 1200 CONTINUE
281 :
282 : itero = 1 + itero
283 : if (itero.ne.1) THEN
284 : c write(6,*) 'itero ib igau', itero,ib,igau
285 : dtlibr = .true.
286 : precis = precis * 7.
287 : c write(6,*) ' precision modifiée ', precis
288 : if (itero.gt.3) then
289 : **** kerre = 460
290 : kerre = 268
291 : return
292 : endif
293 : endif
294 : C**** fin ajout Eloi
295 : DTLEFT = DT
296 : TAU = DTLEFT
297 : ASIG = SQRT(PROCON(SIG0,SIG0,NSTRS1))
298 : ERRABS = PRECIS*ASIG
299 : IF (XMAX.GT.ASIG) ERRABS = PRECIS*XMAX
300 : DO 90 I=1,NSTRS1
301 : SIG(I) = SIG0(I)
302 : EPSV(I) = EPIN0(I)
303 : DSPT(I) = DSIGT(I)/DT
304 : 90 CONTINUE
305 : C
306 : DO 100 I=1,NVARI
307 : VAR(I)=VAR0(I)
308 : 100 CONTINUE
309 : C
310 : C ---------------------------------------------------------------------
311 : INV = 0
312 : NSSINC = 0
313 : C ---------------------------------------------------------------------
314 : C DEBUT DES ITERATIONS EN SSINCREMENTS /FIN SI DTLEFT = 0
315 : C ---------------------------------------------------------------------
316 : 1400 CONTINUE
317 : NSSINC = NSSINC + 1
318 : IF (NSSINC.GT.MSOUPA) THEN
319 : C**** debut ajout Eloi
320 : DTLIBR=.FALSE.
321 : GOTO 1200
322 : C**** fin ajout Eloi
323 : C*** KERRE = 460
324 : C*** GOTO 1100
325 : ENDIF
326 : C
327 : C---------------------------------------------------------------------
328 : C START OF CALCULATIONS
329 : C_____________________________________________________________________
330 :
331 : IF(MFR1.NE.3) THEN
331 : CALL INCR11(TAU,SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,MFR1)
332 : ELSE IF(MFR1.EQ.3) THEN
332 : CALL INCR33(TAU,SIG,EPSV,VAR,XMAT,EVP1,VARP1,ALFA,NSTRS1,NVARI)
333 : ENDIF
334 : C
335 : NITERA = 0
336 :
337 : C --------------------------------------------------------------------
338 : C DEBUT DES ITERATIONS SUR TAU OPTIMAL /FIN SI RA PETIT
339 : C --------------------------------------------------------------------
340 :
341 : 1300 CONTINUE
342 : NITERA = NITERA + 1
343 : IF (MFR1.NE.3) THEN
344 :
344 : CALL ADVA11(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
345 : & XMAT,NSTRS1,NVARI,MFR1)
345 : CALL INCR11(TAU,SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,NVARI,
346 : & MFR1)
347 : DO 110 I=1,NSTRS1
348 : EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
349 : 110 CONTINUE
350 : CCC Eloi : on peut effectuer la boucle suivante sur
351 : CCC les variables internes propres a Ottosen
352 : CCC qui n'ont pas ete modifiees
353 : DO 120 I=1,NVARI
354 : VARP2(I)= 0.5*(VARP1(I) + VARP2(I))
355 : 120 CONTINUE
355 : CALL ADVA11(TAU,SIG,EPSV,VAR,SIGF,EPINF,VARF,DSPT,EVP2,VARP2,
356 : & XMAT,NSTRS1,NVARI,MFR1)
357 : ELSE
358 : C----------------------------------------------------------------------
359 : C CALCULATIONS FOR GENERALISED STRESS/STRAIN FORMULATIONS
360 : C----------------------------------------------------------------------
360 : CALL ADVA22(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
361 : & XMAT,NSTRS1,NVARI)
361 : CALL INCR33(TAU,SIG1,EPSV1,VAR1,XMAT,EVP2,VARP2,ALFA,NSTRS1,
362 : & NVARI)
363 : DO 140 I=1,NSTRS1,1
364 : EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
365 : 140 CONTINUE
366 : CCC Eloi : on peut effectuer la boucle suivante sur
367 : CCC les variables internes propres a Ottosen
368 : CCC qui n'ont pas ete modifiees
369 : DO 150 I=1,NVARI
370 : VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
371 : 150 CONTINUE
371 : CALL ADVA22(TAU,SIG,EPSV,VAR,SIGF,EPINF,VARF,DSPT,EVP2,VARP2,
372 : & XMAT,NSTRS1,NVARI)
373 : ENDIF
374 :
375 : C ---------------------------------------------------------------------
376 : C CALCUL DU RAPPORT : ERREUR CALCULEE / ERREUR ADMISE
377 : C ---------------------------------------------------------------------
378 : DO 170 I=1,NSTRS1
379 : XX(I) = SIGF(I)-SIG1(I)
380 : 170 CONTINUE
381 : RA = SQRT(PROCON(XX,XX,NSTRS1))/(ERRABS)
382 : SQRA = SQRT(RA)
383 :
384 : C ---------------------------------------------------------------------
385 : C TEST DE FIN D'ITERATIONS / MISE A JOUR DE TAU /OPTION JECHER
386 : C DIV =7 BORNE = 2
387 : C SI SQRA>7 TAU = TAU/7 ET NOUVEL ESSAI
388 : C SI 2RMAX on vise RA = 1
467 : C fin des boucles en ss increments si tau = dtleft
468 : C --------------------------------------------------------------------
469 : C
470 : IF (TAU.LT.DTLEFT) THEN
471 : DTLEFT = DTLEFT - TAU
472 : IF (FAC*SQRA.LT.1.D0) THEN
473 : TAU=TAU*FAC
474 : ELSEIF ( (SQRA.LT.RMIN).OR.(SQRA.GT.RMAX) ) THEN
475 : TAU=TAU/SQRA
476 : ENDIF
477 : IF (TAU.GT.DTLEFT) TAU = DTLEFT
478 : GOTO 1400
479 : ENDIF
480 : C
481 : IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
482 : WRITE ( IOIMP,* ) ' PROBLEME TAU > DTLEFT '
483 : KERRE = 223
484 : ENDIF
485 : C-----------------------------------------------------------------------
486 : 1100 CONTINUE
487 : IF (MFR1.EQ.3) THEN
488 : DO 220 I=1,NSTRS1/2
489 : SIGF( I) =SIGF( I)*THICK
490 : SIGF(NSTRS1/2+I) =SIGF(NSTRS1/2+ I)*THICK*THICK/6.0
491 : * DSIGT( I)=DSIGT( I)*THICK
492 : * DSIGT(NSTRS1/2+I)=DSIGT(NSTRS1/2+I)*THICK*THICK/6.0
493 : 220 CONTINUE
494 : ENDIF
495 :
496 : C
497 : C===========================================================
498 : C RETOUR A LA DEFINITION NORMALE DES DEFORMATIONS
499 : C A SAVOIR: LES DEFORMATIONS DE CISAILLEMENT SONT
500 : C DEFINIES PAR DES GAMA.
501 : C SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR Ceracaro:
502 : C MFR1=1 (MASSIF)
503 : C MFR1=3 (COQUES MINCES)
504 : C
505 : C Cas de la formulation massive
506 : C Les termes de cisaillement apparaissent
507 : C au delà de la troisieme composante
508 : C
509 : IF (MFR1.EQ.1) THEN
510 : DO 230 I=1,NSTRS1
511 : A=1.D0
512 : IF (I.GT.3) A=2.D0
513 : EPIN0(I)=EPIN0(I)*A
514 : EPINF(I)=EPINF(I)*A
515 : 230 CONTINUE
516 : C
517 : C Cas des coques minces
518 : C Les termes de cisaillement apparaissent
519 : C pour la troisieme et la sixieme composante
520 : C uniquement dans les cas de calculs
521 : C tridimensionnels ou d'analyse de Fourier
522 : C
523 : ELSE IF (MFR1.EQ.3) THEN
524 : IF (IFOURB.EQ.2) THEN
525 : DO 240 I=1,NSTRS1
526 : A=1.D0
527 : IF (I.EQ.3) A=2.D0
528 : IF (I.EQ.6) A=2.D0
529 : EPIN0(I)=EPIN0(I)*A
530 : EPINF(I)=EPINF(I)*A
531 : 240 CONTINUE
532 : ENDIF
533 : ENDIF
534 : C
535 : C===========================================================
536 : C
537 : 1000 RETURN
538 : END
539 :
540 :
541 :
542 :
543 :
544 :
545 :
546 :
547 :
548 :
549 :
550 :
551 :
552 :