Download uo2dc2.eso
Back to the list
1 : C UO2DC2 SOURCE STRU 07/05/31 21:15:37 5744
2 : SUBROUTINE UO2DC2(IFLAG,T0,TPOINT,FI0,FPOINT,PRECIS,PRECIZ,
3 : 1 XMAT,NCOMAT,NSIMP,AAD,GS,LEBIL,XINVL,PENTE,
4 : 2 EPSPT,SIG0,EPSV0,VAR0,W0,WMAX0,DX0,
5 : 3 NGAT,NC1,NCA,NDIM,NN,TAU,TAUNEX,SIGF,EPSVF,
6 : 4 VARF,WF,DXF,WMAXF,TF,FIF,KERRE)
7 : C----------------------------------------------------------------------
8 : C ECOULEMENT MODELE UO2 (OTTOSEN ET GATT_MONERIE) SIMPLIFIE
9 : C SCHEMA EULER-CAUCHY AVEC DECOUPAGE DU SOUS PAS D INTEGRATION EN 2
10 : C CONTROLE AVEC METHODE DE SIMPSON
11 : C RECHERCHE D UN ETAT CONVERGE POUR TAU <= DT SI IFLAG=1
12 : C----------------------------------------------------------------------
13 : C
14 : C ENTREES
15 : C -------
16 : C IFLAG = 1 RECHERCHE D UN ETAT CONVERGE POUR TAU <= DT
17 : C 0 PAS DE RECHERCHE D UN ETAT CONVERGE (1 ITERATION)
18 : C T0 = TEMPERATURE AU DEBUT DU SOUS PAS D INTEGRATION
19 : C TPOINT = VITESSE DE TEMPERATURE SUR LE PAS D INTEGRATION
20 : C FI0 = DENSITE DE FISSION AU DEBUT DU SOUS PAS D INTEGRATION
21 : C FPOINT = VITESSE DE DENSITE DE FISSION SUR LE PAS D INTEGRATION
22 : C PRECIS = PRECISION POUR LA CONVERGENCE DU SCHEMA NUMERIQUE
23 : C PRECIZ = PRECISION POUR L INVERSION DE LA MATRICE
24 : C XMAT(NCOMAT) = CARACTERISTIQUES THERMOMECANIQUES DU MATERIAU
25 : C NSIMP = POINTE SUR LA CARACTERISTIQUE FACULTATIVE 'SIMP' DE XMAT
26 : C AAD = COEFFICIENT INTERVENANT DS LE CALCUL DE LA VITESSE
27 : C DE LA DEFORMATION DE DENSIFICATION
28 : C GS(3) = RESISTANCES AU CISAILLEMENT
29 : C LEBIL(NC) = COMPRESSION/TRACTION
30 : C EPSPT(6) = VITESSE DES DEFORM. TOTALES SUR LE PAS D INTEGRATION
31 : C SIG0(6) = CONTRAINTES AU DEBUT DU SOUS PAS D INTEGRATION
32 : C EPSV0(6) = DEFORM. VISCOPLAST. AU DEBUT DU SOUS PAS D'INTEGRATION
33 : C VAR0(NGAT) = VAR. INT. SCAL. DE GATT_MONERIE AU DEB. DU SS PAS
34 : C W0(3) = OUVERTURES DE FISS. AU DEB. DU SS PAS D'INTEGRATION
35 : C WMAX0(3) = OUVERTURES MAXIMALES DES FISSURES AU DEB. DU SOUS PAS
36 : C XINVL(3) = PARAMETRES DE TAILLE
37 : C PENTE(NC) = PENTES DES CRITERES
38 : C DX0(NC) = DEF. DE FISSURATION (OUV.) AU DEB. DU SS PAS
39 : C NC1 = NC+1 AVEC NC(=3) NBR. TOTAL DE DIRECTIONS DE FISS.
40 : C NCA = NBR. DE DIRECTIONS DE FISS. OU UN CRITERE EST ATTEINT
41 : C NDIM = NCA+1 SI CP, NCA SINON
42 : C NN(NCA) = NUMEROS DES DIRECTIONS DE FISS. OU UN CRIT. EST ATTEINT
43 : C
44 : C SORTIES
45 : C -------
46 : C TAU = CONVERGENCE POUR t0+TAU
47 : C TAUNEX = ESTIMATION DE LA VALEUR SUIVANTE DE TAU
48 : C SIGF(6) = CONTRAINTES A t0+TAU
49 : C EPSVF(6) = DEFORM. VISCOPLAST. A t0+TAU
50 : C VARF(NGAT) = VAR. INT. SCAL. DE GATT_MONERIE A t0+TAU
51 : C WF(3) = OUVERTURES DE FISS. A t0+TAU
52 : C DXF(NC) = DEF. DE FISSURATION (OUV.) A t0+TAU
53 : C WMAXF(3) = OUVERTURES MAX. DES FISSURES A t0+TAU
54 : C (OUVERTURES + GLISSEMENTS) A t0+TAU
55 : C TF = TEMPERATURE A t0+TAU
56 : C FIF = DENSITE DE FISSION A t0+TAU
57 : C KERRE = GESTION DES ERREURS
58 : C----------------------------------------------------------------------
59 : C
60 : IMPLICIT INTEGER(I-N)
61 : IMPLICIT REAL*8(A-H,O-Z)
62 : -INC CCOPTIO
63 : C
64 : PARAMETER (XZER=0.D0,UNDEMI=.5D0,UN=1.D0,NC=3,NGATT=4)
65 : C
66 : DIMENSION EPSPT(*),SIG0(*),EPSV0(*),VAR0(*),W0(*),XINVL(*),DX0(*)
67 : DIMENSION WMAX0(*)
68 : DIMENSION SIGF(*),EPSVF(*),VARF(*),DXF(*),WF(*),WMAXF(*)
69 : DIMENSION PENTE(*),XMAT(*),NN(*),GS(*),LEBIL(*)
70 : DIMENSION DDE(18),DFDS(6,NC),DGDS(6,NC),HDFDQ(NC)
71 : DIMENSION EPSVP1(6),VARP1(NGATT),EPSDP1(3),EPSGP1(3),
72 : & DDX1(4),SIGP1(6)
73 : DIMENSION EPSVP2(6),VARP2(NGATT),EPSDP2(3),EPSGP2(3),
74 : & DDX2(4),SIGP2(6)
75 : DIMENSION EPSVP3(6),VARP3(NGATT),EPSDP3(3),EPSGP3(3),
76 : & DDX3(4),SIGP3(6)
77 : DIMENSION EPSVP4(6),VARP4(NGATT),EPSDP4(3),EPSGP4(3),
78 : & DDX4(4),SIGP4(6)
79 : DIMENSION SIG1(6),EPSV1(6),VAR1(NGATT),DX1(3),W1(3)
80 : DIMENSION SIG12(6),EPSV12(6),VAR12(NGATT),DX12(3),W12(3)
81 : DIMENSION SIG13(6),EPSV13(6),VAR13(NGATT),DX13(3),W13(3),XX(6)
82 : C
83 : C-------------------------------------
84 : C PARAMETRES POUR TESTS DE CONVERGENCE
85 : C-------------------------------------
86 : BORNE = 2.0
87 : RMAX = 1.3
88 : RMIN = 0.7
89 : DIV = 7.0
90 : FAC = 3.0
91 : DIV2 = 3.*DIV
92 : C
93 : C -------------------------
94 : C COEFFICIENTS D ELASTICITE
95 : C -------------------------
96 : YOUN=XMAT(1)
97 : XNU =XMAT(2)
98 : G = UNDEMI*YOUN/(UN+XNU)
99 : C
100 : C ---------------------------------------------------------------
101 : C COEF. POUR CALCULER LE BURNUP A PARTIR DE LA DENSITE DE FISSION
102 : C ---------------------------------------------------------------
103 : RHO=XMAT(3)
104 : RHO0=1.D0-XMAT(25)
105 : EFIS=XMAT(27)
106 : AMU238= 238.D0
107 : AMUO2 = 238.D0 + 16.D0 + 16.D0
108 : AKBU = (AMUO2/AMU238)*EFIS/(RHO*RHO0)
109 : C
110 : C---------------------------------
111 : C UTILE POUR CALCULER LA PRECISION
112 : C---------------------------------
113 : XMAX=YOUN*1.D-3
114 : C
115 : C ----------------
116 : C INITIALISATION
117 : C ----------------
118 : C
119 : KERRE=0
120 : C
121 : C --- Taux de combustion
122 : BU0=VAR0(2)
123 : C
124 : ASIG = SQRT(PROCON(SIG0,SIG0,6))
125 : ERRABS = PRECIS*ASIG
126 : IF (XMAX.GT.ASIG) ERRABS = PRECIS*XMAX
127 : C
128 : C
129 : C -- derivees 1
130 : BUPT=AKBU*FI0
130 : CALL UO2VI1(SIG0,EPSV0,VAR0,EPSVP1,VARP1,EPSDP1,EPSGP1,
131 : & XMAT,AAD,NGAT,NCOMAT,NSIMP,FI0,T0,BU0,BUPT)
131 : CALL OTTODD(DDE,YOUN,XNU,G,GS,XINVL,W0)
131 : CALL UO2FL(NCA,NN,XINVL,PENTE,SIG0,GS,DFDS,DGDS,HDFDQ,KERRE)
132 : IF (KERRE.NE.0) THEN
133 : GOTO 1100
134 : ENDIF
134 : CALL UO2VI2(NC1,NCA,NDIM,NN,PRECIZ,EPSPT,EPSVP1,EPSDP1,
135 : & EPSGP1,DDE,DFDS,DGDS,HDFDQ,DDX1,SIGP1,KERRE)
136 : IF (KERRE.NE.0) THEN
137 : GOTO 1100
138 : ENDIF
139 : C
140 : C
141 : NITERA = 0
142 : C -----------------------------------------------------------
143 : C DEBUT DES ITERATIONS SUR TAU OPTIMAL / FIN SI RA PETIT
144 : C -----------------------------------------------------------
145 : 1300 CONTINUE
146 : NITERA = NITERA + 1
147 : TAU2=TAU*UNDEMI
148 : C -- maj 1
148 : CALL UO2MAJ(IFOUR,NC,NCA,NDIM,NN,NGAT,XINVL,TAU2,SIG0,EPSV0,
149 : & VAR0,DX0,SIG1,EPSV1,VAR1,DX1,W1,SIGP1,EPSVP1,
150 : & VARP1,DDX1)
151 : C
152 : C -- derivees 2
153 : FI12=FI0+(TAU2*FPOINT)
154 : TI12=T0+(TAU2*TPOINT)
155 : BUPT=0.5D0*AKBU*(FI0+FI12)
156 : BU12=BU0+(TAU2*BUPT)
156 : CALL UO2VI1(SIG1,EPSV1,VAR1,EPSVP2,VARP2,EPSDP2,EPSGP2,
157 : & XMAT,AAD,NGAT,NCOMAT,NSIMP,FI12,TI12,BU12,BUPT)
157 : CALL OTTODD(DDE,YOUN,XNU,G,GS,XINVL,W1)
157 : CALL UO2FL(NCA,NN,XINVL,PENTE,SIG1,GS,DFDS,DGDS,HDFDQ,KERRE)
158 : IF (KERRE.NE.0) THEN
159 : GOTO 1100
160 : ENDIF
160 : CALL UO2VI2(NC1,NCA,NDIM,NN,PRECIZ,EPSPT,EPSVP2,EPSDP2,
161 : & EPSGP2,DDE,DFDS,DGDS,HDFDQ,DDX2,SIGP2,KERRE)
162 : IF (KERRE.NE.0) THEN
163 : GOTO 1100
164 : ENDIF
165 : C
166 : DO I=1,6
167 : SIGP2(I) =UNDEMI*(SIGP1(I)+SIGP2(I))
168 : EPSVP2(I)=UNDEMI*(EPSVP1(I)+EPSVP2(I))
169 : ENDDO
170 : VARP2(1)=UNDEMI*(VARP1(1)+VARP2(1))
171 : VARP2(2)=UNDEMI*(VARP1(2)+VARP2(2))
172 : DO I=1,NDIM
173 : DDX2(I)=UNDEMI*(DDX1(I)+DDX2(I))
174 : ENDDO
175 : C
176 : C -- maj 2
176 : CALL UO2MAJ(IFOUR,NC,NCA,NDIM,NN,NGAT,XINVL,TAU2,SIG0,EPSV0,
177 : & VAR0,DX0,SIG12,EPSV12,VAR12,DX12,W12,SIGP2,EPSVP2,
178 : & VARP2,DDX2)
179 : C ---------------------------------------------------------------------
180 : C Calcul du rapport : erreur calculee / erreur admise a mi-pas (TAU2)
181 : C pour eviter les demarrages trop brutaux avec une contrainte
182 : C elastique a t=0 tres grande si bien que la vitesse de deformation
183 : C viscoplastique diverge a la fin du pas defini par TAU
184 : C ---------------------------------------------------------------------
185 : IF(IIMPI.EQ.42) THEN
186 : WRITE(IOIMP,77002) NITERA
187 : *
188 : WRITE(IOIMP,18010) TAU,ERRABS
189 : 18010 FORMAT('0 UO2DC2 - TAU ERRABS '/2(1X,1PE12.5)/)
190 : *
191 : WRITE(IOIMP,74432) (SIG12(I),I=1,6)
192 : 74432 FORMAT('0 UO2DC2 - SIG12 ',6(1X,1PE12.5)/)
193 : WRITE(IOIMP,74434) (SIG1(I),I=1,6)
194 : *
195 : WRITE(IOIMP,74431) (DX12(I),I=1,3)
196 : 74431 FORMAT('0 UO2DC2 - DX12 ',3(1X,1PE12.5)/)
197 : WRITE(IOIMP,74436) (DX1(I),I=1,3)
198 : *
199 : WRITE(IOIMP,74430) (VAR12(I),I=1,NGAT)
200 : 74430 FORMAT('0 UO2DC2 - VAR12 ',3(1X,1PE12.5)/)
201 : WRITE(IOIMP,74438) (VAR1(I),I=1,NGAT)
202 : *
203 : WRITE(IOIMP,74429) (EPSV12(I),I=1,6)
204 : 74429 FORMAT('0 UO2DC2 - EPSV12 ',6(1X,1PE12.5)/)
205 : WRITE(IOIMP,74440) (EPSV1(I),I=1,6)
206 : ENDIF
207 : *
208 : DO 150 I=1,6
209 : XX(I) = SIG12(I)-SIG1(I)
210 : 150 CONTINUE
211 : RA = SQRT(PROCON(XX,XX,6))/ERRABS
212 : IF (RA.GT.DIV2*DIV2) THEN
213 : TAU = TAU/DIV2
214 : GOTO 1300
215 : ENDIF
216 : C
217 : C -- derivees 3
217 : CALL UO2VI1(SIG12,EPSV12,VAR12,EPSVP3,VARP3,EPSDP3,EPSGP3,
218 : & XMAT,AAD,NGAT,NCOMAT,NSIMP,FI12,TI12,BU12,BUPT)
218 : CALL OTTODD(DDE,YOUN,XNU,G,GS,XINVL,W12)
218 : CALL UO2FL(NCA,NN,XINVL,PENTE,SIG12,GS,DFDS,DGDS,HDFDQ,KERRE)
219 : IF (KERRE.NE.0) THEN
220 : GOTO 1100
221 : ENDIF
221 : CALL UO2VI2(NC1,NCA,NDIM,NN,PRECIZ,EPSPT,EPSVP3,EPSDP3,
222 : & EPSGP3,DDE,DFDS,DGDS,HDFDQ,DDX3,SIGP3,KERRE)
223 : IF (KERRE.NE.0) THEN
224 : GOTO 1100
225 : ENDIF
226 : C
227 : C -- maj 3
227 : CALL UO2MAJ(IFOUR,NC,NCA,NDIM,NN,NGAT,XINVL,TAU2,SIG12,
228 : & EPSV12,VAR12,DX12,SIG13,EPSV13,VAR13,DX13,W13,
229 : & SIGP3,EPSVP3,VARP3,DDX3)
230 : C
231 : C -- derivees 4
232 : FIF=FI0+(TAU*FPOINT)
233 : TF=T0+(TAU*TPOINT)
234 : BUPT=0.5D0*AKBU*(FI0+FIF)
235 : BUF=BU0+(TAU*BUPT)
235 : CALL UO2VI1(SIG13,EPSV13,VAR13,EPSVP4,VARP4,EPSDP4,EPSGP4,
236 : & XMAT,AAD,NGAT,NCOMAT,NSIMP,FIF,TF,BUF,BUPT)
236 : CALL OTTODD(DDE,YOUN,XNU,G,GS,XINVL,W13)
236 : CALL UO2FL(NCA,NN,XINVL,PENTE,SIG13,GS,DFDS,DGDS,HDFDQ,KERRE)
237 : IF (KERRE.NE.0) THEN
238 : GOTO 1100
239 : ENDIF
239 : CALL UO2VI2(NC1,NCA,NDIM,NN,PRECIZ,EPSPT,EPSVP4,EPSDP4,
240 : & EPSGP4,DDE,DFDS,DGDS,HDFDQ,DDX4,SIGP4,KERRE)
241 : IF (KERRE.NE.0) THEN
242 : GOTO 1100
243 : ENDIF
244 : C
245 : DO I=1,6
246 : SIGP4(I) =UNDEMI*(SIGP3(I)+SIGP4(I))
247 : EPSVP4(I)=UNDEMI*(EPSVP3(I)+EPSVP4(I))
248 : ENDDO
249 : VARP4(1)=UNDEMI*(VARP3(1)+VARP4(1))
250 : VARP4(2)=UNDEMI*(VARP3(2)+VARP4(2))
251 : DO I=1,NDIM
252 : DDX4(I)=UNDEMI*(DDX3(I)+DDX4(I))
253 : ENDDO
254 : C
255 : C -- maj 4
255 : CALL UO2MAJ(IFOUR,NC,NCA,NDIM,NN,NGAT,XINVL,TAU2,SIG12,
256 : & EPSV12,VAR12,DX12,SIGF,EPSVF,VARF,DXF,WF,SIGP4,EPSVP4,
257 : & VARP4,DDX4)
258 : IF (IFLAG.EQ.0) THEN
259 : GOTO 200
260 : ENDIF
261 : C
262 : C -- derivees 5
262 : CALL UO2VI1(SIGF,EPSVF,VARF,EPSVP4,VARP4,EPSDP4,EPSGP4,
263 : & XMAT,AAD,NGAT,NCOMAT,NSIMP,FIF,TF,BUF,BUPT)
263 : CALL OTTODD(DDE,YOUN,XNU,G,GS,XINVL,WF)
263 : CALL UO2FL(NCA,NN,XINVL,PENTE,SIGF,GS,DFDS,DGDS,HDFDQ,KERRE)
264 : IF (KERRE.NE.0) THEN
265 : GOTO 1100
266 : ENDIF
266 : CALL UO2VI2(NC1,NCA,NDIM,NN,PRECIZ,EPSPT,EPSVP4,EPSDP4,
267 : & EPSGP4,DDE,DFDS,DGDS,HDFDQ,DDX4,SIGP4,KERRE)
268 : IF (KERRE.NE.0) THEN
269 : GOTO 1100
270 : ENDIF
271 : C
272 : DO I=1,6
273 : SIGP2(I) =(SIGP1(I)+SIGP4(I))/6.d0+SIGP3(I)*2.d0/3.d0
274 : EPSVP2(I)=(EPSVP1(I)+EPSVP4(I))/6.d0+EPSVP3(I)*2.d0/3.d0
275 : ENDDO
276 : VARP2(1)=(VARP1(1)+VARP4(1))/6.d0+VARP3(1)*2.d0/3.d0
277 : VARP2(2)=(VARP1(2)+VARP4(2))/6.d0+VARP3(2)*2.d0/3.d0
278 : DO I=1,NDIM
279 : DDX2(I)=(DDX1(I)+DDX4(I))/6.d0+DDX3(I)*2.d0/3.d0
280 : ENDDO
281 : C
282 : C -- maj 5
282 : CALL UO2MAJ(IFOUR,NC,NCA,NDIM,NN,NGAT,XINVL,TAU,SIG0,EPSV0,
283 : & VAR0,DX0,SIG1,EPSV1,VAR1,DX1,W1,SIGP2,EPSVP2,
284 : & VARP2,DDX2)
285 : C ---------------------------------------------------------------------
286 : C CALCUL DU RAPPORT ERREUR CALCULEE/ERREUR ADMISE A LA FIN DU PAS (TAU)
287 : C ---------------------------------------------------------------------
288 : DO 170 I=1,6
289 : XX(I) = SIGF(I)-SIG1(I)
290 : 170 CONTINUE
291 : RA = SQRT(PROCON(XX,XX,6))/ERRABS
292 : SQRA = SQRT(RA)
293 : *
294 : IF(IIMPI.EQ.42) THEN
295 : WRITE(IOIMP,77002) NITERA
296 : 77002 FORMAT('0 UO2DC2 - NITERA ',I3/)
297 : *
298 : WRITE(IOIMP,18011) TAU,ERRABS,RA
299 : 18011 FORMAT('0 UO2DC2 - TAU ERRABS RA '/3(1X,1PE12.5)/)
300 : *
301 : WRITE(IOIMP,74433) (SIGF(I),I=1,6)
302 : 74433 FORMAT('0 UO2DC2 - SIGF ',6(1X,1PE12.5)/)
303 : WRITE(IOIMP,74434) (SIG1(I),I=1,6)
304 : 74434 FORMAT('0 UO2DC2 - SIG1 ',6(1X,1PE12.5)/)
305 : *
306 : WRITE(IOIMP,74435) (DXF(I),I=1,3)
307 : 74435 FORMAT('0 UO2DC2 - DXF ',3(1X,1PE12.5)/)
308 : WRITE(IOIMP,74436) (DX1(I),I=1,3)
309 : 74436 FORMAT('0 UO2DC2 - DX1 ',3(1X,1PE12.5)/)
310 : *
311 : WRITE(IOIMP,74437) (VARF(I),I=1,NGAT)
312 : 74437 FORMAT('0 UO2DC2 - VARF ',3(1X,1PE12.5)/)
313 : WRITE(IOIMP,74438) (VAR1(I),I=1,NGAT)
314 : 74438 FORMAT('0 UO2DC2 - VAR1 ',3(1X,1PE12.5)/)
315 : *
316 : WRITE(IOIMP,74439) (EPSVF(I),I=1,6)
317 : 74439 FORMAT('0 UO2DC2 - EPSVF ',6(1X,1PE12.5)/)
318 : WRITE(IOIMP,74440) (EPSV1(I),I=1,6)
319 : 74440 FORMAT('0 UO2DC2 - EPSV1 ',6(1X,1PE12.5)/)
320 : ENDIF
321 : C -----------------------------------------------
322 : C TEST DE FIN D'ITERATIONS / MISE A JOUR DE TAU
323 : C DIV =7 BORNE = 2
324 : C SI SQRA>7 TAU = TAU/7 ET NOUVEL ESSAI
325 : C SI 2RMAX on vise RA = 1
348 : C-------------------------------------
349 : TAUNEX=TAU
350 : IF (FAC*SQRA.LT.1.D0) THEN
351 : TAUNEX=TAU*FAC
352 : ELSEIF ( (SQRA.LT.RMIN).OR.(SQRA.GT.RMAX) ) THEN
353 : TAUNEX=TAU/SQRA
354 : ENDIF
355 : C
356 : 200 CONTINUE
357 : C
358 : C---------------------------------------------------------------
359 : C Maj du burnup - des ouvertures max. - des limites de fermeture
360 : C---------------------------------------------------------------
361 : VARF(2)=BUF
362 : C
363 : DO 400 I=1,3
364 : WMAXF(I)=WMAX0(I)
365 : IF (WF(I).GT.WMAX0(I)) THEN
366 : C IF(LEBIL(I).NE.1.OR.WMAX0(I).GE.WRUPT(I)) THEN
367 : WMAXF(I)=WF(I)
368 : C ENDIF
369 : ENDIF
370 : 400 CONTINUE
371 : C
372 : C
373 : 1100 CONTINUE
374 : C
375 : RETURN
376 : END
377 :
378 :
379 :
380 :
381 :
382 :
383 :