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