Télécharger la plaquette Cast3M


fr en

Télécharger exccab.eso

Retour à la liste

Numérotation des lignes : 

exccab

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

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