1 : C EXCPHA SOURCE CHAT 11/03/16 21:21:36 6902 2 : subroutine excpha 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : -INC CCOPTIO 6 : -INC SMCHAML 7 : -INC SMCHPOI 8 : -INC SMMODEL 9 : -INC SMRIGID 10 : -INC CCREEL 11 : -INC SMCOORD 12 : -INC SMELEME 13 : segment icpr(xcoor(/1)/(idim+1)) 14 : segment icpr1(xcoor(/1)/(idim+1)) 15 : segment icpr2(xcoor(/1)/(idim+1)) 16 : segment tdeb(node) 17 : segment tprop(node) 18 : segment qrest(node) 19 : segment qtot(node) 20 : segment qsecon(node) 21 : segment timpo(node) 22 : segment iact(node) 23 : segment prop(node) 24 : * lecture du modele 24 : call lirobj('MMODEL',mmodel,1,iretou) 25 : if(ierr.ne.0) return 26 : * lecture du materiau ( pour les TPHAS) 26 : call lirobj('MCHAML',mchelm,1,iretou) 27 : if(ierr.ne.0) return 28 : * lecture des proportions de phase 28 : call lirobj('MCHAML',mchel1,1,iretou) 29 : if(ierr.ne.0) return 30 : * lecture des chaleur latentes reduies aux noeuds 30 : call lirobj('MCHAML',mchel2,1,iretou) 31 : if(ierr.ne.0) return 32 : * lecture du chpoint de temperatures initiales 32 : call lirobj('CHPOINT',mchpoi,1,iretou) 33 : if(ierr.ne.0) return 34 : * lecture du chpoint de temperatures proposées 34 : call lirobj('CHPOINT',mchpo1,1,iretou) 35 : * write(6,*) ' entree dans excpha temp proposee' 35 : * call ecchpo(mchpo1) 36 : if(ierr.ne.0) return 37 : * lecture des rigidites de blocages de phases. 37 : call lirobj('RIGIDITE',mrigid,1,iretou) 38 : segact mrigid 39 : mpt = irigel(1,1) 40 : * write(6,*) ' maillage sous tendant les rigidites' 40 : * call ecmail ( mpt) 41 : if(ierr.ne.0) return 42 : * fin des lectures 43 : xpre=1.d-8 44 : unxpre=1.D0-xpre 45 : * 46 : * algorithmes on boucle sur les rigidite s de blocages ou sur les 47 : * modeles ( ils sont de structure identiques par construction des 48 : * rigidites) 49 : * puis si la rigidite avait ete prise en compte on teste le sens 50 : * et la valeur de la chaleur latente proposee 51 : * sinon on teste la non traverse de la temperature de changement de 52 : * phase et au besoin on active le blocage ( temperature imposee et 53 : *flux de chaleur latente limite) 54 : * 55 : * Au préalable repérage des champs de temperatures initiales et 56 : * proposees d'abord champs de temperatures debut on melange les 57 : * multiplicateurs 58 : 59 : segini icpr 60 : segact mchpoi 61 : ib = 0 62 : do 80 I=1,ipchp(/1) 63 : msoupo=ipchp(i) 64 : segact msoupo 65 : ipt1=igeoc 66 : segact ipt1 67 : do 81 j=1,ipt1.num(/2) 68 : ia = ipt1.num(1,j) 69 : if(icpr(ia).eq.0) then 70 : ib = ib +1 71 : icpr(ia)=ib 72 : endif 73 : 81 continue 74 : 80 continue 75 : node=ib 76 : segini tdeb 77 : do 82 I = 1 , ipchp(/1) 78 : msoupo=ipchp(i) 79 : ipt1=igeoc 80 : mpoval=ipoval 81 : segact mpoval 82 : j = 1 83 : if(nocomp(/2) .ne.1) then 84 : do 83 j=1,nocomp(/2) 85 : if( nocomp(J).eq.'T') go to 84 86 : 83 continue 87 : moterr(1:4)=' T' 87 : call erreur (181) 88 : return 89 : 84 continue 90 : endif 91 : do 85 k=1,ipt1.num(/2) 92 : ia = ipt1.num(1,k) 93 : ib = icpr(ia) 94 : tdeb(ib)=vpocha(k,j) 95 : 85 continue 96 : 82 continue 97 : * maintenant champs de temperatures proposees 98 : segact mchpo1 99 : ib = 0 100 : segini icpr1 101 : do 90 I=1,mchpo1.ipchp(/1) 102 : msoupo=mchpo1.ipchp(i) 103 : segact msoupo 104 : ipt1=igeoc 105 : segact ipt1 106 : do 91 j=1,ipt1.num(/2) 107 : ia = ipt1.num(1,j) 108 : if(icpr1(ia).eq.0) then 109 : ib = ib +1 110 : icpr1(ia)=ib 111 : endif 112 : 91 continue 113 : 90 continue 114 : node=max(node,ib) 115 : mnode=node 116 : segini tprop 117 : do 92 I = 1 , mchpo1.ipchp(/1) 118 : msoupo=mchpo1.ipchp(i) 119 : ipt1=igeoc 120 : mpoval=ipoval 121 : segact mpoval 122 : j = 1 123 : if(nocomp(/2) .ne.1) then 124 : do 93 j=1,nocomp(/2) 125 : if( nocomp(J).eq.'T') go to 94 126 : 93 continue 127 : moterr(1:4)=' T' 127 : call erreur (181) 128 : return 129 : 94 continue 130 : endif 131 : do 95 k=1,ipt1.num(/2) 132 : ia = ipt1.num(1,k) 133 : ib = icpr1(ia) 134 : tprop(ib)=vpocha(k,j) 135 : 95 continue 136 : 92 continue 137 : * au prealable encore, fabrication et reperage de Qsecon qui est le 138 : * champ de LX passé a resou( remplace les chaleurs latentes) 139 : segact mmodel,mrigid,mchelm,mchel1,mchel2 140 : segini icpr2 141 : ib=0 142 : do 200 meri = 1,irigel(/2) 143 : ipt3=irigel(1,meri) 144 : segact ipt3 145 : do 201 nel=1,ipt3.num(/2) 146 : do 201 nod = 1,1 147 : ia = ipt3.num(nod,nel) 148 : if(icpr2(ia).eq.0) then 149 : ib=ib+1 150 : icpr2(ia)=ib 151 : endif 152 : 201 continue 153 : 200 continue 154 : node = ib 155 : segini qsecon,timpo,iact 156 : * 157 : * debut du travail on boucle sur les sous zones du modeles 158 : * 159 : do 100 mo=1,kmodel(/1) 160 : imodel = kmodel(mo) 161 : segact imodel 162 : meleme = imamod 163 : segact meleme 164 : * recherche des chamelems correspondants 165 : * d'abord le champ de materiau recherche de TPHA 166 : do51 mchm=1,imache(/1) 167 : if( imache(mchm) . eq. meleme) then 168 : mchaml=ichaml(mchm) 169 : go to 52 170 : endif 171 : 51 continue 171 : call erreur ( 472) 172 : return 173 : 52 continue 174 : segact mchaml 175 : do 56 n2=1,nomche(/2) 176 : if( nomche(n2).eq.'TPHA') then 177 : melval=ielval(n2) 178 : go to 57 179 : endif 180 : 56 continue 181 : moterr(1:8) = 'TPHA' 181 : call erreur ( 677) 182 : return 183 : 57 continue 184 : * recherche de la proportion de phase PPHA 185 : do61 mchm=1,mchel1.imache(/1) 186 : if(mchel1. imache(mchm) . eq. meleme) then 187 : mcham1=mchel1.ichaml(mchm) 188 : go to 62 189 : endif 190 : 61 continue 190 : call erreur ( 472) 191 : return 192 : 62 continue 193 : segact mcham1 194 : do 66 n2=1,mcham1.nomche(/2) 195 : if(mcham1. nomche(n2).eq.'PPHA') then 196 : melva1=mcham1.ielval(n2) 197 : go to 67 198 : endif 199 : 66 continue 200 : moterr(1:8) = 'PPHA' 200 : call erreur ( 677) 201 : return 202 : 67 continue 203 : * recherche des chaleurs latentes reduites aux noeuds 204 : do71 mchm=1,mchel2.imache(/1) 205 : if(mchel2. imache(mchm) . eq. meleme) then 206 : mcham2=mchel2.ichaml(mchm) 207 : go to 72 208 : endif 209 : 71 continue 209 : call erreur ( 472) 210 : return 211 : 72 continue 212 : segact mcham2 213 : do 76 n2=1,mcham2.nomche(/2) 214 : if(mcham2. nomche(n2).eq.'Q ') then 215 : melva2=mcham2.ielval(n2) 216 : go to 77 217 : endif 218 : 76 continue 219 : moterr(1:8) = 'QLAT' 219 : call erreur ( 677) 220 : return 221 : 77 continue 222 : segact melval,melva1,melva2 223 : * par sous-modele on suppose TPHA constant. 224 : * on extrait la Temperature de changement de phase. 225 : * fabrication de qtot : chaleur latente additionnée par point 226 : * A l'aide de PPHA et de QLAR on fabrique un QREST qui est la chaleur 227 : * latente qui reste à consommer pour ce sous modele. 228 : * 229 : node=mnode 230 : segini qrest,qtot,prop 231 : do 110 mel=1,num(/2) 232 : do 110 nod = 1,num(/1) 233 : ia = num(nod,mel) 234 : ib = icpr(ia) 235 : pp = melva1.velche(nod,mel) 236 : if( pp.gt.1.d0) then 237 : * write(6,*) ' excp : erreur 110 1' 238 : pp = 1.d0 239 : endif 240 : if( pp. lt.0.d0) then 241 : * write(6,*) ' excp : erreur 110 2' 242 : pp = 0.D0 243 : endif 244 : prop(ib)=pp 245 : ql=melva2.velche(nod,mel) 246 : qr=ql* ( 1.d0 - pp) 247 : qrest(ib) = qrest(ib) + qr 248 : qtot(ib) = qtot(ib) + ql 249 : prop(ib)=1.d0 - qrest(ib)/qtot(ib) 250 : 110 continue 251 : tt = velche(1,1) 252 : ipt2= irigel(1,mo) 253 : segact ipt2 254 : do 120 meri=1,ipt2.num(/2) 255 : ip = ipt2.num(2,meri) 256 : il1=ipt2.num(1,meri) 257 : * ce blocage existe -t-il? 258 : iex= icpr1(il1) 259 : tde = tdeb(icpr(ip)) 260 : * tpro est la valeur proposée 261 : tpro=tprop(icpr1(ip)) 262 : if(iex.ne.0) then 263 : * if (iimpi.eq.1947) then 264 : * write(6,*) 'ip iex qtot(icpr(ip)',ip,iex,qtot(icpr(ip)) 265 : * write(6,*) 'ip qres prop' ,qrest(icpr(ip)),prop(icpr(ip)) 266 : * write(6,*) 'tde tpro ' , tde,tpro 267 : * write(6,*) 'il1 tprop(iex) ', il1,tprop(iex) 268 : * write(6,*) ' tt - tde ', tt-tde 269 : * write(6,*) 'tprop(iex)*qtot(icpr(ip)',tprop(iex)*qtot(icpr(ip)) 270 : * endif 271 : * le blocage existe on regarde si le point etait a la temperature de 272 : * changement de phase dans ce cas il faut garder le blocage si pas tout 273 : * deja dans un etat 274 : if( abs(tt-tde)/max(abs(tt),1.d0).le.1.e-8 ) then 275 : * if( iimpi.eq.1947) then 276 : * write(6,*) ' detection d un point a la temp de changement' 277 : * endif 278 : if( (tprop(iex).gt.0.d0.and . prop(icpr(ip)).lT.unxpre).OR. 279 : $ (tprop(iex).lt.0.D0. and. prop(icpr(ip)). gt.xpre)) then 280 : * write(6,*) ' on garde la condition ' 281 : ipl1=icpr2(il1) 282 : iact(ipl1)=il1 283 : if(tprop(iex)*qtot(icpr(ip)).gt.0.D0) then 284 : * la tempe cherche à monter 285 : qsecon(ipl1)=-qrest(icpr(ip))+qsecon(ipl1) 286 : else 287 : qsecon(ipl1)=qtot(icpr(ip))-qrest(icpr(ip))+qsecon(ipl1) 288 : endif 289 : endif 290 : * maintenant il faut regarder si le blocage est dans le bon sens 291 : elseif(tprop(iex)*qtot(icpr(ip)).gt.0.D0) then 292 : * la temperature cherche à monter 293 : * la temperature initiale etait-elle superieure a tt ? 294 : * si non il faut prendre -qrest si oui relacher la condition 295 : if(tde . le . tt) then 296 : * if( iimpi.eq.1947) then 297 : * write(6,*) 'on garde le blocage' 298 : * endif 299 : ipl1=icpr2(il1) 300 : qsecon(ipl1)=-qrest(icpr(ip))+qsecon(ipl1) 301 : * on impose aussi la temperature 302 : iact(ipl1)=il1 303 : endif 304 : else 305 : * la temperature cherche a descendre 306 : * la temperature initiale etait-elle inferieure a tt? 307 : * si oui relacher la condition sinon prendre qtot - qrest 308 : if(tde.ge. tt) then 309 : * if(iimpi.eq.1947) then 310 : * write(6,*) ' temperature chute on solidifie' 311 : * endif 312 : ipl1=icpr2(il1) 313 : qsecon(ipl1)=qtot(icpr(ip))-qrest(icpr(ip))+qsecon(ipl1) 314 : * on impose aussi la temperature 315 : iact(ipl1)=il1 316 : endif 317 : endif 318 : else 319 : * le blocage n'existe pas 320 : * verification que tdeb et tprop son du meme cote que tt 321 : * si c'est vrai rien a faire sinon impose le blocage et +ou- qtot 322 : aa = (tt - tde)* ( tt - tpro) 323 : tsens= tpro-tde 324 : if( abs((tt-tde)/(tt+1.d-30)).le.1.e-9 ) aa=0.d0 325 : if(aa.eq.0.D0) then 326 : * if(aa.le.0.D0) then 327 : * write(6,*) 'ip iex qtot(icpr(ip)',ip,iex,qtot(icpr(ip)) 328 : * write(6,*) ' tt - tde ', tt-tde 329 : * write(6,*) 'tde tpro aa ', tde,tpro, aa 330 : * write(6,*) 'tt - tde' , 'tt - tpro',tt - tde,tt - tpro 331 : * uuut=tprop(iex) 332 : * write(6,*)' tprop(iex)' , uuut 333 : * endif 334 : if( (tsens.gt.0.d0.and . prop(icpr(ip)).lT.unxpre).OR. 335 : $ (tsens.lt.0.D0. and. prop(icpr(ip)). gt.xpre)) then 336 : * write(6,*) ' on passe ici' 337 : * if( iimpi.eq.1947) then 338 : * write(6,*) ' un point changeait de phase' 339 : * endif 340 : xde=1.d0 341 : if( tpro.gt.tt*1.0000001) xde=-1.D0 342 : ipl1=icpr2(il1) 343 : 344 : * if( iimpi.eq.1947) then 345 : * write(6,*) ' on garde le point ip ',ip 346 : * write(6,*) '2 il1,il2,ipl1,ipl2 xde',il1,il2,ipl1,ipl2,xde 347 : * endif 348 : qsecon(ipl1)=qtot(icpr(ip))* xde+qsecon(ipl1) 349 : iact(ipl1)=il1 350 : * write(6,*) ' qsecon (ipl1), ip,qtot(icpr(ip))', 351 : * $ qsecon (ipl1), ip,qtot(icpr(ip)) 352 : endif 353 : elseif( aa . lt . 0.D0) then 354 : * tous les deux de chaque cote impose le blocage et la temperature 355 : * impose 356 : xde=1.d0 357 : if((tt-tde*1.000000001). GT.0.) xde = -1.d0 358 : ipl1=icpr2(il1) 359 : * if( iimpi.eq.1947) then 360 : * write(6,*) 'il1,il2,ipl1,ipl2 xde',il1,il2,ipl1,ipl2,xde 361 : * endif 362 : qsecon(ipl1)=qtot(icpr(ip))* xde+qsecon(ipl1) 363 : iact(ipl1)=il1 364 : * write(6,*) 'qtot qsec',qtot(icpr(ip)),qsecon(ipl1) 365 : endif 366 : endif 367 : * fin de la boucle sur les rigidites de blocages 368 : 120 continue 369 : * fin de la boucle sur les sous paquets de rigidite 370 : 100 continue 371 : * on dispose maintenant d'un tableau qsecon et timpo et iact qui 372 : * donnent 373 : * les second membre de chaleur latentes a mettre en LX et les 374 : * temperatures de fusion a mettre en FLX 375 : * on cree le champoint 376 : nbelem=0 377 : * write(6,*) (iact(ij),ij=1,iact(/1)) 378 : do 300 ia=1,iact(/1) 379 : if(iact(ia).ne.0) nbelem = nbelem + 1 380 : 300 continue 381 : nbsous=0 382 : nbref=0 383 : nbnn=1 384 : if(nbelem.eq.0) nbnn=0 385 : segini ipt4 386 : ipt4.itypel=1 387 : nat=1 388 : * write(6,*) ' nbnn nbelem ' , nbnn,nbelem 389 : IF(nbelem.ne.0) then 390 : nsoupo=1 391 : segini mchpo3 392 : mchpo3.mochde(1:72)=' creation par excpha' 393 : mchpo3.mtypoi='contmult' 394 : nc=1 395 : segini msoupo 396 : mchpo3. jattri(1)=2 397 : mchpo3.ipchp(1)=msoupo 398 : nocomp(1) = 'LX' 399 : igeoc=ipt4 400 : n=nbelem 401 : segini mpoval 402 : ipoval = mpoval 403 : ipo=0 404 : do 301,ia=1,iact(/1) 405 : if (iact(ia).ne.0) then 406 : ipo=ipo+1 407 : ipt4.num(1,ipo) = iact(ia) 408 : vpocha(ipo,1)=qsecon(icpr2(iact(ia))) 409 : endif 410 : 301 continue 411 : * write(6,*) ' nbelem ipo', nbelem ,ipo 412 : segdes mpoval,ipt4,msoupo,mchpo3 413 : else 414 : nsoupo=0 415 : segini mchpo3 416 : mchpo3.mochde(1:72)='creation par excpha ' 417 : mchpo3.jattri(1)=2 418 : mchpo3.mtypoi='contmult' 419 : segdes mchpo3 420 : endif 420 : call ecrobj('CHPOINT',mchpo3) 421 : * write(6,*) ' ecriture du chpoint, dans excpha, aprs ecrobj' 421 : * call ecchpo ( mchpo3) 422 : segdes melvaL,melva1,melva2 423 : do 303 ir=1,irigel(/2) 424 : meleme = irigel(1,ir) 425 : segdes meleme 426 : imodel = kmodel(ir) 427 : meleme = imamod 428 : segdes meleme,imodel 429 : 303 continue 430 : segdes mrigid,mmodel 431 : segdes mchelm,mchel1,mchel2 432 : * desactivation champ de temperature initiale 433 : do 304 ia=1,ipchp(/1) 434 : msoupo=ipchp(ia) 435 : meleme = igeoc 436 : segdes meleme 437 : mpoval = ipoval 438 : segdes mpoval 439 : segdes msoupo 440 : 304 continue 441 : segdes,mchpoi 442 : * desactivation champ de temperature propose 443 : IF (mchpoi.NE.mchpo1) THEN 444 : mchpoi= mchpo1 445 : do 305 ia=1,ipchp(/1) 446 : msoupo=ipchp(ia) 447 : meleme = igeoc 448 : segdes meleme 449 : mpoval = ipoval 450 : segdes mpoval 451 : segdes msoupo 452 : 305 continue 453 : segdes,mchpoi 454 : ENDIF 455 : segsup icpr,icpr1,icpr2,tdeb,qsecon,timpo,tprop,iact 456 : segsup qrest,qtot,prop 457 : return 458 : end 459 : 460 : 461 : 462 :
© Cast3M 2003 - All rights reserved.
Disclaimer