Télécharger excpha.eso

Retour à la liste

Numérotation des lignes :

  1. C EXCPHA SOURCE CB215821 19/11/07 21:15:01 10364
  2. subroutine excpha
  3.  
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8(A-H,O-Z)
  6.  
  7. -INC CCOPTIO
  8. -INC SMCHAML
  9. -INC SMCHPOI
  10. -INC SMMODEL
  11. -INC CCREEL
  12. -INC SMCOORD
  13. -INC SMELEME
  14.  
  15. segment icpr (NCPR)
  16. segment icpr1 (NCPR)
  17. segment icpr2 (NCPR)
  18.  
  19. SEGMENT SEXCP
  20. REAL*8 tdebu (node)
  21. REAL*8 tprop (node)
  22. REAL*8 qrest (node)
  23. REAL*8 qtot (node)
  24. REAL*8 qsecon(node)
  25. INTEGER iact (node)
  26. ENDSEGMENT
  27.  
  28. CHARACTER*4 MOT4a
  29. LOGICAL LOG1
  30.  
  31. C Un peu de marge sur la precision des calculs (pb aix)
  32. XZPREL=XZPREC*1.D2
  33.  
  34. C lecture du modele
  35. call lirobj('MMODEL ',mmodel,1,iretou)
  36. call actobj('MMODEL ',mmodel,1)
  37. if(ierr.ne.0) return
  38.  
  39. C lecture du materiau
  40. call lirobj('MCHAML ',IPIN,1,iretou)
  41. call actobj('MCHAML ',IPIN,1)
  42. if(ierr.ne.0) return
  43. CALL REDUAF(IPIN,mmodel,mchelm,0,IR,KER)
  44. IF(IR .NE. 1) CALL ERREUR(KER)
  45. IF(IERR .NE. 0) RETURN
  46.  
  47. C lecture des proportions de phase
  48. call lirobj('MCHAML ',IPIN,1,iretou)
  49. call actobj('MCHAML ',IPIN,1)
  50. if(ierr.ne.0) return
  51. CALL REDUAF(IPIN,mmodel,mchel1,0,IR,KER)
  52. IF(IR .NE. 1) CALL ERREUR(KER)
  53. IF(IERR .NE. 0) RETURN
  54.  
  55. C lecture des chaleur latentes reduies aux noeuds
  56. call lirobj('MCHAML ',IPIN,1,iretou)
  57. call actobj('MCHAML ',IPIN,1)
  58. if(ierr.ne.0) return
  59. CALL REDUAF(IPIN,mmodel,mchel2,0,IR,KER)
  60. IF(IR .NE. 1) CALL ERREUR(KER)
  61. IF(IERR .NE. 0) RETURN
  62.  
  63. C lecture du chpoint de temperatures initiales
  64. call lirobj('CHPOINT ',mchpoi,1,iretou)
  65. CALL ACTOBJ('CHPOINT ',mchpoi,1)
  66. if(ierr.ne.0) return
  67.  
  68. C lecture du chpoint d'increment de temperatures proposées
  69. call lirobj('CHPOINT ',mchpo1,1,iretou)
  70. CALL ACTOBJ('CHPOINT ',mchpo1,1)
  71. if(ierr.ne.0) return
  72.  
  73. C fin des lectures
  74. XUN =1.D0
  75. SEGACT,MCOORD
  76. NCPR=xcoor(/1)/(idim+1)
  77. SEGDES,MCOORD
  78. segini,icpr,icpr1,icpr2
  79.  
  80. C algorithmes on boucle sur les sous-modeles
  81. C puis si la rigidite avait ete prise en compte on teste le sens
  82. C et la valeur de la chaleur latente proposee
  83. C sinon on teste la non traverse de la temperature de changement de
  84. C phase et au besoin on active le blocage (temperature imposee et
  85. C flux de chaleur latente limite)
  86.  
  87. * +--------------------------------------------------------------------+
  88. * | Preparation de QSECOND : |
  89. * | CHPOINT de 'LX' |
  90. * | Remplace les chaleurs latentes dans RESO : UNILATER |
  91. * +--------------------------------------------------------------------+
  92. N1 = mmodel.kmodel(/1)
  93. ib = 0
  94. do 200 meri = 1,N1
  95. imodel = mmodel.kmodel(meri)
  96.  
  97. nfor=imodel.formod(/2)
  98. call place(imodel.formod,nfor,iplac,'CHANGEMENT_PHASE')
  99. if (iplac .eq. 0) goto 200
  100. if(tymode(2) .ne. 'MAILLAGE')then
  101. call erreur(5)
  102. endif
  103. ipt2=ivamod(2)
  104. do 201 nel=1,ipt2.num(/2)
  105. C Noeud 1 : 'LX'
  106. C Noeud 2 : 'Inconnues classiques'
  107. ia = ipt2.num(1,nel)
  108. if(icpr2(ia).eq.0) then
  109. ib = ib+1
  110. icpr2(ia)= ib
  111. endif
  112. 201 continue
  113. 200 continue
  114. node = ib
  115.  
  116. * +--------------------------------------------------------------------+
  117. * | NOEUDS du CHPOINT debu et proposé |
  118. * +--------------------------------------------------------------------+
  119. segini,icpr
  120. ib = 0
  121. do 80 I=1,mchpoi.ipchp(/1)
  122. msoupo=mchpoi.ipchp(i)
  123. C Le 'LX' du champ debut ne nous interesse pas
  124. if (msoupo.nocomp(1) .EQ. 'LX ') goto 80
  125. ipt1 =msoupo.igeoc
  126. mpoval=msoupo.ipoval
  127. do 81 j=1,ipt1.num(/2)
  128. ia = ipt1.num(1,j)
  129. if(icpr(ia).eq.0) then
  130. ib = ib +1
  131. icpr(ia)=ib
  132. endif
  133. 81 continue
  134. 80 continue
  135. node=max(node,ib)
  136.  
  137. ib = 0
  138. segini,icpr1
  139. do 90 I=1,mchpo1.ipchp(/1)
  140. msoupo=mchpo1.ipchp(i)
  141. ipt1 =msoupo.igeoc
  142. do 91 j=1,ipt1.num(/2)
  143. ia = ipt1.num(1,j)
  144. if(icpr1(ia).eq.0) then
  145. ib = ib +1
  146. icpr1(ia)= ib
  147. endif
  148. 91 continue
  149. 90 continue
  150. node=max(node,ib)
  151.  
  152. segini,SEXCP
  153.  
  154. * +--------------------------------------------------------------------+
  155. * | Boucle sur les SOUS-ZONES du MODELE |
  156. * +--------------------------------------------------------------------+
  157. do 100 mo=1,N1
  158. imodel = mmodel.kmodel(mo)
  159. meleme = imodel.imamod
  160.  
  161. nfor=imodel.formod(/2)
  162. call place(imodel.formod,nfor,iplac,'CHANGEMENT_PHASE')
  163. if (iplac .eq. 0) goto 100
  164.  
  165. * +------------------------------------------------------------------------+
  166. * | Recherche des composantes PRIMALE + 'LX' du CHPOINT debu et propose |
  167. * +------------------------------------------------------------------------+
  168. nomid=imodel.lnomid(1)
  169. MOT4a=nomid.lesobl(1)
  170.  
  171. C remise a zero qui va bien
  172. IF(mo .GT. 1)THEN
  173. CALL ZERO(tdebu(1),node,1)
  174. CALL ZERO(tprop(1),node,1)
  175. ENDIF
  176.  
  177. do 82 I=1,mchpoi.ipchp(/1)
  178. msoupo=mchpoi.ipchp(i)
  179. C Le 'LX' du champ debut ne nous interesse pas
  180. if (msoupo.nocomp(1) .EQ. 'LX ') goto 82
  181. call place(msoupo.nocomp,msoupo.nocomp(/2),iplac,MOT4a)
  182. if (iplac .eq. 0) then
  183. moterr = MOT4a
  184. call erreur(181)
  185. return
  186. endif
  187. ipt1 =msoupo.igeoc
  188. mpoval=msoupo.ipoval
  189. do 85 k=1,ipt1.num(/2)
  190. ia = ipt1.num(1,k)
  191. tdebu(icpr(ia))= mpoval.vpocha(k,iplac)
  192. 85 continue
  193. 82 continue
  194.  
  195. do 92 I=1,mchpo1.ipchp(/1)
  196. msoup1= mchpo1.ipchp(i)
  197. iplac = 1
  198. if (msoup1.nocomp(1) .EQ. 'LX ') goto 94
  199. call place(msoup1.nocomp,msoup1.nocomp(/2),iplac,MOT4a)
  200. if (iplac .eq. 0) then
  201. moterr = MOT4a
  202. call erreur(181)
  203. return
  204. endif
  205. 94 continue
  206. ipt1 = msoup1.igeoc
  207. mpova1= msoup1.ipoval
  208. do 95 k=1,ipt1.num(/2)
  209. ia = ipt1.num(1,k)
  210. if(icpr(ia) .EQ. 0)then
  211. C Cas des 'LX'
  212. tprop(icpr1(ia))= mpova1.vpocha(k,iplac)
  213. else
  214. tprop(icpr1(ia))= mpova1.vpocha(k,iplac) + tdebu(icpr(ia))
  215. endif
  216. 95 continue
  217. 92 continue
  218.  
  219. * +--------------------------------------------------------------+
  220. * | Recherche du MCHAML correspondant a 'PRIM' |
  221. * +--------------------------------------------------------------+
  222. do 51 mchm=1,imache(/1)
  223. if(conche(mchm) .ne. conmod) goto 51
  224. if(imache(mchm) .eq. meleme) then
  225. mchaml=ichaml(mchm)
  226. goto 52
  227. endif
  228. 51 continue
  229. call erreur(472)
  230. return
  231. 52 continue
  232.  
  233. do 56 n2=1,nomche(/2)
  234. if(nomche(n2).eq.'PRIM') then
  235. if (typche(n2).ne.'REAL*8') then
  236. moterr(1:16) = typche(n2)
  237. moterr(17:20) ='PPHA'
  238. moterr(21:36) = mchelm.titche
  239. call erreur(552)
  240. RETURN
  241. endif
  242. melval=ielval(n2)
  243. goto 57
  244. endif
  245. 56 continue
  246. moterr(1:8) = 'PRIM'
  247. call erreur(677)
  248. return
  249.  
  250. 57 continue
  251. C On suppose 'PRIM' constant par sous-modele : verification
  252. IF(velche(/1) .NE. 1 .AND. velche(/2).NE. 1)THEN
  253. MOTERR(1:4)='PRIM'
  254. L1 =mchelm.TITCHE(/1)
  255. L2 =MIN(20,L1)
  256. MOTERR(5:L2)=mchelm.TITCHE
  257. CALL ERREUR(1065)
  258. RETURN
  259. ENDIF
  260.  
  261. T_chph = melval.velche(1,1)
  262.  
  263. * +--------------------------------------------------------------+
  264. * | Recherche du MCHAML correspondant a 'PPHA' |
  265. * +--------------------------------------------------------------+
  266. do 61 mchm=1,mchel1.imache(/1)
  267. if(mchel1.conche(mchm) .ne. conmod) goto 61
  268. if(mchel1.imache(mchm) .eq. meleme) then
  269. mcham1=mchel1.ichaml(mchm)
  270. goto 62
  271. endif
  272. 61 continue
  273. call erreur(472)
  274. return
  275.  
  276. 62 continue
  277. do 66 n2=1,mcham1.nomche(/2)
  278. if(mcham1.nomche(n2).eq.'PPHA') then
  279. if (typche(n2).ne.'REAL*8') then
  280. moterr(1:16) = mcham1.typche(n2)
  281. moterr(17:20) ='PPHA'
  282. moterr(21:36) = mchel1.titche
  283. call erreur(552)
  284. RETURN
  285. endif
  286. melva1=mcham1.ielval(n2)
  287. goto 67
  288. endif
  289. 66 continue
  290. moterr(1:8) = 'PPHA'
  291. call erreur(677)
  292. return
  293. 67 continue
  294.  
  295. * +---------------------------------------------------------------------+
  296. * | Recherche du MCHAML correspondant au flux equivalent (aux NOEUDS)|
  297. * +---------------------------------------------------------------------+
  298. do 71 mchm=1,mchel2.imache(/1)
  299. if(mchel2.conche(mchm) .ne. conmod) goto 71
  300. if(mchel2. imache(mchm) .eq. meleme) then
  301. mcham2=mchel2.ichaml(mchm)
  302. goto 72
  303. endif
  304. 71 continue
  305. call erreur(472)
  306. return
  307. 72 continue
  308.  
  309. nomid = lnomid(2)
  310. do 76 n2=1,mcham2.nomche(/2)
  311. if(mcham2.nomche(n2).eq.nomid.lesobl(1)) then
  312. if (mcham2.typche(n2).ne.'REAL*8') then
  313. moterr(1:16) = mcham2.typche(n2)
  314. moterr(17:20) = nomid.lesobl(1)
  315. moterr(21:36) = mchel2.titche
  316. call erreur(552)
  317. RETURN
  318. endif
  319. melva2=mcham2.ielval(n2)
  320. goto 77
  321. endif
  322. 76 continue
  323. moterr(1:8) =nomid.lesobl(1)
  324. call erreur(677)
  325. return
  326. 77 continue
  327.  
  328. C remise a zero qui va bien
  329. IF(mo .GT. 1)THEN
  330. CALL ZERO(qrest(1),node,1)
  331. CALL ZERO(qtot(1) ,node,1)
  332. ENDIF
  333.  
  334. C Fabrication de qtot : chaleur latente totale cumulee par point
  335. C Fabrication de qrest: chaleur latente restante cumulee par point
  336. noMAX1=melva1.velche(/1)
  337. meMAX1=melva1.velche(/2)
  338. noMAX2=melva2.velche(/1)
  339. meMAX2=melva2.velche(/2)
  340. do 110 mel=1,meleme.num(/2)
  341. do 111 nod = 1,meleme.num(/1)
  342. ia = meleme.num(nod,mel)
  343. ib = icpr(ia)
  344. pp = melva1.velche(MIN(nod,noMAX1),MIN(mel,meMAX1))
  345. umpp = MAX(MIN(XUN-pp,XUN),XZERO)
  346. ql = melva2.velche(min(nod,noMAX2),min(mel,meMAX2))
  347. IF(ql .EQ. XZERO)THEN
  348. C Protection pour les chaleurs latentes strictement egales a 0.D0 !!
  349. ql = XPETIT/XZPREL
  350. ENDIF
  351. qr = ql*umpp
  352.  
  353. qrest(ib) = qrest(ib) + qr
  354. qtot (ib) = qtot (ib) + ql
  355. 111 continue
  356. 110 continue
  357.  
  358. ipt2 = imodel.ivamod(2)
  359. do 120 iel=1,ipt2.num(/2)
  360. C Noeud 1 : 'LX'
  361. C Noeud 2 : 'Inconnues classiques'
  362. No_LX = ipt2.num(1,iel)
  363. No_INC= ipt2.num(2,iel)
  364.  
  365. iex = icpr1(No_LX)
  366. iri = icpr2(No_LX)
  367. ideb = icpr(No_INC)
  368. ipro = icpr1(No_INC)
  369.  
  370. tdeb = tdebu(ideb)
  371. IF(ipro .NE. 0)THEN
  372. tpro = tprop(ipro)
  373. ELSE
  374. C Cas a la 1ere iteration
  375. tpro = tdeb
  376. ENDIF
  377.  
  378. q_tot = qtot(ideb)
  379. q_res = qrest(ideb)
  380. q_dej = q_tot-q_res
  381. q_pre = q_tot*XZPREL
  382.  
  383. C Pour eviter les soucis de denormalisation si T_chph ~ XZERO
  384. LOG1 =(tdeb.EQ.T_chph) .OR.
  385. & ABS(T_chph-tdeb).LT.MAX((XZPREL*T_chph),XPETIT)
  386. ibloc = 0
  387.  
  388. if(iex .ne. 0) then
  389. * +----------------------------------------------------------+
  390. * | Le 'LX' existe dans T proposée sur le noeud No_LX |
  391. * +----------------------------------------------------------+
  392. rea_LX= tprop(iex)
  393. if(LOG1) then
  394. C Temperature initiale = T_chph
  395. if (rea_LX.GT.XZERO .AND. q_res.GT.q_pre) then
  396. C la temperature cherche a monter
  397. C PRINT*,'-E:',No_LX,tdeb,T_chph,tpro
  398. ibloc =-1
  399. elseif(rea_LX.LT.XZERO .AND. q_dej.GT.q_pre) then
  400. C la temperature cherche a descendre
  401. C PRINT*,'-F:',No_LX,tdeb,T_chph,tpro
  402. ibloc = 1
  403. endif
  404.  
  405. elseif(rea_LX.GT.XZERO) then
  406. C La temperature cherche a monter
  407. if(tdeb.LT.T_chph) then
  408. C PRINT*,'-G:',No_LX,tdeb,tpro,q_tot,q_res,q_pre
  409. ibloc =-1
  410. endif
  411.  
  412. elseif(rea_LX.LT.XZERO) then
  413. C la temperature cherche a descendre
  414. if(tdeb .GT. T_chph) then
  415. C PRINT*,'-G:',No_LX,tdeb,tpro,q_tot,q_dej,q_pre
  416. ibloc = 1
  417. endif
  418. endif
  419.  
  420. else
  421. * +-------------------------------------------------------------+
  422. * | Le 'LX' n''existe pas dans T proposée sur le noeud No_LX |
  423. * +-------------------------------------------------------------+
  424. if(LOG1) then
  425. C Temperature initiale = T_chph
  426. tsens = tpro - T_chph
  427. if (tsens.GT.XZERO .AND. q_res.GT.q_pre) then
  428. C La temperature cherche a monter
  429. ibloc =-1
  430. C PRINT*,'-AA:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  431. elseif(tsens.LT.XZERO .AND. q_dej.GT.q_pre) then
  432. C la temperature cherche a descendre
  433. ibloc = 1
  434. C PRINT*,'-BB:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  435. endif
  436.  
  437. elseif((T_chph-tdeb) * (T_chph-tpro) .LT. XZERO) then
  438. C tdeb et tprop de chaque cote de T_chph : blocage actif
  439. if (tdeb .LT. T_chph) then
  440. ibloc =-1
  441. C PRINT*,'-C:',No_LX,(T_chph-tdeb),(T_chph-tpro)
  442. elseif(tdeb .GT. T_chph) then
  443. ibloc = 1
  444. C PRINT*,'-D:',No_LX,tdeb,T_chph,tpro
  445. else
  446. CALL ERREUR(5)
  447. endif
  448.  
  449. elseif(tdeb.LT.T_chph .AND. q_dej.GT.q_pre)then
  450. C apparemment il reste de la phase
  451. ibloc = 1
  452. C PRINT*,'-A:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  453.  
  454. elseif(tdeb.GT.T_chph .AND. q_res.GT.q_pre)then
  455. C apparemment on n'a pas tout consomme
  456. ibloc =-1
  457. C PRINT*,'-B:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  458. endif
  459. endif
  460.  
  461. C Les conditions mises sont unilaterales, RESO s'en sort en ne
  462. C gardant que la premiere 'excitee' dans le sens en question
  463. IF (ibloc.EQ. 1)THEN
  464. iact(iri) = No_LX
  465. qsecon(iri) = q_dej
  466. C PRINT*,'iact(iri)',No_LX,ibloc,iact(iri),qsecon(iri)
  467. ELSEIF(ibloc.EQ.-1)THEN
  468. iact(iri) = No_LX
  469. qsecon(iri) =-q_res
  470. C PRINT*,'iact(iri)',No_LX,ibloc,iact(iri),qsecon(iri)
  471. ELSE
  472. iact(iri) = 0
  473. qsecon(iri) = XZERO
  474. ENDIF
  475.  
  476. C fin de la boucle sur les elements du MAILLAGE de la rigidite
  477. 120 continue
  478. C fin de la boucle sur les sous-modeles
  479. 100 continue
  480.  
  481. * +--------------------------------------------------------------+
  482. * | Creation du resultat : CHPOINT de 'LX' |
  483. * | on dispose maintenant du tableau qsecon qui donne la |
  484. * | chaleur latente a mettre en 'LX' sur les noeuds de iact |
  485. * +--------------------------------------------------------------+
  486. nbelem = 0
  487. do 300 ia=1,node
  488. if(iact(ia) .ne. 0) nbelem = nbelem + 1
  489. 300 continue
  490.  
  491. nbref = 0
  492. nbsous = 0
  493. nat = 1
  494. if(nbelem .eq. 0) then
  495. nbnn = 0
  496. nsoupo = 0
  497. segini,mchpo3
  498.  
  499. else
  500. nbnn = 1
  501. segini,ipt4
  502. ipt4.itypel = 1
  503.  
  504. nsoupo = 1
  505. nc = 1
  506. n = nbelem
  507. segini,mchpo3,msoupo,mpoval
  508. mchpo3.mochde(1:72)='chpoint cree par EXCP'
  509. mchpo3.mtypoi ='contmult'
  510. mchpo3.jattri(1) = 2
  511. mchpo3.ipchp(1) = msoupo
  512. msoupo.nocomp(1) ='LX'
  513. msoupo.igeoc = ipt4
  514. msoupo.ipoval = mpoval
  515.  
  516. ipo=0
  517. do 301,ia=1,node
  518. iact1 = iact(ia)
  519. if (iact1 .eq. 0) goto 301
  520. ipo = ipo + 1
  521. ipt4.num(1,ipo) = iact1
  522. vpocha(ipo,1) = qsecon(icpr2(iact1))
  523. 301 continue
  524. endif
  525.  
  526. mchpo3.mochde(1:72)='chpoint cree par EXCP'
  527. mchpo3.mtypoi ='contmult'
  528. mchpo3.jattri(1) = 2
  529.  
  530. segsup icpr,icpr1,icpr2,SEXCP
  531.  
  532. call actobj('CHPOINT ',mchpo3,1)
  533. call ecrobj('CHPOINT ',mchpo3)
  534.  
  535. end
  536.  
  537.  

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