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

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