Télécharger la plaquette Cast3M


fr en

Télécharger excfro.eso

Retour à la liste

Numérotation des lignes : 

excfro

  1. C EXCFRO    SOURCE    CB215821  24/04/12    21:15:50     11897          
  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. *  tableau d'indirection generale
  32.       SEGMENT ncpr(nbpts)
  33.       SEGMENT icpr(nbptl)
  34.       SEGMENT xcpr(nbptl)
  35.       SEGMENT jcpr(nbptl)
  36.       SEGMENT ycpr(nbptl)
  37. *  segments 3D
  38.       SEGMENT kcpr(nbptl)
  39.       SEGMENT zcprd(nbptl)
  40.       SEGMENT zcprf(nbptl)
  41. *  tableau donnant la sous matrice associee a un multiplicateur
  42.       segment irilx(nbptl)
  43. *  tableau donnant le descripteur associe a un multiplicateur
  44.       segment irides(nbptl)
  45. *  tableau donnant l'element dans la sous raideur
  46.       segment iellx(nbptl)
  47. *  tableau donnant le maillage du modele associe a un multiplicateur
  48.       segment imolx(nbptl)
  49. *  tableau donnant l'element dans le maillage
  50.       segment imllx(nbptl)
  51. *  tableau donnant le coefficient  a partir du LX
  52.       segment xdcp(nbptl)
  53.  
  54. *
  55. C
  56. *
  57. *  pour tester le 3D non symetrique, mettre NS3D a vrai
  58. *
  59.  
  60.  
  61.  
  62. C
  63. *
  64.       logical ns3d,ns2d
  65. *
  66.       ns3d=.false.
  67.       ns2d=.true.
  68. *
  69. C=======================================================================
  70. C= LECTURE DES DONNEES DE L'OPERATEUR
  71. C=======================================================================
  72. **    write(6,*) ' entree dans excfro '
  73.       call lirobj('MMODEL  ',modini,1,iretou)
  74.       call actobj('MMODEL  ',modini,1)
  75.       if (ierr.ne.0) return
  76.       call lirobj('MCHAML  ',icham1,1,iretou)
  77.       call actobj('MCHAML  ',icham1,1)
  78.       if (ierr.ne.0) return
  79.       call lirobj('CHPOINT ',mchpin,1,iretou)
  80.       call actobj('CHPOINT ',mchpin,1)
  81.       if (ierr.ne.0) return
  82.       call lirobj('RIGIDITE',mrigid,1,iretou)
  83.       if (ierr.ne.0) return
  84.       icham2 = 0
  85.       mchpi2 = 0
  86.       call lirent(noncon,1,iretou)
  87.       if (ierr.ne.0) return
  88. C- Tri des MCHAML donnes en entree
  89.       CALL rngcha(icham1,icham2,'CARACTERISTIQUES','CONTRAINTES',
  90.      &            ipcar,ipcon)
  91.       if (ierr.ne.0) return
  92. C- Le champ de caracteristiques est toujours obligatoire !
  93.       if (ipcar.eq.0) then
  94.         moterr(1:16)='CARACTERISTIQUES'
  95.         call erreur(565)
  96.         return
  97.       endif
  98.  
  99. C=======================================================================
  100. C= QUELQUES INITIALISATIONS
  101. C=======================================================================
  102.       idimp1 = IDIM + 1
  103. *      write(ioimp,*) ' Entree dans excfro jcpr(ncpr(/1) ',nbpts
  104.  
  105.  
  106. C=======================================================================
  107. C= ANALYSE DU MODELE DU FROTTEMENT et DECOUPAGE EN MODELES DE MEME TYPE
  108. C=======================================================================
  109.       mmodel=modini
  110.       SEGACT,mmodel
  111.       nsoum = kmodel(/1)
  112.  
  113.       n1 = nsoum
  114.       icoulo = 0
  115.       SEGINI,mmode2
  116.       mocoul = mmode2
  117.  
  118.       DO isoum = 1, nsoum
  119.         imodel = kmodel(isoum)
  120.         SEGACT,imodel
  121. *        write(ioimp,*) isoum,' formod(1) matmod(1)',formod(1),matmod(1)
  122.         IF (formod(1).ne.'CONTACT') then
  123.            moterr(1:16)='CONTACT'
  124.            call erreur(719)
  125.            goto 9999
  126.         ENDIF
  127.         icoul=0
  128.         do inma=1,matmod(/2)
  129.           if ( matmod(inma).eq.'COULOMB') icoul=1
  130.         enddo
  131.         IF (icoul.eq.1) THEN
  132. *          IF (idim.eq.2.and.NEFMOD.NE.107) call erreur(16)
  133. *          IF (idim.eq.3.and.NEFMOD.NE.165) call erreur(16)
  134.           icoulo = icoulo + 1
  135.           mmode2.kmodel(icoulo) = imodel
  136.         ELSE
  137. **        write(ioimp,*) 'Modele de FROTTEMENT non implemente'
  138. **        call erreur(21)
  139.         ENDIF
  140.         IF (ierr.ne.0) goto 9999
  141.       ENDDO
  142.       if(n1.ne.icoulo) then
  143.          n1=icoulo
  144.          segadj mmode2
  145.       endif
  146.  
  147.  
  148. C=======================================================================
  149. *  noter les contacts actifs dans le champ de deplacement
  150. C=======================================================================
  151. *  ncpr tableau des multiplicateurs de lagrange
  152. *  recuperes dans le model
  153.       segini ncpr
  154.       nbptl=0
  155.       mmodel=mocoul
  156.       segact mmodel
  157.       do m1=1,kmodel(/1)
  158.         imodel=kmodel(m1)
  159.         segact imodel
  160.         ipt7=imamod
  161.         ipt1=ipt7
  162.         segact ipt7
  163.         do is=1,max(1,ipt7.lisous(/1))
  164.          if(ipt7.lisous(/1).ne.0) ipt1=ipt7.lisous(is)
  165.          segact ipt1
  166.          nbele1=ipt1.num(/2)
  167.          nbnoe1=ipt1.num(/1)
  168.          do iel=1,nbele1
  169.           impt=ipt1.num(1,iel)
  170.           if(ncpr(impt).eq.0) then
  171.            nbptl=nbptl+1
  172.            ncpr(impt)=nbptl
  173.           endif
  174. *  cas des cables mult en 2
  175. **        impt=ipt1.num(2,iel)
  176. **        if(ncpr(impt).eq.0) then
  177. **         nbptl=nbptl+1
  178. **         ncpr(impt)=nbptl
  179. **        endif
  180.           impt=ipt1.num(nbnoe1,iel)
  181.           if(ncpr(impt).eq.0) then
  182.            nbptl=nbptl+1
  183.            ncpr(impt)=nbptl
  184.           endif
  185.           if(idim.eq.3) then
  186.           impt=ipt1.num(nbnoe1-1,iel)
  187.           if(ncpr(impt).eq.0) then
  188.            nbptl=nbptl+1
  189.            ncpr(impt)=nbptl
  190.           endif
  191.           endif
  192.          enddo
  193.         enddo
  194.       enddo
  195.  
  196. ** maintenant on cree les tableaux permettant de retrouver la raideur et le modele a partir du lx
  197.       segini irilx,iellx,irides
  198.       segact mrigid
  199.       do 11 ir=1,irigel(/2)
  200.         if (irigel(6,ir).eq.0) goto 11
  201.         descr=irigel(3,ir)
  202.         segact descr
  203.         if (lisinc(1).ne.'LX  ') goto 11
  204.         meleme=irigel(1,ir)
  205.         segact meleme
  206.         do iel=1,num(/2)
  207.          impt=ncpr(num(1,iel))
  208.          if(impt.ne.0) then
  209.            if (irilx(impt).ne.0) call erreur(5)
  210.            irilx(impt)=irigel(4,ir)
  211.            irides(impt)=irigel(3,ir)
  212.            iellx(impt)=iel
  213.          endif
  214.         enddo
  215.   11  continue
  216.       segini imolx,imllx
  217.       do m1=1,kmodel(/1)
  218.         imodel=kmodel(m1)
  219.         ipt7=imamod
  220.         meleme=ipt7
  221.         segact meleme
  222.         do is=1,max(1,ipt7.lisous(/1))
  223.          if(ipt7.lisous(/1).ne.0) meleme=ipt7.lisous(is)
  224.          segact meleme
  225.          do iel=1,num(/2)
  226.            impt=ncpr(num(1,iel))
  227.            if (imolx(impt).ne.0) call erreur(5)
  228.            imolx(impt)=meleme
  229.            imllx(impt)=iel
  230.          enddo
  231.         enddo
  232.       enddo
  233. *
  234.  
  235.  
  236.  
  237. **    write(6,*) ' nbptl nbpts dans excfro ',nbptl,nbpts
  238.       segini icpr,xcpr,jcpr,ycpr,xdcp
  239.                      segini zcprd,zcprf,kcpr
  240.  
  241. *      call ecchpo(mchpin,0)
  242.       mchpo1=mchpin
  243.       segact mchpo1
  244.       xcprmf=xpetit/xzprec/xzprec
  245.       xcprmd=xpetit/xzprec/xzprec
  246.       do 100 isoupo=1,mchpo1.ipchp(/1)
  247.         msoupo=mchpo1.ipchp(isoupo)
  248.         segact msoupo
  249.         do 110 ic=1,nocomp(/2)
  250.           if (nocomp(ic).ne.'LX  ') go to 110
  251.           ipt2=igeoc
  252.           mpoval=ipoval
  253.           segact ipt2,mpoval
  254.           nbel2 = ipt2.num(/2)
  255.           do 130 iel=1,nbel2
  256.             imptr = ncpr(ipt2.num(1,iel))
  257.             if (imptr.ne.0) then
  258.              icpr(imptr)=1
  259.              xcpr(imptr)=xcpr(imptr)+vpocha(iel,ic)
  260.              xcprmf=max(xcprmf,abs(xcpr(imptr)))
  261. *            write(6,*) 'impt imptr icpr xcpr ',impt,imptr,xcpr(imptr)
  262.             endif
  263.  130      continue
  264.  110    continue
  265.  100  continue
  266.  
  267.  
  268. C=======================================================================
  269. C= TRAITEMENT DES MODELES DE TYPE 'COULOMB'
  270. C=======================================================================
  271.  1000 CONTINUE
  272.       IF (icoulo.EQ.0) GOTO 2000
  273.       mmodel = mocoul
  274.       mchelm=ipcar
  275.       segact,mchelm
  276.       segact,mmodel
  277.       DO 200 m1 = 1, icoulo
  278.         imodel = kmodel(m1)
  279.         SEGACT,imodel
  280.         ipt1 = imamod
  281.         segact ipt1
  282. *  recherche rapide d'un maillage correspondant dans le mchaml
  283.         DO 210 n2 = 1,imache(/1)
  284. *          write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,imache(n2),ipt1
  285.           if (imache(n2).eq.ipt1) goto 220
  286.  210    continue
  287.         call erreur(472)
  288.         goto 9200
  289.  220    continue
  290.         mchaml=ichaml(n2)
  291.         ipt7=imache(n2)
  292.         segact,mchaml,ipt7
  293. *  recherche des bonnes composantes dans le chamelem
  294. *  nouveau elem frot PV
  295.         imu=0
  296.         icohe=0
  297.         melva1=mmodel
  298.         melva2=mmodel
  299.         melva3=mmodel
  300.         do 225 iche=1,nomche(/2)
  301.           if (nomche(iche).eq.'MU  ') then
  302.             imu=iche
  303.             melva1=ielval(imu)
  304.             segact melva1
  305.             icva1 = melva1.velche(/2)
  306.             do i=1,ipt7.num(/2)
  307. *            write(6,*) ' creation xdcp ',ipt7.num(1,i),
  308. *    >       ncpr(ipt7.num(1,i)),
  309. *    >          melva1.velche(1,min(i,icva1))
  310.              xdcp(ncpr(ipt7.num(1,i)))=melva1.velche(1,min(i,icva1))
  311.             enddo
  312.           else if (nomche(iche).eq.'COHE') then
  313.             icohe=iche
  314.             melva2=ielval(icohe)
  315.             segact melva2
  316.             icva2 = melva2.velche(/2)
  317.           endif
  318.  225    continue
  319.         ipt7=imamod
  320.         SEGACT,ipt7
  321. *  attention ipt7 peut etre compose
  322.         ipt1=ipt7
  323.         do 231 is=1,max(1,ipt7.lisous(/1))
  324.         if (ipt7.lisous(/1).ne.0) ipt1=ipt7.lisous(is)
  325.         segact ipt1
  326.         nbele1 = ipt1.num(/2)
  327.         nbnoe1 = ipt1.num(/1)
  328.         if (nbnoe1.eq.0) goto 231
  329. *  on met la raideur additionnelle de frottement meme si contact inactif
  330.         do 230 iel=1, nbele1
  331.           impt = ipt1.num(1,iel)
  332.           imptr=ncpr(impt)
  333. **        if (icpr(imptr).eq.0) goto 230
  334.           ipasp = 0
  335.           xpres= xcprmf*xzprec
  336.           ypres = 0.
  337. *  bon on a un frottement actif ou non. On a le multiplicateur de lagrange
  338. *  associe. Il faut le transformer en pression
  339.           if (icpr(imptr).ne.0) then
  340.               xpres = xcpr(imptr)
  341.           endif
  342.           if (idim.eq.2) then
  343.             if (imu.ne.0) then
  344.               ypres = ypres + xpres * melva1.velche(1,min(iel,icva1))
  345.             endif
  346.             if (icohe.ne.0.and.icpr(imptr).ne.0) then
  347. *  cohesion uniquement si contact
  348.               ypres = ypres + melva2.velche(1,min(iel,icva2))
  349.             endif
  350.             ip1 = ipt1.num(nbnoe1,iel)
  351.             jcpr(ncpr(ip1)) = ip1
  352.             ycpr(ncpr(ip1)) = ypres
  353. *  cas 3D
  354.           else
  355. *  on ne traite que mu pour le moment quelle que soit la formulation
  356.             if (imu.ne.0) then
  357.               ypres = xpres * melva1.velche(1,min(iel,icva1))
  358.             endif
  359.             if (icohe.ne.0.and.icpr(imptr).ne.0) then
  360.               ypres = ypres + melva2.velche(1,min(iel,icva2))
  361.             endif
  362. * ypres est le multiplicateur donnant la force de frottement totale
  363. * on la repartira ulterieurement entre les deux conditions
  364. * elements de contact a nbnoe1 noeuds en 3D (les 2 derniers etant les
  365. * points supports du frottement)
  366. * on met la moitié ici car
  367.             ip1 = ipt1.num(nbnoe1-1,iel)
  368.             ip2 = ipt1.num(nbnoe1,iel)
  369.             jcpr(ncpr(ip1)) = ip2
  370.             ycpr(ncpr(ip1)) = ypres
  371.             jcpr(ncpr(ip2)) = ip1
  372.             ycpr(ncpr(ip2)) = ypres
  373.           endif
  374.  230    continue
  375.  231    continue
  376.  200  CONTINUE
  377.  9200 CONTINUE
  378.       IF (ierr.ne.0) goto 9999
  379.  
  380. C=======================================================================
  381. C= TRAITEMENT DES MODELES DE TYPE '........'
  382. C=======================================================================
  383.  2000 CONTINUE
  384.  
  385. C=======================================================================
  386. *  RECHERCHE DU SENS DU FROTTEMENT (TOUS MODELES)
  387. C=======================================================================
  388. * 1 - On applique la raideur de frottement sur le chpt de deplacement.
  389. *     On regarde le signe des multiplicateurs.
  390. * 2 - Si le blocage est actif, il faut maintenir la condition.
  391. *     On utilise le signe de la reaction.
  392. *
  393.       ymcritd = xpetit/xzprec
  394.       ymcritf = xpetit/xzprec
  395.       ymcritf=max(xpetit,xcprmf)
  396. *  on remplit dans tous les cas les reactions
  397.       segact mcoord
  398. *  si le blocage est actif on s'en servira
  399. **    call ecchpo(mchpo1,0)
  400.       segact mchpo1
  401.       do 450 isoupo=1,mchpo1.ipchp(/1)
  402.         msoupo=mchpo1.ipchp(isoupo)
  403.         segact msoupo
  404.         mpoval=ipoval
  405.         segact mpoval
  406.         ipt2=igeoc
  407.         segact ipt2
  408.           do  il=1,ipt2.num(/2)
  409.           xg=-xgrand
  410.           do  ip=2,ipt2.num(/1)
  411.             it1=ipt2.num(ip-1,il)
  412.             it2=ipt2.num(ip,il)
  413.             xg=max(xg,abs(xcoor((it1-1)*(idim+1)+1)-
  414.      >            xcoor((it2-1)*(idim+1)+1)))
  415.             xg=max(xg,abs(xcoor((it1-1)*(idim+1)+2)-
  416.      >            xcoor((it2-1)*(idim+1)+2)))
  417.             if(idim.eq.3)
  418.      >      xg=max(xg,abs(xcoor((it1-1)*(idim+1)+3)-
  419.      >            xcoor((it2-1)*(idim+1)+3)))
  420.           enddo
  421.             ymcritd=max(ymcritd,xg)
  422.           enddo
  423.         do 465 ic=1,nocomp(/2)
  424.         if (nocomp(ic).eq.'LX  ') then
  425.           do  iel=1,vpocha(/1)
  426.             ymcritf=max(ymcritf,abs(vpocha(iel,ic)))
  427.           enddo
  428.         else
  429.           do  iel=1,vpocha(/1)
  430.             ymcritd=max(ymcritd,abs(vpocha(iel,ic)))
  431.           enddo
  432.         endif
  433.  465    continue
  434.  450  continue
  435.       ymcritf=1d-10 * ymcritf
  436. **    write (ioimp,*) ' 3 ymcritf ',ymcritf
  437.  
  438.       do 400 isoupo=1,mchpo1.ipchp(/1)
  439.         msoupo=mchpo1.ipchp(isoupo)
  440.         mpoval=ipoval
  441.         do 405 ic=1,nocomp(/2)
  442.         if (nocomp(ic).ne.'LX  ') goto 405
  443.         ipt2=igeoc
  444.         do 430 iel=1,ipt2.num(/2)
  445.           impt = ipt2.num(1,iel)
  446.           imptr=ncpr(impt)
  447.           if (imptr.eq.0) goto 430
  448. *  pas de contact => pas de reaction on passe
  449.           if (jcpr(imptr).eq.0) goto 430
  450.           if (abs(vpocha(iel,ic)).lt.ymcritf)  then
  451.           else
  452.                kcpr(imptr)=0
  453.                zcprf(imptr)= vpocha(iel,ic)
  454.           endif
  455.  430    continue
  456.  405    continue
  457.  400  continue
  458.  
  459. *
  460. *
  461. *  etude des deplacements
  462. *  si le deplacement est non nul, c'est que le frottement glisse
  463. *  et on met alors la direction du glissement
  464.       mchpo1=mchpin
  465.       call mucpri(mchpo1,mrigid,mchpoi)
  466.       segact mchpoi
  467.       do 290 isoupo=1,ipchp(/1)
  468.         msoupo=ipchp(isoupo)
  469.         segact msoupo
  470.         mpoval=ipoval
  471.         segact mpoval
  472.         do ic=1,nocomp(/2)
  473.         if (nocomp(ic).ne.'FLX ') goto 290
  474.         ipt2=igeoc
  475.         segact ipt2
  476.         do 293 iel=1,vpocha(/1)
  477.           ymcritd=max(ymcritd,abs(vpocha(iel,ic)))
  478.  293    continue
  479.         enddo
  480.  290  continue
  481.       ymcritd=ymcritd*1d-10
  482.  
  483.       do 300 isoupo=1,ipchp(/1)
  484.         msoupo=ipchp(isoupo)
  485.         mpoval=ipoval
  486.         do 305 ic=1,nocomp(/2)
  487.         if (nocomp(ic).ne.'FLX ') goto 300
  488.         ipt2=igeoc
  489.         do 330 iel=1,ipt2.num(/2)
  490.           impt = ipt2.num(1,iel)
  491.           imptr=ncpr(impt)
  492.           if (imptr.eq.0) goto 331
  493. ** si jcpr nul, pas de contact donc pas de force de frottement
  494.           if (jcpr(imptr).eq.0) goto 331
  495.           if (abs(vpocha(iel,ic)).gt.ymcritd) then
  496.             kcpr(imptr)=1
  497.             zcprd(imptr) = vpocha(iel,ic)
  498.           else
  499.           endif
  500.           goto 330
  501.  331    continue
  502. **           write(6,*) 'point saute ',impt,imptr
  503.  330    continue
  504.  305    continue
  505.         segsup,mpoval
  506.         segsup,msoupo
  507.  300  continue
  508.       segsup mchpoi
  509.  
  510. *
  511. * On cree un CHPOINT de MULT de FROTTEMENT qui nous donnera les forces
  512. * de frottement en le multipliant par la raideur
  513. * CHPOINT s'appuyant sur un seul maillage (nsoupo = 1) avec une seule
  514. * composante (nc = 1) de nom 'LX'
  515.       jg=jcpr(/1)
  516.       nat=1
  517.       nsoupo=1
  518.       segini mchpoi
  519.       jattri(1)=2
  520.       mtypoi='LX'
  521.       mochde='cree par excfro'
  522.       ifopoi=ifour
  523.       nc=1
  524.       segini msoupo
  525.       ipchp(1)=msoupo
  526.       nocomp(1)='LX  '
  527.       nbnn=1
  528.       nbelem=2*jg
  529.       nbref=0
  530.       nbsous=0
  531.       segini ipt7
  532.       ipt7.itypel=1
  533.       igeoc=ipt7
  534.       n=nbelem
  535.       segini mpoval
  536.       ipoval=mpoval
  537.       ip=0
  538. ** boucle sur le modele pour recuperer les conditions
  539.       mmodel=mocoul
  540.       DO 610 m1 = 1, icoulo
  541.         imodel = kmodel(m1)
  542.         ipt5=imamod
  543.         SEGACT,ipt5
  544. *  attention ipt5 peut etre compose
  545.         ipt1=ipt5
  546.         do 631 is=1,max(1,ipt5.lisous(/1))
  547.         if (ipt5.lisous(/1).ne.0) ipt1=ipt5.lisous(is)
  548.         segact ipt1
  549.         if(ipt1.itypel.ne.22) call erreur(5)
  550.         nbnn1=ipt1.num(/1)
  551.       do 635 iel=1,ipt1.num(/2)
  552.         ipr=ipt1.num(nbnn1,iel)
  553.         lxf1=ncpr(ipr)
  554. *  test si on a bien contact
  555.         if (idim.eq.2) then
  556.           if(lxf1.eq.0) goto 635
  557.           xxd=zcprd(lxf1)
  558.           xxf=zcprf(lxf1)
  559. *  deplacements non nuls, on les prend
  560.           ssd=abs(xxd)
  561.           ssf=abs(xxf)
  562.           if (ssd.lt.xpetit) ssd=1.
  563.           if (ssf.lt.xpetit) ssf=1.
  564.           if (abs(xxd).gt.ymcritd) then
  565.              xx= xxd/ssd
  566. **           write(6,*) ' point glisse ',ipr,xx
  567. *  sinon les forces
  568.           elseif(abs(xxf).gt.ymcritf) then
  569.                 xx= xxf/ssf
  570. **           write(6,*) ' point bloque ',ipr,xx
  571.           else
  572. *  pas de deplacements, pas de forces
  573.               jcpr(lxf1)=0
  574.               xdcp(lxf1)=0.d0
  575.               xx=0.d0
  576.               goto 635
  577.           endif
  578.  601      continue
  579.           ip=ip+1
  580.           ipt7.num(1,ip)=ipr
  581. *  en 2d utilisation de la formulation non symetrique
  582. *  on laisse juste une force residuelle pour passer le signe a excite
  583.           vpocha(ip,1)=-ycpr(lxf1)*xx
  584.           if(ns2d) vpocha(ip,1)=vpocha(ip,1)*xzprec
  585.           xdcp(lxf1)= xx
  586.         else
  587. * en 3d il faut corriger en fonction de la direction du depl et des reac
  588.           idu=ipt1.num(nbnn1-1,iel)
  589.           lxf2=ncpr(idu)
  590.           if(lxf1.eq.0.and.lxf2.eq.0) goto 635
  591.           xxd=0.d0
  592.           xxf=0.d0
  593.           yyd=0.d0
  594.           yyf=0.d0
  595.           if (lxf1.ne.0) then 
  596.            xxd=zcprd(lxf1)
  597.            xxf=zcprf(lxf1)
  598.           endif
  599.           if (lxf2.ne.0) then
  600.            yyd=zcprd(lxf2)
  601.            yyf=zcprf(lxf2)
  602.           endif
  603. *  deplacements non nuls, on les prend
  604.              ssd=sqrt(xxd**2+yyd**2)
  605.              ssf=sqrt(xxf**2+yyf**2)
  606.              if (ssd.lt.ymcritd) then
  607. **pv          write(6,*) 'point bloque ',ipr,idu
  608.               ssd=xgrand
  609.              endif
  610.              if (ssf.lt.ymcritf) then
  611. **pv           write(6,*) 'point bloque ',ipr,idu
  612.                ssf=xgrand
  613.              endif
  614. *  on moyenne la direction entre force et deplacement
  615. *  sauf en non convergence ou on garde la direction forces
  616. *  noncon est le compteur d'iteration
  617. *  les 5 premieres, on ne relaxe pas beaucoup pour obtenir la bonne direction de glissement
  618. *  ensuite on relaxe de plus en plus pour converger sur la direction (meme fausse)
  619.              if(ssd.lt.xsgran) then
  620.               if(noncon.gt.6) ssd=ssd*2*(noncon-6)
  621.               if(noncon.gt.25) ssd=ssd*(2.**(noncon-25))
  622.              endif
  623.              xx= xxd/ssd*3  + xxf/ssf*2
  624.              yy= yyd/ssd*3  + yyf/ssf*2
  625. *  normalisation
  626.           ss=sqrt(xx**2+yy**2)
  627.           if (ss.lt.xpetit/xzprec.or.(ssd.ge.xgrand.and.
  628.      >          ssf.ge.xgrand)) then
  629. *  pas de deplacements, pas de forces
  630.               jcpr(lxf1)=0
  631.               jcpr(lxf2)=0
  632.               xdcp(lxf1)=0.d0
  633.               xdcp(lxf2)=0.d0
  634.               goto 635
  635.           endif
  636. ** on met une grande excitation a la premiere iteration pour bloquer les frottements
  637. ** autant qu'on peut
  638.           if(noncon.le.1) ss=ss/16
  639.           xx=xx/(ss+xpetit)
  640.           yy=yy/(ss+xpetit)
  641. **        write(6,*) 'depl forc ',xxd,yyd,xxf,yyf,xx,yy
  642.           ip=ip+1
  643.           ipt7.num(1,ip)=ipr
  644.           vpocha(ip,1)= -ycpr(lxf1)*xx
  645. ** rajouter xzprec pour passer en NS
  646.           if (ns3d) vpocha(ip,1)=vpocha(ip,1)*xzprec
  647.           xdcp(lxf1)= xx
  648.           ip=ip+1
  649.           ipt7.num(1,ip)=idu
  650. ** rajouter xzprec pour passer en NS
  651.           vpocha(ip,1)= -ycpr(lxf2)*yy
  652.           if (ns3d) vpocha(ip,1)=vpocha(ip,1)*xzprec
  653.           xdcp(lxf2)= yy
  654.           if(ycpr(lxf1).ne.ycpr(lxf2)) write(6,*) 'ycpr differents',
  655.      >      ycpr(lxf1),ycpr(lxf2)
  656. *          write(ioimp,*) ' force ',vpocha(ip,1),ipr,idu
  657.         endif
  658.  600  continue
  659.  635  continue
  660.  631  continue
  661.  610  continue
  662.       nbelem=ip
  663.       n=ip
  664.       segadj ipt7,mpoval
  665.  
  666. **    goto 1954
  667. *
  668. *      maintenant, en cas de frottement, on construit la partie non symetrique de la raideur.
  669. *
  670.       segact mrigid
  671.       nrigel=irigel(/2)
  672.       iforis=iforig
  673.       segini,ri1
  674.       ri1.iforig=iforis
  675.       ri1.mtymat='FROTTEME'   
  676.       ir2=0
  677.       do 1200 ir=1,irigel(/2)
  678.        meleme=irigel(1,ir)
  679.        segact meleme
  680.        descr=irigel(3,ir)
  681.        nbsous=0
  682.        nbref=0
  683.        nbnn=num(/1)
  684.        nbelem=num(/2)
  685.        segini,ipt1
  686.        ipt1.itypel=itypel
  687.        iel2=0
  688.        xmatri=irigel(4,ir)
  689.        segact xmatri
  690.        nligrd=re(/1)
  691.        nligrp=re(/2)
  692.        nelrig=re(/3)
  693.        xmatr1=xmatri
  694. **       segini,xmatr1
  695.        nbnn=num(/1)
  696.        do 1100 iel=1,num(/2)
  697.         lxc=ncpr(num(1,iel))
  698.         if (lxc.eq.0) goto 1100
  699.         ipt7=imolx(lxc)
  700.         if (ipt7.eq.0) goto 1100
  701.         segact ipt7
  702.         nbnn7=ipt7.num(/1)
  703.         ilm=imllx(lxc)
  704.         if(ilm.eq.0) goto 1100
  705.         if(xmatr1.eq.xmatri) segini xmatr1
  706. *  verif que c'est bien un multiplicateur de contact
  707.         if (ncpr(ipt7.num(1,ilm)).ne.lxc) goto 1100
  708. *  et qu'on a determine le sens de frottement
  709.         lxf1=ncpr(ipt7.num(nbnn7,ilm))
  710.         if (idim.eq.2) then
  711.           if (.not. ns2d) goto 1100
  712. *pv        if(jcpr(lxf1).eq.0) then
  713. *pv         goto 1100
  714. *pv        endif
  715.         else
  716. *  test si NS
  717. *  pas de matrice non symetrique en 3D pour le moment
  718.           if (.not. ns3d) goto 1100
  719.           lxf2=ncpr(ipt7.num(nbnn7-1,ilm))
  720. **        if(jcpr(lxf1).eq.0.and.jcpr(lxf2).eq.0) goto 1100
  721.         endif
  722.         iel2=iel2+1
  723.         do ip=1,num(/1)
  724.          ipt1.num(ip,iel2)=num(ip,iel)
  725.         enddo
  726.         xmatr6=irilx(lxf1)
  727.         des6  =irides(lxf1)
  728.         segact xmatr6
  729.         if(xmatr6.re(/1).ne.re(/1)) call erreur(5)
  730.         if(des6.ne.descr) call erreur(5)
  731.         nelri1=iellx(lxf1)
  732. ** xdcp contient le coef de frottement pour les multiplicateurs de contact
  733. ** et l'angle (cos sin) pour les multiplicateurs de frottement
  734.         do id=1,re(/1)
  735.          xmatr1.re(id,1,iel2)=
  736.      >        +xdcp(lxc)*xdcp(lxf1)*xmatr6.re(id,1,nelri1)
  737.         enddo
  738.         if(idim.eq.3) then
  739.           xmatr6=irilx(lxf2)
  740.           segact xmatr6
  741.           nelri2=iellx(lxf2)
  742.           do id=1,re(/1)
  743.            xmatr1.re(id,1,iel2)=xmatr1.re(id,1,iel2)
  744.      >         +xdcp(lxc)*xdcp(lxf2)*xmatr6.re(id,1,nelri2)
  745.           enddo
  746.         endif
  747. **      write(6,*) 'mu ',lxc,xdcp(lxc),'angle',lxf,xdcp(lxf)
  748.  1100  continue
  749.        nbnn=ipt1.num(/1)
  750.        nbelem=iel2
  751.        nsous=0
  752.        nbref=0
  753.        segadj ipt1
  754.        nelrig=iel2
  755.        nligrd=xmatr1.re(/1)
  756.        nligrp=xmatr1.re(/2)
  757.        if(nelrig.ne.0) then
  758.        segadj xmatr1
  759.         ir2=ir2+1
  760.         ri1.coerig(ir2)=coerig(ir)
  761.         ri1.irigel(1,ir2)=ipt1
  762.         ri1.irigel(2,ir2)=irigel(2,ir)
  763.         ri1.irigel(3,ir2)=irigel(3,ir)
  764.         ri1.irigel(4,ir2)=xmatr1
  765.         ri1.irigel(5,ir2)=irigel(5,ir)
  766.         ri1.irigel(6,ir2)=irigel(6,ir)
  767.         ri1.irigel(7,ir2)=2
  768.         ri1.irigel(8,ir2)=irigel(8,ir)
  769.         xmatr1.symre=2
  770.        endif
  771.  1200 continue
  772.       nrigel=ir2
  773.       segadj ri1
  774. 1954  continue
  775. ** ecriture de la raideur de contact non associee
  776. **    call prrigi(ri1,0)
  777.       call ecrobj('RIGIDITE',ri1)
  778. *
  779. *  on construit maintenant la raideur de passage des multiplicateurs de contact aux multiplicateurs de frottement
  780. *  elle permettre en fin de unilater de recalculer les multiplicateurs de frottement du probleme sans la raideur additionnelle
  781.        nligrp=idim
  782.        nligrd=idim
  783.        nelrig=nbptl
  784.        segini xmatri
  785.        nbsous=0
  786.        nbref=0
  787.        nbnn=idim
  788.        nbelem=nbptl
  789.        iel2=0
  790.        segini meleme
  791.        itypel=idim
  792. *  Test si 3D NS
  793.        if ((.not.ns3d).and.idim.eq.3) goto 2101
  794.        if ((.not.ns2d).and.idim.eq.2) goto 2101
  795.        do 2100 lxc=1,nbptl
  796.         ipt7=imolx(lxc)
  797.         if (ipt7.eq.0) goto 2100
  798.         segact ipt7
  799.         ilm=imllx(lxc)
  800.         if(ilm.eq.0) goto 2100
  801. *  verif que c'est bien un multiplicateur de contact
  802.         if (ncpr(ipt7.num(1,ilm)).ne.lxc) goto 2100
  803.         nbnn7=ipt7.num(/1)
  804.         iel2=iel2+1
  805.          num(1,iel2)=ipt7.num(1,ilm)
  806.          num(2,iel2)=ipt7.num(nbnn7,ilm)
  807.          re(1,1,iel2)=0.d0
  808.          lxf1=ncpr(num(2,iel2))
  809. * verifier qu'on a bien mis la matrice additionnelle
  810.         if(idim.eq.2) then
  811. *pv       if(jcpr(lxf1).eq.0) then
  812. **         write(6,*) ' pas de matrice additionnelle '
  813. *pv        iel2=iel2-1
  814. *pv        goto 2100
  815. *pv       endif
  816. **        write (6,*) ' matrice de transfert pour ',
  817. **   >    num(2,iel2),jcpr(lxf1)
  818.         else
  819.           num(3,iel2)=ipt7.num(nbnn7-1,ilm)
  820.           lxf2=ncpr(num(3,iel2))
  821. *pv       if(jcpr(lxf1).eq.0.and.jcpr(lxf2).eq.0) then
  822. *pv        iel2=iel2-1
  823. *pv        goto 2100
  824. *pv       endif
  825.         endif
  826. ** xdcp contient le coef de frottement pour les multiplicateurs de contact
  827. ** et l'angle (cos sin) pour les multiplicateurs de frottement
  828. **       write(6,*) 'excfro transfert',
  829. **   >         num(2,iel2),xdcp(lxf1),num(1,iel2),xdcp(lxc)
  830.          re(2,1,iel2)= xdcp(lxf1) *  xdcp(lxc)
  831.          if(idim.eq.3) then
  832.           re(3,1,iel2)= xdcp(lxf2)*xdcp(lxc)
  833. **       write(6,*) ' angle du frottement ',ipt7.num(1,ilm),
  834. **   >     atan2(xdcp(lxf1),xdcp(lxf2))
  835.          endif
  836.  2100  continue
  837.  2101  continue
  838.        segini descr
  839.         lisinc(1)='LX'
  840.         lisinc(2)='LX'
  841.         lisdua(1)='LX'
  842.         lisdua(2)='LX'
  843.         noelep(1)=1
  844.         noelep(2)=2
  845.         noeled(1)=1
  846.         noeled(2)=2
  847.         if(idim.eq.3) then 
  848.          lisinc(3)='LX'
  849.          lisdua(3)='LX'
  850.          noelep(3)=3
  851.          noeled(3)=3
  852.         endif
  853.        nelrig=iel2
  854.        segadj xmatri
  855.        nbelem=iel2
  856.        segadj meleme
  857.        nrigel=0
  858.        if(nelrig.ne.0) nrigel=1
  859.        segini mrigid
  860.        iforig=iforis
  861.        if(nrigel.eq.1) then
  862.        coerig(1)=1.d0
  863.        irigel(1,1)=meleme
  864.        irigel(3,1)=descr
  865.        irigel(4,1)=xmatri
  866.        irigel(7,1)=2
  867.        symre=2
  868.        endif
  869.  
  870. ** ecriture de la raideur de transfert du contact en frottement
  871. **    call prrigi(mrigid,0)
  872.       call ecrobj('RIGIDITE',mrigid)
  873.  
  874.  
  875.  
  876.  9000 CONTINUE
  877.       call actobj('CHPOINT ',mchpoi,1)
  878. **    write(6,*) 'sortie de excfro '
  879. **    call ecchpo(mchpoi,0)
  880.       call ecrobj('CHPOINT ',mchpoi)
  881.  
  882.  
  883.  9900 CONTINUE
  884.       segsup icpr,jcpr,xcpr,ycpr
  885.                      segsup zcprd,zcprf,kcpr
  886.       segsup ncpr,irilx,imolx,iellx,imllx,irides
  887.  
  888.  
  889.  9999 CONTINUE
  890.       mmode2 = mocoul
  891.       SEGSUP mmode2
  892.  
  893.       END
  894.  
  895.  
  896.  
  897.  
  898.  
  899.  
  900.  
  901.  
  902.  
  903.  
  904.  
  905.  
  906.  
  907.  
  908.  
  909.  
  910.  
  911.  
  912.  
  913.  
  914.  
  915.  
  916.  
  917.  
  918.  

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