Télécharger excfro.eso

Retour à la liste

Numérotation des lignes :

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

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