Télécharger excpha.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales