Télécharger excfro.eso

Retour à la liste

Numérotation des lignes :

  1. C EXCFRO SOURCE PV 16/06/03 21:15:03 8942
  2.  
  3. SUBROUTINE EXCFRO
  4.  
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7.  
  8. * on lit un modele de contact qui contient des frottement
  9. * un materiau de frottement, un chpoint
  10. * de deplacement, la raideur de frottement et on ecrit un chpoint de
  11. * forces de frottement
  12.  
  13. -INC CCOPTIO
  14. -INC CCREEL
  15.  
  16. -INC SMCHAML
  17. -INC SMCHPOI
  18. -INC SMCOORD
  19. -INC SMELEME
  20. -INC SMMODEL
  21. POINTEUR mmode3.mmodel,mmode4.mmodel,mmode5.mmodel,mmode6.mmodel
  22. POINTEUR mmode7.mmodel,mmode8.mmodel,imode6.imodel,imode7.imodel
  23.  
  24.  
  25. -INC SMRIGID
  26.  
  27. * icpr xcpr lx du deplacement
  28. * jcpr ycpr flx du frottement
  29. SEGMENT icpr(nbtpts)
  30. SEGMENT xcpr(nbtpts)
  31. SEGMENT jcpr(nbtpts)
  32. SEGMENT ycpr(nbtpts)
  33. SEGMENT kcpr(nbtpts)
  34. SEGMENT zcpr(nbtpts)
  35. SEGMENT icablr(nbtpts)
  36. *
  37. * Segments pour le frottement cable
  38. *
  39. SEGMENT siezo
  40. INTEGER iezon(nszo)
  41. ENDSEGMENT
  42. C
  43. SEGMENT altrav
  44. C*NU* REAL*8 ANG(NAM+1),ACUR(NAM+1)
  45. REAL*8 DANG(NAM),DACUR1(NAM),DACUR2(NAM)
  46. ENDSEGMENT
  47. C
  48. SEGMENT mwrk3
  49. C*NU* REAL*8 X1(3),X2(3),X3(3)
  50. REAL*8 RL(3),RS(3)
  51. ENDSEGMENT
  52. C
  53. SEGMENT icprg(nbtpts)
  54. C*NU* SEGMENT icprl(nbtpts)
  55. C
  56. SEGMENT idmul(nbnds)
  57. C*NU* SEGMENT idcp(nbnds)
  58. C*NU* SEGMENT idcpl(nbndsl)
  59.  
  60. C=======================================================================
  61. C= LECTURE DES DONNEES DE L'OPERATEUR
  62. C=======================================================================
  63. call lirobj('MMODEL ',modini,1,iretou)
  64. if (ierr.ne.0) return
  65. call lirobj('MCHAML ',icham1,1,iretou)
  66. if (ierr.ne.0) return
  67. call lirobj('CHPOINT ',mchpin,1,iretou)
  68. if (ierr.ne.0) return
  69. call lirobj('RIGIDITE',mrigid,1,iretou)
  70. if (ierr.ne.0) return
  71.  
  72. icham2 = 0
  73. call lirobj('MCHAML ',icham2,0,iretou)
  74. if (ierr.ne.0) return
  75. mchpi2 = 0
  76. call lirobj('CHPOINT ',mchpi2,0,iretou)
  77. if (ierr.ne.0) return
  78.  
  79. C- Tri des MCHAML donnes en entree
  80. CALL rngcha(icham1,icham2,'CARACTERISTIQUES','CONTRAINTES',
  81. & ipcar,ipcon)
  82. if (ierr.ne.0) return
  83. C- Le champ de caracteristiques est toujours obligatoire !
  84. if (ipcar.eq.0) then
  85. moterr(1:16)='CARACTERISTIQUES'
  86. call erreur(565)
  87. return
  88. endif
  89.  
  90. C=======================================================================
  91. C= QUELQUES INITIALISATIONS
  92. C=======================================================================
  93. idimp1 = IDIM + 1
  94. nbtpts = XCOOR(/1) / idimp1
  95. * write(ioimp,*) ' Entree dans excfro jcpr(/1) ',nbtpts
  96.  
  97. icablr = 0
  98. mwrk3 = 0
  99. ** ymcrit1 = 0.d0
  100.  
  101. C=======================================================================
  102. C= ANALYSE DU MODELE DU FROTTEMENT et DECOUPAGE EN MODELES DE MEME TYPE
  103. C=======================================================================
  104. mmodel=modini
  105. SEGACT,mmodel
  106. nsoum = kmodel(/1)
  107.  
  108. n1 = nsoum
  109. icable = 0
  110. SEGINI,mmode1
  111. mocabl = mmode1
  112. icoulo = 0
  113. SEGINI,mmode2
  114. mocoul = mmode2
  115.  
  116. DO isoum = 1, nsoum
  117. imodel = kmodel(isoum)
  118. SEGACT,imodel
  119. * write(ioimp,*) isoum,' formod(1) matmod(1)',formod(1),matmod(1)
  120. IF (formod(1).ne.'CONTACT') then
  121. moterr(1:16)='CONTACT'
  122. call erreur(719)
  123. goto 9999
  124. ENDIF
  125. icabl=0
  126. icoul=0
  127. do inma=1,matmod(/2)
  128. if ( matmod(inma).eq.'FROCABLE' ) icabl=1
  129. if ( matmod(inma).eq.'COULOMB') icoul=1
  130. enddo
  131. IF (icabl.eq.1) THEN
  132. * IF (idim.eq.2.and.NEFMOD.NE.261) call erreur(16)
  133. * IF (idim.eq.3.and.NEFMOD.NE.262) call erreur(16)
  134. * IF (tymode(1).NE.'MMODEL') call erreur(975)
  135. icable = icable + 1
  136. mmode1.kmodel(icable) = imodel
  137. ELSE IF (icoul.eq.1) THEN
  138. * IF (idim.eq.2.and.NEFMOD.NE.107) call erreur(16)
  139. * IF (idim.eq.3.and.NEFMOD.NE.165) call erreur(16)
  140. icoulo = icoulo + 1
  141. mmode2.kmodel(icoulo) = imodel
  142. ELSE
  143. write(ioimp,*) 'Modele de FROTTEMENT non implemente'
  144. call erreur(21)
  145. ENDIF
  146. IF (ierr.ne.0) goto 9999
  147. ENDDO
  148. ** write(6,*) ' icable icoulo ' , icable ,icoulo
  149. n1 = icable
  150. SEGADJ,mmode1
  151. mocabl = mmode1
  152. n1 = icoulo
  153. SEGADJ,mmode2
  154. mocoul = mmode2
  155.  
  156. C- Le champ de contraintes est obligatoire pour les modeles FROCABLE !
  157. if (icable.ne.0 .and. ipcon.eq.0) then
  158. moterr(1:16)='CONTRAINTES '
  159. call erreur(565)
  160. goto 9999
  161. endif
  162.  
  163. C=======================================================================
  164. * noter les contacts actifs dans le champ de deplacement
  165. C=======================================================================
  166. segact mcoord
  167. segini icpr,xcpr,jcpr,ycpr
  168. if (IDIM.eq.3) segini zcpr,kcpr
  169. segini,icablr
  170.  
  171. * call ecchpo(mchpin,0)
  172. mchpo1=mchpin
  173. segact mchpo1
  174. do 100 isoupo=1,mchpo1.ipchp(/1)
  175. msoupo=mchpo1.ipchp(isoupo)
  176. segact msoupo
  177. do 110 ic=1,nocomp(/2)
  178. if (nocomp(ic).ne.'LX ') go to 110
  179. ipt2=igeoc
  180. mpoval=ipoval
  181. segact ipt2,mpoval
  182. nbel2 = ipt2.num(/2)
  183. do 130 iel=1,nbel2
  184. impt = ipt2.num(1,iel)
  185. icpr(impt)=1
  186. xcpr(impt)=xcpr(impt)+vpocha(iel,ic)
  187. 130 continue
  188. segdes ipt2,mpoval
  189. 110 continue
  190. segdes msoupo
  191. 100 continue
  192. segdes mchpo1
  193.  
  194. C=======================================================================
  195. C= TRAITEMENT DES MODELES DE TYPE 'FROCABLE'
  196. C=======================================================================
  197. IF (icable.EQ.0) GOTO 1000
  198.  
  199. C --- Modèle global de frottement des cables
  200. * write(ioimp,*) 'Modèle global: ',mocabl
  201. * write(ioimp,*)'Excfro - Nbre de sous modèles: ',icable
  202. mmode1 = mocabl
  203.  
  204. C --- Renumerotation locale sur les noeuds du maillage
  205. segini icprg,mwrk3
  206. C*NU* segini,icprl
  207.  
  208. nbnds=0
  209. segact mrigid
  210. do 12 ir=1,irigel(/2)
  211. if (irigel(6,ir).ne.2) goto 12
  212. ipt3=irigel(1,ir)
  213. segact ipt3
  214. do ip=1,ipt3.num(/2)
  215. impt=ipt3.num(2,ip)
  216. if (icprg(impt).eq.0) then
  217. nbnds=nbnds+1
  218. icprg(impt)=nbnds
  219. endif
  220. enddo
  221. c* segdes,ipt3
  222. 12 continue
  223. c* segdes,mrigid
  224. * write(ioimp,*) ' nbnds ' , nbnds
  225.  
  226. *C*NU*C --- Tableau inverse de ICPRG
  227. *C*NU* segini idcp
  228. *C*NU* do 20 i=1, nbtpts
  229. *C*NU* if (icprg(i).ne.0) idcp(icprg(i))=i
  230. *C*NU* 20 continue
  231. *C*NU** write(ioimp,*) ' idcp', (idcp(i),i=1,nbnds)
  232. *C*NU*
  233. C --- Tableau des multiplicateurs correspondants aux points des cables
  234. segini idmul
  235. c* segact mrigid
  236. do 22 ir=1,irigel(/2)
  237. if (irigel(6,ir).ne.2) go to 22
  238. ipt3=irigel(1,ir)
  239. c* segact,ipt3
  240. do ip=1,ipt3.num(/2)
  241. idmul(icprg(ipt3.num(2,ip)))=ipt3.num(1,ip)
  242. enddo
  243. segdes,ipt3
  244. 22 continue
  245. c* segdes,mrigid
  246.  
  247. C*NU* Le CHPOINT mchpo2 n'est pas utilise actuellement :
  248. C*NU*C --- Préparation du champ par points de sortie (en norme)
  249. C*NU* nsoupo=1
  250. C*NU* nc=1
  251. C*NU* nat=1
  252. C*NU* segini mchpo2
  253. C*NU* mchpo2.mochde='SEUILS DE FROTTEMENT AUX NOEUDS'
  254. C*NU* mchpo2.mtypoi='FORCES'
  255. C*NU* mchpo2.jattri(1)=2
  256. C*NU* mchpo2.ifopoi=ifour
  257. C*NU** write(ioimp,*)' Segini de MCHPO2',mchpo2
  258. C*NU*c
  259. C*NU* segini msoupo
  260. C*NU* mchpo2.ipchp(nsoupo)=msoupo
  261. C*NU* nocomp(1)='LX'
  262. C*NU* NOHARM(1)=nifour
  263. C*NU** write(ioimp,*)' Segini de MSOUPO,mchpo2.ipchp(nsoupo)',msoupo
  264. C*NU*C --- Création du meleme support du chpoint
  265. C*NU* nbsous=0
  266. C*NU* nbref=0
  267. C*NU* nbnn=1
  268. C*NU* nbelem=nbnds
  269. C*NU* segini meleme
  270. C*NU* do i=1,nbnds
  271. C*NU* num(1,i)=idmul(i)
  272. C*NU* enddo
  273. C*NU** write(ioimp,*)' num', (num(1,i),i=1,nbnds)
  274. C*NU* segdes meleme
  275. C*NU* igeoc= meleme
  276. C*NU* n=nbnds
  277. C*NU* segini mpoval
  278. C*NU* ipoval=mpoval
  279. C*NU* segdes msoupo
  280. C*NU* segdes mchpo2
  281.  
  282. C --- Création d'un modele élémentaire local (a supprimer a la fin)
  283. n1=1
  284. segini mmode3
  285. mn3=0
  286. nmat=0
  287. nfor=0
  288. nobmod=0
  289. segact mmode1
  290. segini,mmode6=mmode1
  291. do isouf=1,icable
  292. imodel=mmode6.kmodel(isouf)
  293. segini,imode6=imodel
  294. mmode6.kmodel(isouf)=imode6
  295. segini mmode7
  296. segini imode7
  297. mmode7.kmodel(1)=imode7
  298. imode7.imamod= imode6.ivamod(1)
  299. imode6.ivamod(1)=mmode7
  300. imode6.tymode(1)='MMODEL'
  301. enddo
  302. segdes mmode1
  303. mmode1=mmode6
  304. mchel5=ipcar
  305. mchel6=ipcon
  306. segact mchel5,mchel6
  307. n35=mchel5.infche(/2)
  308. n36=mchel6.infche(/2)
  309.  
  310. C --- Boucle sur les sous modèles de frottement
  311. C (uniquement les groupes de câbles)
  312. *
  313. * magouille pour ce mettre dans la meme situation que cell que l'on avait
  314. * avec les modeles frottements
  315. *
  316.  
  317. C=======================================
  318. DO 500 isouf = 1, icable
  319. C=======================================
  320. imode1 = mmode1.kmodel(isouf)
  321. C* segact imode1
  322. ipt1 = imode1.imamod
  323. mmode2 = imode1.ivamod(1)
  324.  
  325. C --- On repère dans les sous modèles le modèle groupe de câble
  326. segact mmode2
  327. imode2 = mmode2.kmodel(1)
  328. segact imode2
  329. ipt2 = imode2.imamod
  330.  
  331. * write(ioimp,*)
  332. * write(ioimp,*) '*** Travail sur le sous modele de frottement ***'
  333. * write(ioimp,*) 'Sous modèle numero ',isouf,':',imode1
  334. * write(ioimp,*) 'Sous modèle de type câble de ce modèle',mmode2
  335.  
  336. C --- Création d'un matériau élémentaire
  337. do iou=1,mchel5.imache(/1)
  338. if ( mchel5.imache(iou).eq. ipt1 .or.
  339. $ mchel5.imache(iou).eq. ipt2 ) goto 7755
  340. enddo
  341. call erreur (560)
  342. GOTO 9100
  343. 7755 continue
  344. n1=1
  345. l1=40
  346. n3=n35
  347. segini mchel3
  348. mchel3.imache(1)=ipt1
  349. mchel3.ichaml(1)=mchel5.ichaml(iou)
  350. mchel3.conche(1)=mchel5.conche(iou)
  351. do i=1,n3
  352. mchel3.infche(1,i)=mchel5.infche(iou,i)
  353. enddo
  354.  
  355. C --- Creation des contraintes elementaires
  356. do iou=1,mchel6.imache(/1)
  357. if ( mchel6.imache(iou).eq.ipt1 .or.
  358. $ mchel6.imache(iou).eq.ipt2 ) go to 7756
  359. enddo
  360. call erreur(560)
  361. goto 9100
  362. 7756 continue
  363. n1=1
  364. l1=40
  365. n3=n36
  366. segini mchel4
  367. mchel4.imache(1)=ipt1
  368. mchel4.ichaml(1)=mchel6.ichaml(iou)
  369. mchel4.conche(1)=mchel6.conche(iou)
  370. do i=1,n3
  371. mchel4.infche(1,i)=mchel6.infche(iou,i)
  372. enddo
  373.  
  374. segact,mmode3*mod
  375. mmode3.kmodel(1)=imode2
  376. siezo = 0
  377. call splitag(mmode3,mchel3,mchel4,mmodel,mchel2,
  378. & mchel1,siezo)
  379. segsup,siezo
  380. * call zpmode(mmode3,0)
  381. * call zpmode(mmodel,0)
  382. * write(ioimp,*) ' lister de mchel1 contraintes splitées'
  383. * call zpchel(mchel1,0)
  384. * write(ioimp,*) ' lister de mchel2 materiaux splités'
  385. * call zpchel(mchel2,0)
  386.  
  387. segact mmodel
  388. nszo = kmodel(/1)
  389. segact mchel1,mchel2
  390.  
  391. C --- Boucle sur les sous zones du nouveau modèle, donc sur les cables
  392. C==============================
  393. DO 900 isouc=1,nszo
  394. C==============================
  395. * write(ioimp,*)'*** Travail sur le cable ',isouc,'***'
  396.  
  397. C --- Recherche des noms de composantes des contraintes EFFX
  398. C --- et de coefficients de frottements FF et PHIF
  399. C --- (mchel1=contraintes et mchel2=materiau)
  400. C --- On s'occupe d'abord des contraintes
  401. mcham1=mchel1.ichaml(isouc)
  402. segact,mcham1
  403. C --- Il ne doit y avoir qu'une contrainte de nom EFFX
  404. if (mcham1.nomche(1).ne.'EFFX') then
  405. call erreur(976)
  406. goto 9100
  407. endif
  408. melva1=mcham1.ielval(1)
  409. segdes,mcham1
  410.  
  411. C --- Au tour du materiau
  412. C --- On doit prendre FF puis PHIF
  413. mcham2=mchel2.ichaml(isouc)
  414. segact,mcham2
  415. melva2 = 0
  416. melva3 = 0
  417. do i=1,mcham2.nomche(/2)
  418. * write(ioimp,*)'Nom de la composante',i,'du mchaml materiau du cable',
  419. * & isouc,':',mcham2.nomche(i)
  420. if (mcham2.nomche(i).eq.'FF' ) then
  421. melva2=mcham2.ielval(i)
  422. goto 901
  423. endif
  424. enddo
  425. * write(ioimp,*) 'Composante FF non trouvee dans mchel2'
  426. call erreur(977)
  427. goto 9100
  428. 901 continue
  429. do i=1,mcham2.nomche(/2)
  430. if (mcham2.nomche(i).eq.'PHIF') then
  431. * write(ioimp,*)'Nom de la composante',i,'du mchaml materiau du cable',
  432. * & isouc,':',mcham2.nomche(i)
  433. melva3=mcham2.ielval(i)
  434. go to 902
  435. endif
  436. enddo
  437. * write(ioimp,*) 'Composante PHIF non trouvee dans mchel2'
  438. call erreur(977)
  439. goto 9100
  440. 902 continue
  441. segdes,mcham2
  442.  
  443. * write(ioimp,*) 'Debut des calculs'
  444. C --- Début des calculs a proprement parler
  445. imodel = mmodel.kmodel(isouc)
  446. segact imodel
  447. ipt2 = imodel.imamod
  448. segact ipt2
  449. nbelem = ipt2.num(/2)
  450. nbnn = ipt2.num(/1)
  451.  
  452. C*NU*C --- Renumérotation locale sur les noeuds du maillage
  453. C*NU* do i=1,icprl(/1)
  454. C*NU* icprl(i)=0
  455. C*NU* enddo
  456. C*NU* nbndsl = 0
  457. C*NU* do 25 iel=1,nbelem
  458. C*NU* do 25 i=1,nbnn
  459. C*NU* impt=ipt2.num(i,iel)
  460. C*NU* if (icprl(impt).eq.0) then
  461. C*NU* nbndsl=nbndsl+1
  462. C*NU* icprl(impt)=nbndsl
  463. C*NU** write(ioimp,*)'Ancien noeud : ',impt
  464. C*NU** write(ioimp,*)'Nouveau noeud: ',icprl(impt)
  465. C*NU* endif
  466. C*NU* 25 continue
  467. C*NU*C --- Tableau inverse de ICPRL
  468. C*NU* segini idcpl
  469. C*NU* do 26 i=1,icprl(/1)
  470. C*NU* if (icprl(i).ne.0) idcpl(icprl(i))=i
  471. C*NU* 26 continue
  472.  
  473. C --- 1-ere Boucle sur les éléments pour trouver les grandeurs géométriques
  474. nam = nbelem
  475. segini altrav
  476. C* slon= 0.
  477. C* fai = 0.
  478. iam = 0
  479. DO 5005 iel=1,nbelem
  480. * write(ioimp,*)'Element n°',iel
  481. NC2=ipt2.NUM(1,iel)
  482. NC3=ipt2.NUM(2,iel)
  483. JS2=(NC2-1)*IDIMP1
  484. JS3=(NC3-1)*IDIMP1
  485. IF (iel.EQ.1) THEN
  486. NC1=NC2
  487. JS1=JS2
  488. ELSE
  489. NC1=IPT2.NUM(1,iel-1)
  490. JS1=(NC1-1)*IDIMP1
  491. ENDIF
  492. * PRINT *,'NC1=',NC1,' NC2=',NC2,' NC3=',NC3
  493. C --- Distance entre points
  494. DS1=0.
  495. DS2=0.
  496. DO i=1,IDIM
  497. X2 = XCOOR(JS2+i)
  498. RL(i) = XCOOR(JS3+i) - X2
  499. RS(i) = X2 - XCOOR(JS1+i)
  500. C*NU* X1(i) = XCOOR(JS1+i)
  501. C*NU* X2(i) = XCOOR(JS2+i)
  502. C*NU* X3(i) = XCOOR(JS3+i)
  503. C*NU* RL(i) = X3(i) - X2(i)
  504. C*NU* RS(i) = X2(i) - X1(i)
  505. DS1 = DS1 + RL(i)*RL(i)
  506. DS2 = DS2 + RS(i)*RS(i)
  507. ENDDO
  508. DS1=SQRT(DS1)
  509. IF (DS2.NE.0.) THEN
  510. DS2=SQRT(DS2)
  511. CS1=0.
  512. DO i=1,IDIM
  513. RL(i)=RL(i)/DS1
  514. RS(i)=RS(i)/DS2
  515. CS1=CS1+RL(i)*RS(i)
  516. ENDDO
  517. IF (CS1.GT.1.) CS1=1.
  518. ALFA=ACOS(CS1)
  519. if (abs(alfa).lt.xpetit) alfa=xpetit
  520.  
  521. ELSE
  522. ALFA=0.d0
  523. ENDIF
  524. IAM=IAM+1
  525. C*NU* ANG(IAM)=FAI
  526. C*NU* ACUR(IAM)=SLON
  527. DANG(IAM)=ALFA
  528. DACUR1(IAM)=DS1
  529. DACUR2(IAM)=DS2
  530. C* SLON=SLON+DS1
  531. C* FAI=FAI+ALFA
  532. 5005 CONTINUE
  533. C*NU* ANG(NAM+1) = SLON
  534. C*NU* ACUR(NAM+1) = FAI
  535. * write(ioimp,*)'Longueur du câble: ',slon
  536. * write(ioimp,*)'Déviation du câble: ',fai
  537. * on cherche pour chaque element le Rayon du cercle et donc l'angle
  538. * que l'on met dans dang et on a la longueur dans dacur1
  539. DO 5006 iel=2,NBELEM
  540. dacur2(iel)=(dacur1(iel-1)+dacur1(iel))/(2.d0*dang(iel))
  541. 5006 continue
  542. dacur2(1)=dacur2(2)
  543. do 5007 iel=nbelem-1,2,-1
  544. dacur2(iel)= 0.5 * (dacur2(iel+1)+dacur2(iel))
  545. 5007 continue
  546. do 5008 iel=1,nbelem
  547. dang(iel) =2. * asin( dacur1(iel)/ (2.*dacur2(iel)))
  548. 5008 continue
  549.  
  550. C --- 2-eme Boucle sur les noeuds pour trouver les forces
  551. C --- de frottement, et remplir le champ par points.
  552. segact,melva1,melva2,melva3
  553. icva1 = melva1.velche(/2)
  554. icva2 = melva2.velche(/2)
  555. icva3 = melva3.velche(/2)
  556. jcva1 = min(2,melva1.velche(/1))
  557. jcva2 = min(2,melva2.velche(/1))
  558. jcva3 = min(2,melva3.velche(/1))
  559.  
  560. DO 4005 iel = 1,NBELEM
  561.  
  562. C --- Grandeurs géométriques concernant l'élément
  563. ALFA = DANG(iel)
  564. DST = DACUR1(iel)
  565. C* write(ioimp,*)'DANG n°',iel,DANG(iel)
  566. C* write(ioimp,*)'DACUR1 n°',iel,DACUR1(iel)
  567.  
  568. C --- Calcul de la contrainte moyenne au noeud
  569. i_z = min(iel,icva1)
  570. SIGIC = 0.5 * ( melva1.velche( 1 ,i_z) +
  571. & melva1.velche(jcva1,i_z) )
  572.  
  573. C --- Récuperation des coefficient F1 et F2 moyens
  574. i_z = min(iel,icva2)
  575. F1 = 0.5 * ( melva2.velche( 1 ,i_z) +
  576. & melva2.velche(jcva2,i_z) )
  577. i_z = min(iel,icva3)
  578. F2 = 0.5 * ( melva3.velche( 1 ,i_z) +
  579. & melva3.velche(jcva3,i_z) )
  580. C --- Calcul des forces de frottement
  581. * write(ioimp,*) ' f1 alfa f2 dst ' ,f1,alfa,f2,dst
  582. * write(ioimp,*) ' sigic ' , sigic
  583. SP = F1*ALFA*2.1 + F2*DST
  584. VALFRO = SIGIC * (1.D0-EXP(-SP))
  585. FNORM = 0.5D0 * MAX(VALFRO,XZERO)
  586.  
  587. C --- Numéros globaux des noeuds considérés
  588. ip1 = ipt2.num(1,iel)
  589. ip2 = ipt2.num(2,iel)
  590. lip1 = idmul(icprg(ip1))
  591. lip2 = idmul(icprg(ip2))
  592. jcpr(lip1) = 1
  593. jcpr(lip2) = 1
  594. icablr(lip1) = 1
  595. icablr(lip2) = 1
  596. ycpr(lip1) = FNORM + ycpr(lip1)
  597. ycpr(lip2) = FNORM + ycpr(lip2)
  598. 4005 CONTINUE
  599. segdes,melva1,melva2,melva3
  600. SEGSUP altrav
  601. C*NU* SEGSUP idcpl
  602. C====================
  603. 900 CONTINUE
  604. segdes,mchel1,mchel2
  605. C====================
  606. 500 CONTINUE
  607. C====================
  608. 9100 CONTINUE
  609. segsup,icprg,idmul
  610. segsup,mmode3
  611. if (IERR.ne.0) goto 9900
  612.  
  613. C=======================================================================
  614. C= TRAITEMENT DES MODELES DE TYPE 'COULOMB'
  615. C=======================================================================
  616. 1000 CONTINUE
  617. IF (icoulo.EQ.0) GOTO 2000
  618. mmodel = mocoul
  619. mchelm=ipcar
  620. segact,mchelm
  621. DO 200 m1 = 1, icoulo
  622. segact,mmodel
  623. imodel = kmodel(m1)
  624. SEGACT,imodel
  625. ipt1 = imamod
  626. * recherche rapide d'un maillage correspondant dans le mchaml
  627. DO 210 n2 = 1,imache(/1)
  628. * write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,imache(n2),ipt1
  629. if (imache(n2).eq.ipt1) goto 220
  630. 210 continue
  631. call erreur(472)
  632. goto 9200
  633. 220 continue
  634. mchaml=ichaml(n2)
  635. segact,mchaml
  636. * recherche des bonnes composantes dans le chamelem
  637. * nouveau elem frot PV
  638. imu=0
  639. icohe=0
  640. iadhe=0
  641. melva1=mmodel
  642. melva2=mmodel
  643. melva3=mmodel
  644. do 225 iche=1,nomche(/2)
  645. if (nomche(iche).eq.'MU ') then
  646. imu=iche
  647. melva1=ielval(imu)
  648. segact melva1
  649. icva1 = melva1.velche(/2)
  650. else if (nomche(iche).eq.'COHE') then
  651. icohe=iche
  652. melva2=ielval(icohe)
  653. segact melva2
  654. icva2 = melva2.velche(/2)
  655. else if (nomche(iche).eq.'ADHE') then
  656. iadhe = iche
  657. melva3=ielval(iadhe)
  658. segact melva3
  659. icva3 = melva3.velche(/2)
  660. endif
  661. 225 continue
  662. segdes,mchaml
  663. ipt7=ivamod(1)
  664. SEGACT,ipt7
  665. * attention ipt7 peut être compose
  666. ipt1=ipt7
  667. do 231 is=1,max(1,ipt7.lisous(/1))
  668. if (ipt7.lisous(/1).ne.0) ipt1=ipt7.lisous(is)
  669. segact ipt1
  670. nbele1 = ipt1.num(/2)
  671. nbnoe1 = ipt1.num(/1)
  672. if (nbnoe1.eq.0) goto 231
  673. ** write (6,*) ' excfro ipt1 ',ipt1,nbele1,nbnoe1
  674. do 230 iel=1, nbele1
  675. impt = ipt1.num(1,iel)
  676. * write (6,*) ' excfro impt icpr ',impt,icpr(impt)
  677. if (icpr(impt).eq.0) goto 230
  678. ipasp = 0
  679. xpres = 0.
  680. ypres = 0.
  681. * bon on a un frottement actif. On a le multiplicateur de lagrange
  682. * associe. Il faut le transformer en pression
  683. if (idim.eq.2) then
  684. xpres = xcpr(impt)
  685. * cas de l'adherence
  686. if (iadhe.ne.0) then
  687. if (-xpres+xpetit.gt.
  688. > melva3.velche(1,min(iel,icva3))*0.99999) then
  689. icpr(impt)=0
  690. jcpr(impt)=0
  691. goto 230
  692. endif
  693. if (xpres.lt.xpetit) then
  694. xpres=xpetit
  695. ipasp=1
  696. endif
  697. endif
  698. if (imu.ne.0) then
  699. ypres = ypres + xpres * melva1.velche(1,min(iel,icva1))
  700. endif
  701. if (icohe.ne.0) then
  702. ypres = ypres + melva2.velche(1,min(iel,icva2))
  703. ipasp = 0
  704. endif
  705. if (ipasp.eq.0) then
  706. ip1 = ipt1.num(nbnoe1,iel)
  707. jcpr(ip1) = 1
  708. ycpr(ip1) = ypres
  709. * write (6,*) ' excfro ip1 ypres ',ip1,ypres
  710. endif
  711. * traitement de l'adherence
  712. if (iadhe.ne.0) then
  713. jcpr(impt) = 2
  714. * ycpr(impt) = melva3.velche(1,min(iel,icva3))*xlong
  715. ycpr(impt) = melva3.velche(1,min(iel,icva3))
  716. * write (6,*) ' excfro 2d impt ycpr',impt, ycpr(impt)
  717. endif
  718. * cas 3D
  719. else
  720. * on ne traite que mu pour le moment quelle que soit la formulation
  721. xpres = xcpr(impt)
  722. if (xpres.lt.xpetit) then
  723. xpres = xpetit
  724. ipasp = 1
  725. endif
  726. * cas de l'adherence
  727. if (iadhe.ne.0) then
  728. * write (6,*)'excfro3d,iadhe,impt,xcpr',iadhe,impt
  729. * > ,xcpr(impt)
  730. if (-xpres+xpetit.gt.
  731. > melva3.velche(1,min(iel,icva3))*0.99999) then
  732. icpr(impt)=0
  733. jcpr(impt)=0
  734. goto 230
  735. endif
  736. endif
  737. if (imu.ne.0) then
  738. ypres = xpres * melva1.velche(1,min(iel,icva1))
  739. endif
  740. if (icohe.ne.0) then
  741. ypres = ypres + melva2.velche(1,min(iel,icva2))
  742. ipasp = 0
  743. endif
  744. * ypres est le multiplicateur donnant la force de frottement totale
  745. * on la repartira ultérieurement entre les deux conditions
  746. * elements de contact a nbnoe1 noeuds en 3D (les 2 derniers etant les
  747. * points supports du frottement)
  748. if (ipasp.eq.0) then
  749. ip1 = ipt1.num(nbnoe1-1,iel)
  750. ip2 = ipt1.num(nbnoe1,iel)
  751. jcpr(ip1) = ip2
  752. ycpr(ip1) = ypres
  753. jcpr(ip2) = ip1
  754. ycpr(ip2) = ypres
  755. ** ymcrit1 = max(ymcrit1,abs(ypres))
  756. endif
  757. * traitement de l'adherence
  758. if (iadhe.ne.0) then
  759. jcpr(impt) = -2
  760. * ycpr(impt) = melva3.velche(1,min(iel,icva3))*xlong
  761. ycpr(impt) = melva3.velche(1,min(iel,icva3))
  762. * write (6,*)'excfro3d,impt,ycpr',impt,ycpr(impt)
  763. endif
  764. endif
  765. 230 continue
  766. if (ipt7.lisous(/1).ne.0) segdes ipt1
  767. 231 continue
  768. segdes ipt7,melva1,melva2,melva3
  769. 200 CONTINUE
  770.  
  771. 9200 CONTINUE
  772. segdes mchelm
  773. IF (ierr.ne.0) goto 9999
  774.  
  775. C=======================================================================
  776. C= TRAITEMENT DES MODELES DE TYPE '........'
  777. C=======================================================================
  778. 2000 CONTINUE
  779.  
  780. C=======================================================================
  781. * RECHERCHE DU SENS DU FROTTEMENT (TOUS MODELES)
  782. C=======================================================================
  783. * 1 - On applique la raideur de frottement sur le chpt de deplacement.
  784. * On regarde le signe des multiplicateurs.
  785. * 2 - Si le blocage est actif, il faut maintenir la condition.
  786. * On utilise le signe de la reaction.
  787. *
  788. C*NU* xmcrit = 0.
  789. ymcrit = xpetit
  790. *
  791. mchpo1=mchpin
  792. *
  793. call mucpri(mchpo1,mrigid,mchpoi)
  794. segact mchpoi
  795. * cas du blocage actif etudier les reactions
  796. do 300 isoupo=1,ipchp(/1)
  797. msoupo=ipchp(isoupo)
  798. segact msoupo
  799. mpoval=ipoval
  800. do 310 ic=1,nocomp(/2)
  801. if (nocomp(ic).eq.'FLX ') goto 320
  802. 310 continue
  803. goto 350
  804. 320 continue
  805. ipt2=igeoc
  806. segact ipt2,mpoval
  807. do 330 iel=1,ipt2.num(/2)
  808. C*NU* xmcrit = max(xmcrit,abs(vpocha(iel,ic)))
  809. impt = ipt2.num(1,iel)
  810. if (jcpr(impt).eq.0) goto 330
  811. if (idim.eq.2.or.icablr(impt).eq.1) then
  812. if (jcpr(impt).eq.2) goto 330
  813. jcpr(impt) =-sign(1.1D0,vpocha(iel,ic))
  814. else
  815. zcpr(impt) = vpocha(iel,ic)
  816. endif
  817. 330 continue
  818. segdes,ipt2
  819. 350 continue
  820. segsup,mpoval
  821. segsup,msoupo
  822. 300 continue
  823. segsup mchpoi
  824. C*NU* xmcrit=xmcrit*1.D-10
  825. C*NU** write(ioimp,*) '0 - xmcrit ymcrit dans excfro ',xmcrit
  826.  
  827. * etude des deplacements
  828. segact mchpo1
  829. do 450 isoupo=1,mchpo1.ipchp(/1)
  830. msoupo=mchpo1.ipchp(isoupo)
  831. segact msoupo
  832. do 460 ic=1,nocomp(/2)
  833. if (nocomp(ic).eq.'LX ') goto 470
  834. 460 continue
  835. goto 450
  836. 470 continue
  837. ipt2=igeoc
  838. mpoval=ipoval
  839. segact ipt2,mpoval
  840. do 480 iel=1,ipt2.num(/2)
  841. ymcrit=max(ymcrit,abs(vpocha(iel,ic)))
  842. 480 continue
  843. segdes,ipt2,mpoval
  844. 450 continue
  845. * write(ioimp,*) ' 1 ymcrit1 ',ymcrit1
  846. * write(ioimp,*) ' 2 ymcrit ',ymcrit
  847. ** ymcrit=1.D-6 * max(ymcrit1,ymcrit)
  848. ymcrit=1.D-6 * ymcrit
  849. * write (ioimp,*) ' 3 ymcrit ',ymcrit
  850. C*NU** write(ioimp,*) ' xmcrit ymcrit dans excfro ',xmcrit,ymcrit
  851.  
  852. do 400 isoupo=1,mchpo1.ipchp(/1)
  853. msoupo=mchpo1.ipchp(isoupo)
  854. C* segact msoupo
  855. do 410 ic=1,nocomp(/2)
  856. if (nocomp(ic).eq.'LX ') goto 420
  857. 410 continue
  858. goto 440
  859. 420 continue
  860. ipt2=igeoc
  861. mpoval=ipoval
  862. segact ipt2,mpoval
  863. do 430 iel=1,ipt2.num(/2)
  864. impt = ipt2.num(1,iel)
  865. if (jcpr(impt).eq.0) goto 430
  866. if (idim.eq.2.or.icablr(impt).eq.1) then
  867. if (jcpr(impt).eq.2) goto 430
  868. if (abs(vpocha(iel,ic)).lt.ymcrit) then
  869. jcpr(impt)=0
  870. * write(ioimp,*) ' on met a zero le point', impt
  871. else
  872. jcpr(impt)= -sign(1.1D0,vpocha(iel,ic))
  873. endif
  874. else
  875. if (jcpr(impt).eq.-2) goto 430
  876. kcpr(impt)=1
  877. zcpr(impt)= vpocha(iel,ic)
  878. endif
  879. 430 continue
  880. segdes,ipt2,mpoval
  881. 440 continue
  882. segdes,msoupo
  883. 400 continue
  884. segdes mchpo1
  885. *
  886. * On cree un CHPOINT de MULT de FROTTEMENT qui nous donnera les forces
  887. * de frottement en le multipliant par la raideur
  888. * CHPOINT s'appuyant sur un seul maillage (nsoupo = 1) avec une seule
  889. * composante (nc = 1) de nom 'LX'
  890. jg=jcpr(/1)
  891. nat=1
  892. nsoupo=1
  893. segini mchpoi
  894. jattri(1)=1
  895. mtypoi='LX'
  896. mochde='créé par excfro'
  897. ifopoi=ifomod
  898. nc=1
  899. segini msoupo
  900. ipchp(1)=msoupo
  901. nocomp(1)='LX '
  902. nbnn=1
  903. nbelem=2*jg
  904. nbref=0
  905. nbsous=0
  906. segini ipt7
  907. ipt7.itypel=1
  908. igeoc=ipt7
  909. n=nbelem
  910. segini mpoval
  911. ipoval=mpoval
  912. ip=0
  913. do 600 ipr=1,jg
  914. if (jcpr(ipr).eq.0) goto 600
  915. if (idim.eq.2.and.jcpr(ipr).eq.2) jcpr(ipr)=1
  916. ip=ip+1
  917. ipt7.num(1,ip)=ipr
  918. if (idim.eq.2.or.icablr(ipr).eq.1) then
  919. vpocha(ip,1)= ycpr(ipr)*jcpr(ipr)
  920. * write(6,*)'excfro 2D ip ipr vpocha(ip,1)'
  921. * > ,ip,ipr,vpocha(ip,1)
  922. * Si la force est trop petite la direction peut etre fausse
  923. if (abs(vpocha(ip,1)).lt.ymcrit) then
  924. ip=ip-1
  925. * write(ioimp,*) ' non prise en compte du point i',ipr
  926. endif
  927. else
  928. if (jcpr(ipr).eq.-2) then
  929. vpocha(ip,1)= ycpr(ipr)
  930. * write(6,*)'excfro 3DD ip ipr vpocha(ip,1)'
  931. * > ,ip,ipr,vpocha(ip,1)
  932. * Si la force est trop petite la direction peut etre fausse
  933. if (abs(vpocha(ip,1)).lt.ymcrit) then
  934. ip=ip-1
  935. * write(ioimp,*) ' non prise en compte du point i',ipr
  936. endif
  937. else
  938. ss=ycpr(ipr)
  939. if (abs(ss).lt.ymcrit) then
  940. * write(ioimp,*) ' pt elimine ',ipr,abs(ss),ymcrit
  941. ip=ip-1
  942. goto 600
  943. endif
  944. * en 3d il faut corriger en fonction de la direction du depl ou des reac
  945. idu=jcpr(ipr)
  946. xx=zcpr(ipr)
  947. yy=zcpr(idu)
  948. * write(ioimp,*) ' xx yy ',xx,yy
  949. * condition sur une direction, pas sur l'autre
  950. if (kcpr(ipr).eq.0.and.kcpr(idu).ne.0) then
  951. * write(ioimp,*) 'prob ',ipr,kcpr(ipr),idu,kcpr(idu)
  952. yy=sign(SQRT(MAX(ss*ss-xx*xx,xzero)),yy)
  953. endif
  954. if (kcpr(ipr).ne.0.and.kcpr(idu).eq.0) then
  955. * write(ioimp,*) 'prob ',ipr,kcpr(ipr),idu,kcpr(idu)
  956. xx=sign(sqrt(max(ss*ss-yy*yy,xzero)),xx)
  957. endif
  958. ss=sqrt(max(xx*xx+yy*yy,xzero))
  959. if (ss.lt.xpetit) ss=1.
  960. * write(ioimp,*) ' xx yy ss ',xx,yy,ss
  961. if (ss.lt.ymcrit) then
  962. if (kcpr(ipr).eq.1.or.kcpr(idu).eq.1) then
  963. * write(ioimp,*) ' 2-pt elimine ',ipr,ss,abs(ycpr(ipr))
  964. ip=ip-1
  965. goto 600
  966. endif
  967. endif
  968. xx=xx/ss
  969. * yy=yy/ss
  970. * write(ioimp,*) 'prim dual ',ipr,idu,xx,yy
  971. vpocha(ip,1)= -ycpr(ipr)*xx
  972. * write(ioimp,*) ' force ',vpocha(ip,1),ipr,idu
  973. if (abs(vpocha(ip,1)).lt.ymcrit) then
  974. * write(ioimp,*) ' condition eliminee ',ipr,kcpr(ipr),kcpr(idu),ss
  975. ip=ip-1
  976. endif
  977. endif
  978. endif
  979. 600 continue
  980. nbelem=ip
  981. n=ip
  982. segadj ipt7,mpoval
  983. segdes,mpoval,ipt7,msoupo,mchpoi
  984.  
  985. * En cas de presence d'un champ, mettre au minimum les valeurs d'icelui
  986. if (mchpi2.eq.0) goto 9000
  987. * on remet a zero icpr et xcpr
  988. do impt=1,icpr(/1)
  989. icpr(impt) = 0
  990. xcpr(impt) = 0.
  991. enddo
  992. mchpo1=mchpi2
  993. segact mchpo1
  994. do 700 isoupo=1,mchpo1.ipchp(/1)
  995. msoupo=mchpo1.ipchp(isoupo)
  996. segact msoupo
  997. do 710 ic=1,nocomp(/2)
  998. if (nocomp(ic).eq.'LX ') goto 720
  999. 710 continue
  1000. goto 740
  1001. 720 continue
  1002. ipt2=igeoc
  1003. mpoval=ipoval
  1004. segact ipt2,mpoval
  1005. do 730 iel=1,ipt2.num(/2)
  1006. impt = ipt2.num(1,iel)
  1007. r_z = vpocha(iel,ic)
  1008. if (r_z.lt.-ymcrit) then
  1009. icpr(impt) = -1
  1010. else if (r_z.gt.ymcrit) then
  1011. icpr(impt) = 1
  1012. endif
  1013. xcpr(impt) = r_z
  1014. 730 continue
  1015. segdes ipt2,mpoval
  1016. 740 continue
  1017. segdes msoupo
  1018. 700 continue
  1019. segdes mchpo1
  1020. * Mise a jour du CHPOINT de MULT de FROTTEMENT (CHPOINT DE 'LX')
  1021. * Rappel : nsoupo = ipchp(/1) = 1 et nc = nocomp(/2) = 1 'LX'
  1022. segact mchpoi
  1023. do 800 isoupo=1,ipchp(/1)
  1024. msoupo=ipchp(isoupo)
  1025. segact msoupo
  1026. do 810 ic=1,nocomp(/2)
  1027. if (nocomp(ic).eq.'LX ') goto 820
  1028. 810 continue
  1029. goto 850
  1030. 820 continue
  1031. ipt2=igeoc
  1032. mpoval=ipoval
  1033. segact ipt2,mpoval*mod
  1034. nbelem=ipt2.num(/2)
  1035. do 830 iel=1,nbelem
  1036. impt=ipt2.num(1,iel)
  1037. if (icpr(impt).eq.-1) then
  1038. * write(ioimp,*) ' excfro ancien ',impt,vpocha(iel,ic),xcpr(impt)
  1039. vpocha(iel,ic)=min(vpocha(iel,ic),xcpr(impt))
  1040. else if (icpr(impt).eq. 1) then
  1041. * write(ioimp,*) ' excfro ancien ',impt,vpocha(iel,ic),xcpr(impt)
  1042. vpocha(iel,ic)=max(vpocha(iel,ic),xcpr(impt))
  1043. endif
  1044. icpr(impt)=0
  1045. 830 continue
  1046. * on rajoute les autres valeurs du champ precedent
  1047. nbeleo=nbelem
  1048. do 835 impt=1,icpr(/1)
  1049. if (icpr(impt).ne.0) nbelem=nbelem+1
  1050. 835 continue
  1051. if (nbelem.ne.nbeleo) then
  1052. nbnn=1
  1053. nbsous=0
  1054. nbref=0
  1055. segadj ipt2
  1056. n=nbelem
  1057. nc=1
  1058. segadj mpoval
  1059. do 840 impt=1,icpr(/1)
  1060. if (icpr(impt).eq.0) goto 840
  1061. nbeleo=nbeleo+1
  1062. ipt2.num(1,nbeleo)=impt
  1063. vpocha(nbeleo,ic)=xcpr(impt)
  1064. * write(ioimp,*) ' excfro nouveau',impt,xcpr(impt)
  1065. 840 continue
  1066. endif
  1067. segdes ipt2,mpoval
  1068. 850 continue
  1069. segdes msoupo
  1070. 800 continue
  1071. segdes,mchpoi
  1072.  
  1073. 9000 CONTINUE
  1074. call ecrobj('CHPOINT',mchpoi)
  1075.  
  1076. 9900 CONTINUE
  1077. segsup icpr,jcpr,xcpr,ycpr
  1078. if (idim.eq.3) segsup zcpr,kcpr
  1079. SEGSUP,icablr
  1080. if (mwrk3.ne.0) SEGSUP,mwrk3
  1081.  
  1082. 9999 CONTINUE
  1083. mmodel = modini
  1084. DO isoum = 1, nsoum
  1085. imodel = kmodel(isoum)
  1086. SEGDES,imodel
  1087. ENDDO
  1088. SEGDES,mmodel
  1089.  
  1090. mmode1 = mocabl
  1091. mmode2 = mocoul
  1092. SEGSUP,mmode1,mmode2
  1093.  
  1094. RETURN
  1095. END
  1096.  
  1097.  
  1098.  
  1099.  
  1100.  
  1101.  
  1102.  
  1103.  
  1104.  

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