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