1!$Id:$ 2 subroutine frame2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 3 4! * * F E A P * * A Finite Element Analysis Program 5 6!.... Copyright (c) 1984-2017: Regents of the University of California 7! All rights reserved 8 9!-----[--.----+----.----+----.-----------------------------------------] 10 11! Two dimensional frame element 12 13! Control Information: 14 15! ndm - Spatial dimension of problem = 2 16! ndf - Number degree-of-freedoms at node >= 3 17! ( 1 = u_1 ; 2 = u_2 ; 3 = theta ) 18! nen - Two node element (nel = 2) >= 2 19 20!-----[--.----+----.----+----.-----------------------------------------] 21 22 implicit none 23 24 include 'bdata.h' 25 include 'cdata.h' 26 include 'eldata.h' 27 include 'eltran.h' 28 include 'hdata.h' 29 include 'iofile.h' 30 include 'prld1.h' 31 32 include 'comblk.h' 33 34 integer ndf,ndm,nst,isw,tdof,i 35 real*8 cs,sn,le 36 37 integer ix(*) 38 real*8 d(*),ul(ndf,nen,*),xl(ndm,*),s(nst,*),p(*) 39 40 save 41 42! INPUT MATERIAL PARAMETERS 43 44 if(isw.eq.1) then 45 if(ior.lt.0) write(*,2000) 46 write(iow,2000) 47 call inmate(d,tdof,3*2,3) 48 49! Set plot sequence 50 51 pstyp = 1 52 53! Check dimensions 54 55 if(ndm.ne.2 .or. ndf.lt.3 .or. nen.lt.2) then 56 write(iow,3000) ndm,ndf,nen 57 call plstop(.true.) 58 endif 59 60! Deactivate dof in element for dof > 3 61 62 do i = 4,ndf 63 ix(i) = 0 64 end do 65 66! Check for geometric storage for finite deformation element 67 68 if(d(18).lt.0.0d0 ) then 69 if(d(79).eq.0.0d0) then 70 call framf2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 71 else 72 call franf2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 73 endif 74 endif 75 76! CHECK ELEMENT FOR ERRORS 77 78 elseif(isw.eq.2) then 79 80 cs = xl(1,2) - xl(1,1) 81 sn = xl(2,2) - xl(2,1) 82 le = sqrt(cs*cs+sn*sn) 83 84 if(ix(1).le.0 .or. ix(2).le.0 .or. ix(1).eq.ix(2)) then 85 write(iow,4000) n,ix(1),ix(2) 86 if(ior.lt.0) write(*,4000) n,ix(1),ix(2) 87 endif 88 if(le.le.0.0d0) then 89 write(iow,4001) n 90 if(ior.lt.0) write(*,4001) n 91 endif 92 93! COMPUTE ELEMENT ARRAYS 94 95 else 96 97! Small deformaion 98 99 if(d(18).gt.0.0d0 ) then 100 101! Shear deformable (2-node: linear interpolations) 102 103 if(d(79).eq.0.0d0) then 104 call frams2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 105 106! Euler-Bernoulli (2-node: cubic interpolations) 107 108 else 109 call frans2d(d,ul,xl,s,p,ndf,ndm,nst,isw) 110 endif 111 112! Finite deformation (2-node: linear interpolations) 113 114 else 115 116! Shear deformable (2-node: linear, finite displacements) 117 118 if(d(79).eq.0.0d0) then 119 120 call framf2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 121 122! No shear case (2-node: cubic, 2-nd order displacements) 123 124 else 125 126 call franf2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 127 128 endif 129 endif 130 endif 131 132! Formats 133 1342000 format(5x,'T w o D i m e n s i o n a l F r a m e'/) 135 1363000 format(/' *ERROR* Problem control values incorrect. Set as:'/ 137 & ' ndm = ',i2,': (Should be 2)'/ 138 & ' ndf = ',i2,': (Minimum = 3)'/ 139 & ' nen = ',i2,': (Minimum = 2)'/) 140 1414000 format(' *ERROR* Element',i7,' has unspecified node: ix =',2i7) 1424001 format(' *ERROR* Element',i7,' has zero length') 143 144 end 145 146 subroutine frams2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 147 148! Purpose: Two dimensional small displacement frame element 149 150 implicit none 151 152 include 'bdata.h' 153 include 'bm2com.h' 154 include 'cdata.h' 155 include 'cdat1.h' 156 include 'eldata.h' 157 include 'eltran.h' 158 include 'hdata.h' 159 include 'iofile.h' 160 include 'prld1.h' 161 include 'prstrs.h' 162 include 'ptdat6.h' 163 include 'tdata.h' 164 165 include 'comblk.h' 166 167 integer ndf,ndm,nst,isw 168 integer i,ii, j,jj, ll,lint, mm,nn 169 real*8 cs,sn,a1,a2,a3,a4,le,b1,b2,dx,dva,dvi,xjac, energy 170 real*8 rhoa, rhoi, ctan1, ctan3 171 172 integer ix(*) 173 174 real*8 aa(4,4,2),shp(2,3),sg(3),wg(3),forc(4,2),xx(2),b(2) 175 real*8 xl(ndm,*),ul(ndf,nen,*) 176 real*8 d(*),p(ndf,*),s(nst,nst) 177 178 save 179 180! Small deformation visco-plastic BEAM element. 181 182! d(21)*d(32) = EA 183! d(37)*d(27)*d(32) = kappa*GA 184! d(21)*d(33) = EI 185! d( 4)*d(32) = rho*A 186! d( 4)*d(33) = rho*I 187! 188 189! Compute element length and direction cosines 190 191 if(isw.ge.2) then 192 cs = xl(1,nel) - xl(1,1) 193 sn = xl(2,nel) - xl(2,1) 194 le = sqrt(cs*cs + sn*sn) 195 cs = cs/le 196 sn = sn/le 197 198 rhoa = d(4)*d(32) 199 rhoi = d(4)*d(33) 200 endif 201 202! Read data 203 204 if(isw.eq.1) then 205 206! Increment history storage if necessary 207 208 elseif(isw.eq.3 .or. isw.eq.6) then 209 210 lint = nel - 1 211 call int1d(lint,sg,wg) 212 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 213 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 214 call bm2trn (ul(1,1,5),cs,sn,ndf*nel,ndf,2) 215 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 216 dva = 0.5d0*rhoa*le 217 dvi =0.5d0*rhoi*d(69)*le 218 do ll = 1,lint 219 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 220 dx = wg(ll)*xjac 221 call b2mods (d,ul,forc,aa,energy,shp,ndf,nen,isw) 222 223! Multiply moduli by solution parameter: ctan(1) 224 225 ctan1 = ctan(1) 226 do jj = 1,4 227 do ii = 1,4 228 aa(ii,jj,1) = aa(ii,jj,1)*ctan1 229 end do ! ii 230 end do ! jj 231 232! Mechanical tangent terms 233 234 mm = 0 235 do ii = 1,nel 236 b1 = shp(1,ii)*dx 237 b2 = shp(2,ii)*dx 238 if(isw.eq.3) then 239 nn = 0 240 do jj = 1,nel 241 a1 = b1*shp(1,jj) 242 a2 = b1*shp(2,jj) 243 a3 = b2*shp(1,jj) 244 a4 = b2*shp(2,jj) 245 do i = 1,3 246 s(i+mm,3+nn) = s(i+mm,3+nn) + a2*aa(i,4,1) 247 s(3+mm,i+nn) = s(3+mm,i+nn) + a3*aa(4,i,1) 248 do j = 1,3 249 s(i+mm,j+nn) = s(i+mm,j+nn) + a1*aa(i,j,1) 250 end do ! j 251 end do ! i 252 s(3+mm,3+nn) = s(3+mm,3+nn) + a4*aa(4,4,1) 253 nn = nn + ndf 254 end do ! jj 255 endif 256 do i = 1,3 257 p(i,ii) = p(i,ii) - forc(i,1)*b1 258 end do ! i 259 p(3,ii) = p(3,ii) - forc(4,1)*b2 260 mm = mm + ndf 261 end do ! ii 262 end do ! ll 263 264! Lumped and consistent inertia contributions 265 266 if(nint(d(7)).eq.1) then 267 b1 = 2.d0/3.d0 268 b2 = 0.5d0 269 else 270 b1 = 1.0d0 271 b2 = 0.0d0 272 endif 273 274! Form diagonal terms 275 ctan3 = ctan(3) + d(77)*ctan(2) 276 jj = 0 277 do ii = 1,nel 278 p(1,ii) = p(1,ii) - dva*(ul(1,ii,5)+d(77)*ul(1,ii,4)) 279 p(2,ii) = p(2,ii) - dva*(ul(2,ii,5)+d(77)*ul(2,ii,4)) 280 p(3,ii) = p(3,ii) - dvi*(ul(3,ii,5)+d(77)*ul(3,ii,4)) 281 s(jj+1,jj+1) = s(jj+1,jj+1) + dva*ctan3*b1 282 s(jj+2,jj+2) = s(jj+2,jj+2) + dva*ctan3*b1 283 s(jj+3,jj+3) = s(jj+3,jj+3) + dvi*ctan3*b1 284 jj = jj + ndf 285 end do 286 287! Off diagonal terms for consistent mass 288 289 do ii = 1,3 290 s(ii,ndf+ii) = b2*s(ii,ii) 291 s(ndf+ii,ii) = b2*s(ii,ii) 292 end do ! ii 293 294! Transform stiffness and residual to global coordinates 295 296 if(isw.eq.3) then 297 call bm2trn (s,cs,sn,nst,ndf,1) 298 endif 299 call bm2trn ( p,cs,-sn,nst,ndf,2) 300 301! Set body loading factors 302 303 if(int(d(74)).gt.0) then 304 b(1) = 0.5*le*(d(11) + prldv(int(d(74)))*d(71)) 305 else 306 b(1) = 0.5*le*d(11)*dm 307 endif 308 if(int(d(75)).gt.0) then 309 b(2) = 0.5*le*(d(12) + prldv(int(d(75)))*d(72)) 310 else 311 b(2) = 0.5*le*d(12)*dm 312 endif 313 314! Add body force components 315 316 p(1,1) = p(1,1) + b(1) 317 p(1,2) = p(1,2) + b(1) 318 p(2,1) = p(2,1) + b(2) 319 p(2,2) = p(2,2) + b(2) 320 321! Output forces 322 323 elseif(isw.eq.4 .or. isw.eq.8) then 324 lint = nel - 1 325 call int1d(lint,sg,wg) 326 327! Loop over quadrature points 328 329 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 330 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 331 call bm2trn (ul(1,1,5),cs,sn,ndf*nel,ndf,2) 332 call bm2trn (xl,cs,sn,ndm*nel,ndm,2) 333 do ll = 1,lint 334 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 335 336! Output forces 337 338 call b2mods (d,ul,forc,aa,energy,shp,ndf,nen,isw) 339 if(isw.eq.4) then 340 do i = 1,ndm 341 xx(i) = 0. 342 do ii = 1,nel 343 xx(i) = xx(i) + xl(i,ii)*shp(2,ii) 344 end do ! ii 345 end do ! i 346 mct = mct - 3 347 if (mct.le.0) then 348 write(iow,2001) o,head,ttim 349 if(ior.lt.0) write(*,2001) o,head,ttim 350 mct = 50 351 endif 352 write(iow,2002) n,ma,(xx(i),i=1,2), 353 & (strs(i,1),i=1,3),(defa(i,1),i=1,3) 354 if(ior.lt.0) then 355 write(*,2002) n,ma,(xx(i),i=1,2), 356 & (strs(i,1),i=1,3),(defa(i,1),i=1,3) 357 endif 358 359! Stress projections save 360 361 else 362 363 dx = wg(ll)*xjac 364 do ii = 1,nel 365 b1 = shp(1,ii)*dx 366 do i = 1,3 367 p(i,ii) = p(i,ii) - forc(i,1)*b1 368 end do ! i 369 p(3,ii) = p(3,ii) - forc(4,1)*shp(2,ii)*dx 370 end do ! ii 371 372 endif 373 374 end do ! ll 375 376 if(isw.eq.8) then 377 do i = 1,3 378 p(i,2) = -p(i,2) 379 end do 380 call frcn2d(ix,p,ndf,hr(nph),hr(nph+numnp),numnp) 381 endif 382 383! Compute mass array 384 385 elseif(isw.eq.5) then 386 387 lint = nel 388 call int1d(lint,sg,wg) 389 390! Compute lumped mass matrix 391 392 call bm2trn (xl,cs,sn,ndm*nel,ndm,2) 393 do ll = 1,lint 394 395! Compute shape functions 396 397 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 398 399 dva = wg(ll)*xjac*rhoa 400 dvi = wg(ll)*xjac*rhoi*d(69) 401 402! For each node j compute db = rho*shape*dv 403 404 do j = 1,nel 405 b1 = shp(2,j)*dva 406 b2 = shp(2,j)*dvi 407 408! Compute a lumped mass 409 410 p(1,j) = p(1,j) + b1 411 p(2,j) = p(2,j) + b1 412 p(3,j) = p(3,j) + b2 413 end do ! j 414 end do ! ll 415 416! Place in consistent mass 417 418 if(nint(d(7)).eq.1) then 419 b1 = 2.d0/3.d0 420 b2 = 0.5d0 421 else 422 b1 = 1.0d0 423 b2 = 0.0d0 424 endif 425 426 jj = 0 427 do j = 1,nel 428 s(jj+1,jj+1) = p(1,j)*b1 429 s(jj+2,jj+2) = p(2,j)*b1 430 s(jj+3,jj+3) = p(3,j)*b1 431 jj = jj + ndf 432 end do ! j 433 434! Off diagonal consistent mass terms 435 436 do ii = 1,3 437 s(ii,ndf+ii) = b2*s(ii,ii) 438 s(ndf+ii,ii) = b2*s(ii,ii) 439 end do ! ii 440 441 elseif(isw.eq.12) then 442 443 lint = nel - 1 444 call int1d(lint,sg,wg) 445 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 446 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 447 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 448 do ll = 1,lint 449 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 450 call b2mods (d,ul,forc,aa,energy,shp,ndf,nen,isw) 451 end do ! ll 452 453 454! Compute energy 455 456 elseif(isw.eq.13) then 457 458 lint = nel - 1 459 call int1d(lint,sg,wg) 460 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 461 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 462 call bm2trn (ul(1,1,4),cs,sn,ndf*nel,ndf,2) 463 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 464 dva = 0.5d0*rhoa*a1 465 dvi =0.5d0*rhoi*d(69)*a1 466 467! Compute internal energy 468 469 do ll = 1,lint 470 471! Compute energy density from stress and deformation 472 473 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 474 dx = wg(ll)*xjac 475 call b2mods (d,ul,forc,aa,energy,shp,ndf,nen,isw) 476 477! Accumulate energy 478 479 epl(8) = epl(8) + 0.5d0*energy*dx 480 481 end do ! ll 482 483! Compute kinetic energy for lumped mass 484 485 epl(7) = epl(7) + 0.5d0*dva*(ul(1,1,4)**2 + ul(1,2,4)**2 486 & + ul(2,1,4)**2 + ul(2,2,4)**2) 487 & + 0.5d0*dvi*(ul(3,1,4)**2 + ul(3,2,4)**2) 488 489 endif 490 491! Formats 492 4932001 format(a1,20a4/5x,'time',e13.5,5x,' element forces '// 494 & 43x,'********* FORCE / STRAIN *********'/ 495 & 3x,'element material', 496 & 3x,'1-coord',3x,'2-coord',6x,'n-dir',8x,'s-dir',8x,'m-dir'/) 497 4982002 format(2i10,0p,2f10.3,1p,3e13.4/40x,1p,3e13.4) 499 500 end 501 502 subroutine b2mods (d,ul,forca,aa,energy,shp,ndf,nen,isw) 503 504 implicit none 505 506 include 'bm2com.h' 507 include 'ddata.h' 508 include 'eldata.h' 509 include 'tdata.h' 510 511 integer ndf,nen, i,ii,isw 512 513 real*8 d(*),ul(ndf,nen,*),cc(3,3,2) 514 real*8 forca(4,2),aa(4,4,2), shp(2,3), energy 515 516 save 517 518! Compute beam strains 519 520 do i = 1,3 521 defa(i,1) = 0.0d0 522 defa(i,2) = 0.0d0 523 do ii = 1,nel 524 defa(i,1) = defa(i,1) + ul(i,ii,1)*shp(1,ii) 525 defa(i,2) = defa(i,2) + ul(i,ii,4)*shp(1,ii) 526 end do ! ii 527 end do ! i 528 529 do ii = 1,nel 530 defa(2,1) = defa(2,1) - ul(3,ii,1)*shp(2,ii) 531 defa(2,2) = defa(2,2) - ul(3,ii,4)*shp(2,ii) 532 end do ! ii 533 534! Compute forces 535 536 call bm2con (d,cc,strs,defa,isw) 537 538! Compute stored energy density 539 540 if(isw.eq.13) then 541 542 energy = strs(1,1)*defa(1,1) 543 & + strs(2,1)*defa(2,1) 544 & + strs(3,1)*defa(3,1) 545 546 elseif(isw.ne.12) then 547 548 do ii = 1,2 549 550! Compute first Piola-material frame 551 552 forca(1,ii) = strs(1,ii) 553 forca(2,ii) = strs(2,ii) 554 forca(3,ii) = strs(3,ii) 555 forca(4,ii) = -forca(2,ii) 556 557! Compute tangent tensor 558 559 aa(1,1,ii) = cc(1,1,ii) 560 aa(1,2,ii) = cc(1,2,ii) 561 aa(1,3,ii) = cc(1,3,ii) 562 563 aa(2,1,ii) = cc(1,2,ii) 564 aa(2,2,ii) = cc(2,2,ii) 565 aa(2,3,ii) = cc(2,3,ii) 566 567 aa(3,1,ii) = cc(1,3,ii) 568 aa(3,2,ii) = cc(2,3,ii) 569 aa(3,3,ii) = cc(3,3,ii) 570 571 aa(1,4,ii) = -cc(1,2,ii) 572 aa(2,4,ii) = -cc(2,2,ii) 573 aa(3,4,ii) = -cc(2,3,ii) 574 575 aa(4,1,ii) = aa(1,4,ii) 576 aa(4,2,ii) = aa(2,4,ii) 577 aa(4,3,ii) = aa(3,4,ii) 578 579 aa(4,4,ii) = cc(2,2,ii) 580 end do ! ii 581 582 endif 583 584 end 585 586 subroutine framf2d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) 587 588! Purpose: 2-d finite deformation frame element 589 590 implicit none 591 592 include 'augdat.h' 593 include 'bdata.h' 594 include 'bm2com.h' 595 include 'cdata.h' 596 include 'cdat1.h' 597 include 'eldata.h' 598 include 'eltran.h' 599 include 'hdata.h' 600 include 'iofile.h' 601 include 'prld1.h' 602 include 'prstrs.h' 603 include 'ptdat6.h' 604 include 'tdata.h' 605 606 include 'comblk.h' 607 608 integer ndf,ndm,nst,isw 609 integer i,ii, j,jj,j1, ll,lint, mm,nn 610 real*8 cs,sn,a1,a2,a3,a4,b1,b2,dx,dva,dvi,xjac, energy 611 real*8 rhoa, rhoi 612 613 integer ix(*) 614 615 real*8 a(4,4),shp(2,3),sg(3),wg(3),forc(4),xx(2) 616 real*8 xl(ndm,*),ul(ndf,nen,*) 617 real*8 d(*),p(ndf,*),s(nst,nst) 618 619 save 620 621! Finite deformation visco-plastic BEAM element. 622 623! d(1)*d(32) = EA 624! d(37)*d(2)*d(32)/2/(1+d(2)) = kappa*GA 625! d(1)*d(33) = EI 626! d(4)*d(32) = rho*A 627! d(4)*d(33) = rho*I 628 629 if(isw.eq.1) then 630 if(nint(d(160)).eq.2) then 631 nh1 = nh1 + 1 ! Augmented storage 632 d(166) = 1 633 else 634 d(166) = 0 635 endif 636 637! Compute element length and direction cosines 638 639 elseif(isw.ge.2) then 640 cs = xl(1,nel) - xl(1,1) 641 sn = xl(2,nel) - xl(2,1) 642 a1 = sqrt(cs*cs + sn*sn) 643 cs = cs/a1 644 sn = sn/a1 645 646 rhoa = d(4)*d(32) 647 rhoi = d(4)*d(33) 648 endif 649 650! Read data 651 652 if(isw.eq.1) then 653 654! History storage for rotation parameters 655 656 nh1 = nh1 + 2 657 658 elseif(isw.eq.3 .or. isw.eq.6) then 659 660 lint = nel - 1 661 call int1d(lint,sg,wg) 662 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 663 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 664 call bm2trn (ul(1,1,5),cs,sn,ndf*nel,ndf,2) 665 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 666 dva = 0.5d0*rhoa*a1 667 dvi =0.5d0*rhoi*d(69)*a1 668 do ll = 1,lint 669 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 670 dx = wg(ll)*xjac 671 call b2modl (d,ul,forc,a,energy,shp,ndf,nen,isw) 672 673! Multiply moduli by solution parameter: ctan(1) 674 do jj = 1,4 675 do ii = 1,4 676 a(ii,jj) = a(ii,jj)*ctan(1) 677 end do ! ii 678 end do ! jj 679 680! Compute residual and tangent 681 682! Mechanical and Geometric tangent terms 683 684 mm = 0 685 do ii = 1,nel 686 b1 = shp(1,ii)*dx 687 b2 = shp(2,ii)*dx 688 if(isw.eq.3) then 689 nn = 0 690 do jj = 1,nel 691 a1 = b1*shp(1,jj) 692 a2 = b1*shp(2,jj) 693 a3 = b2*shp(1,jj) 694 a4 = b2*shp(2,jj) 695 do i = 1,3 696 s(i+mm,3+nn) = s(i+mm,3+nn) + a2*a(i,4) 697 s(3+mm,i+nn) = s(3+mm,i+nn) + a3*a(4,i) 698 do j = 1,3 699 s(i+mm,j+nn) = s(i+mm,j+nn) + a1*a(i,j) 700 end do ! j 701 end do ! i 702 s(3+mm,3+nn) = s(3+mm,3+nn) + a4*a(4,4) 703 nn = nn + ndf 704 end do ! jj 705 endif 706 do i = 1,3 707 p(i,ii) = p(i,ii) - forc(i)*b1 708 end do ! i 709 p(3,ii) = p(3,ii) - forc(4)*b2 710 mm = mm + ndf 711 end do ! ii 712 end do ! ll 713 714! Lumped inertia contributions 715 716 jj = 0 717 do ii = 1,nel 718 p(1,ii) = p(1,ii) - dva*ul(1,ii,5) 719 p(2,ii) = p(2,ii) - dva*ul(2,ii,5) 720 p(3,ii) = p(3,ii) - dvi*ul(3,ii,5) 721 s(jj+1,jj+1) = s(jj+1,jj+1) + dva*ctan(3) 722 s(jj+2,jj+2) = s(jj+2,jj+2) + dva*ctan(3) 723 s(jj+3,jj+3) = s(jj+3,jj+3) + dvi*ctan(3) 724 jj = jj + ndf 725 end do 726 727! Transform stiffness and residual to global coordinates 728 729 if(isw.eq.3) then 730 call bm2trn (s,cs,sn,nst,ndf,1) 731 endif 732 call bm2trn ( p,cs,-sn,nst,ndf,2) 733 734! Set body loading factors 735 736 if(int(d(74)).gt.0) then 737 b1 = d(11) + prldv(int(d(74)))*d(71) 738 else 739 b1 = d(11)*dm 740 endif 741 if(int(d(75)).gt.0) then 742 b2 = d(12) + prldv(int(d(75)))*d(72) 743 else 744 b2 = d(12)*dm 745 endif 746 747! Add body forces 748 749 if(nel.eq.2) then 750 751 p(1,1) = p(1,1) + b1*xjac 752 p(2,1) = p(2,1) + b2*xjac 753 p(1,2) = p(1,2) + b1*xjac 754 p(2,2) = p(2,2) + b2*xjac 755 756 else 757 do ll = 1,lint 758 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 759 dx = dm*wg(ll)*xjac 760 do ii = 1,nel 761 a1 = shp(2,ii)*dx 762 p(1,ii) = p(1,ii) + a1*b1 763 p(2,ii) = p(2,ii) + a1*b2 764 end do ! ii 765 end do ! ll 766 endif 767 768! Output forces 769 770 elseif(isw.eq.4 .or. isw.eq.8) then 771 lint = nel - 1 772 call int1d(lint,sg,wg) 773 774! Loop over quadrature points 775 776 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 777 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 778 call bm2trn (ul(1,1,5),cs,sn,ndf*nel,ndf,2) 779 call bm2trn (xl,cs,sn,ndm*nel,ndm,2) 780 do ll = 1,lint 781 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 782 783! Output forces 784 785 call b2modl (d,ul,forc,a,energy,shp,ndf,nen,isw) 786 if(isw.eq.4) then 787 do i = 1,ndm 788 xx(i) = 0. 789 do ii = 1,nel 790 xx(i) = xx(i) + xl(i,ii)*shp(2,ii) 791 end do ! ii 792 end do ! i 793 mct = mct - 3 794 if (mct.le.0) then 795 write(iow,2001) o,head,ttim 796 if(ior.lt.0) write(*,2001) o,head,ttim 797 mct = 50 798 endif 799 write(iow,2002) n,ma,(xx(i),i=1,ndm), 800 & (strs(i,1),i=1,3),(defa(i,1),i=1,3) 801 if(ior.lt.0) then 802 write(*,2002) n,ma,(xx(i),i=1,ndm), 803 & (strs(i,1),i=1,3),(defa(i,1),i=1,3) 804 endif 805 806! Stress projections save 807 808 else 809 810 dx = wg(ll)*xjac 811 do ii = 1,nel 812 b1 = shp(1,ii)*dx 813 do i = 1,3 814 p(i,ii) = p(i,ii) - forc(i)*b1 815 end do ! i 816 p(3,ii) = p(3,ii) - forc(4)*shp(2,ii)*dx 817 end do ! ii 818 819 endif 820 821 end do ! ll 822 823 if(isw.eq.8) then 824 do i = 1,3 825 p(i,2) = -p(i,2) 826 end do 827 call frcn2d(ix,p,ndf,hr(nph),hr(nph+numnp),numnp) 828 endif 829 830! Compute mass array 831 832 elseif(isw.eq.5) then 833 834 lint = nel 835 call int1d(lint,sg,wg) 836 837! Compute lumped mass matrix 838 839 call bm2trn (xl,cs,sn,ndm*nel,ndm,2) 840 do ll = 1,lint 841 842! Compute shape functions 843 844 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 845 846 dva = wg(ll)*xjac*rhoa 847 dvi = wg(ll)*xjac*rhoi*d(69) 848 849! For each node j compute db = rho*shape*dv 850 851 j1 = 1 852 do j = 1,nel 853 b1 = shp(2,j)*dva 854 b2 = shp(2,j)*dvi 855 856! Compute a lumped mass 857 858 p(1,j) = p(1,j) + b1 859 p(2,j) = p(2,j) + b1 860 p(3,j) = p(3,j) + b2 861 862 s(j1 ,j1 ) = s(j1 ,j1 ) + b1 863 s(j1+1,j1+1) = s(j1+1,j1+1) + b1 864 s(j1+2,j1+2) = s(j1+2,j1+2) + b2 865 866 j1 = j1 + ndf 867 end do ! j 868 end do ! ll 869 870! Augmented update or time update of history variables 871 872 elseif(isw.eq.10 .or. isw.eq.12) then 873 874 lint = nel - 1 875 call int1d(lint,sg,wg) 876 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 877 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 878 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 879 do ll = 1,lint 880 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 881 call b2modl (d,ul,forc,a,energy,shp,ndf,nen,isw) 882 end do ! ll 883 884! Compute energy 885 886 elseif(isw.eq.13) then 887 888 lint = nel - 1 889 call int1d(lint,sg,wg) 890 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 891 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 892 call bm2trn (ul(1,1,4),cs,sn,ndf*nel,ndf,2) 893 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 894 dva = 0.5d0*rhoa*a1 895 dvi =0.5d0*rhoi*d(69)*a1 896 897! Compute internal energy 898 899 do ll = 1,lint 900 901! Compute energy density from stress and deformation 902 903 call shp1db(sg(ll),xl,shp,ndm,nel,xjac) 904 dx = wg(ll)*xjac 905 call b2modl (d,ul,forc,a,energy,shp,ndf,nen,isw) 906 907! Accumulate energy 908 909 epl(8) = epl(8) + 0.5d0*energy*dx 910 911 end do ! ll 912 913! Compute kinetic energy for lumped mass 914 915 epl(7) = epl(7) + 0.5d0*dva*(ul(1,1,4)**2 + ul(1,2,4)**2 916 & + ul(2,1,4)**2 + ul(2,2,4)**2) 917 & + 0.5d0*dvi*(ul(3,1,4)**2 + ul(3,2,4)**2) 918 919 endif 920 921! Formats 922 9232001 format(a1,20a4/5x,'time',e13.5,5x,' element forces '// 924 & 43x,'********* FORCE / STRAIN *********'/ 925 & 3x,'element material', 926 & 3x,'1-coord',3x,'2-coord',6x,'n-dir',8x,'s-dir',8x,'m-dir'/) 927 9282002 format(2i10,0p,2f10.3,1p,3e13.4/40x,1p,3e13.4) 929 930 end 931 932 subroutine b2modl (d,ul,forca,a,energy,shp,ndf,nen,isw) 933 934 implicit none 935 936 include 'bm2com.h' 937 include 'ddata.h' 938 include 'eldata.h' 939 include 'pmod2d.h' 940 include 'tdata.h' 941 942 integer ndf,nen, i,ii,isw 943 real*8 energy, psi, cn,sn, phia1,phia2 944 945 real*8 d(*),ul(ndf,nen,*),fa(3),df(3),cc(3,3,2) 946 real*8 forca(4),a(4,4),b(2,2),shp(2,3) 947 948 save 949 950! Compute gradient 951 952 do i = 1,3 953 fa(i) = 0.0d0 954 df(i) = 0.0d0 955 do ii = 1,nel 956 fa(i) = fa(i) + ul(i,ii,1)*shp(1,ii) 957 df(i) = df(i) + ul(i,ii,2)*shp(1,ii) 958 end do ! ii 959 end do ! i 960 fa(1) = fa(1) + 1.d0 961 962! Generalized mid-point formulation 963 964 psi = 0.d0 965 do ii = 1,nel 966 psi = psi + ul(3,ii,1)*shp(2,ii) 967 end do ! ii 968 cn = cos(psi) 969 sn = sin(psi) 970 971! Compute strains and forces in gaussian frame 972 973 defa(1,1) = fa(1)*cn + fa(2)*sn - 1.d0 974 defa(2,1) = fa(2)*cn - fa(1)*sn 975 defa(3,1) = fa(3) 976 977 call bm2con (d,cc,strs,defa,isw) 978 979! Compute stored energy density 980 981 if(isw.eq.13) then 982 983 energy = strs(1,1)*defa(1,1) 984 & + strs(2,1)*defa(2,1) 985 & + strs(3,1)*defa(3,1) 986 987 elseif(isw.ne.12) then 988 989! Compute first Piola-material frame 990 991 forca(1) = cn*strs(1,1) - sn*strs(2,1) 992 forca(2) = sn*strs(1,1) + cn*strs(2,1) 993 forca(3) = strs(3,1) 994 forca(4) = fa(2)*forca(1) - fa(1)*forca(2) 995 996! Compute tangent elastic tensor 997 998 b(1,1) = cn*cc(1,1,1) - sn*cc(2,1,1) 999 b(1,2) = cn*cc(1,2,1) - sn*cc(2,2,1) 1000 b(2,1) = sn*cc(1,1,1) + cn*cc(2,1,1) 1001 b(2,2) = sn*cc(1,2,1) + cn*cc(2,2,1) 1002 1003 a(1,1) = b(1,1)*cn - b(1,2)*sn 1004 a(1,2) = b(1,1)*sn + b(1,2)*cn 1005 a(1,3) = cn*cc(1,3,1) - sn*cc(2,3,1) 1006 1007 a(2,1) = a(1,2) 1008 a(2,2) = b(2,1)*sn + b(2,2)*cn 1009 a(2,3) = sn*cc(1,3,1) + cn*cc(2,3,1) 1010 1011 a(3,1) = a(1,3) 1012 a(3,2) = a(2,3) 1013 a(3,3) = cc(3,3,1) 1014 1015 phia1 = defa(1,1) + 1.d0 1016 phia2 = defa(2,1) 1017 1018 a(1,4) = b(1,1)*phia2 - b(1,2)*phia1 1019 a(2,4) = b(2,1)*phia2 - b(2,2)*phia1 1020 a(3,4) = cc(1,3,1)*phia2 - cc(2,3,1)*phia1 1021 1022 a(4,4) = phia2*(cc(1,1,1)*phia2 - cc(1,2,1)*phia1) 1023 & + phia1*(cc(2,2,1)*phia1 - cc(1,2,1)*phia2) 1024 1025 if(gflag) then 1026 a(1,4) = a(1,4) - forca(2) 1027 a(2,4) = a(2,4) + forca(1) 1028 a(4,4) = a(4,4) - fa(1)*forca(1) - fa(2)*forca(2) 1029 endif 1030 1031 a(4,1) = a(1,4) 1032 a(4,2) = a(2,4) 1033 a(4,3) = a(3,4) 1034 1035 endif 1036 1037 end 1038 1039 subroutine bm2con(d, cc,strs,def,isw) 1040 1041 implicit none 1042 1043 include 'hdata.h' 1044 include 'comblk.h' 1045 1046 integer isw, i 1047 1048 real*8 d(*), cc(3,3,2),strs(3,2),def(3,2), dl(3) 1049 1050 save 1051 1052! Resultant model 1053 1054 dl(1) = d(1)*d(32) 1055 dl(2) = d(1)*d(32)*d(37)*0.5d0/(1.d0+d(2)) 1056 dl(3) = d(1)*d(33) 1057 1058 do i = 1,9 1059 cc(i,1,1) = 0.0d0 1060 cc(i,1,2) = 0.0d0 1061 end do ! i 1062 1063! Elastic resultant model only 1064 1065 do i = 1,3 1066 cc(i,i,1) = dl(i) 1067 cc(i,i,2) = dl(i) 1068 strs(i,1) = dl(i)*def(i,1) 1069 strs(i,2) = dl(i)*def(i,2) 1070 end do ! i 1071 1072 if(nint(d(160)).eq.2) then ! Add augmented value 1073 strs(1,1) = strs(1,1) + hr(nh2) 1074 if(isw.eq.10) then ! Update for augmentation 1075 hr(nh2) = strs(1,1) 1076 endif 1077 endif 1078 1079 end 1080 1081 subroutine shp1db(s,xl,shp,ndm,nel,xaj) 1082 1083 implicit none 1084 1085 integer ndm,nel 1086 real*8 s,s2, xaj, hx1,hx2,h,h2, sh 1087 real*8 xl(ndm,nel), shp(2,nel) 1088 1089 save 1090 1091 hx1 = xl(1,nel) - xl(1,1) 1092 hx2 = xl(2,nel) - xl(2,1) 1093 h = sqrt(hx1*hx1 + hx2*hx2) 1094 sh = s*0.5 1095 1096! Linear element 1097 1098 if(nel.eq.2) then 1099 1100 shp(2,1) = 0.5d0 - sh 1101 shp(2,2) = 0.5d0 + sh 1102 xaj = h*0.5d0 1103 shp(1,1) = -1.0d0/h 1104 shp(1,2) = -shp(1,1) 1105 1106! Quadratic element 1107 1108 elseif(nel.eq.3) then 1109 1110 hx1 = xl(1,nel) - xl(1,nel-1) 1111 hx2 = xl(2,nel) - xl(2,nel-1) 1112 h2 = sqrt(hx1*hx1 + hx2*hx2) 1113 s2 = s*s*0.5d0 1114 shp(2,1) = s2 - sh 1115 shp(2,2) = 1.0d0 - s2 - s2 1116 shp(2,3) = s2 + sh 1117 xaj = h*0.5d0 + s*(h2+h2-h) 1118 shp(1,1) = (s-0.5)/xaj 1119 shp(1,2) = -(s+s)/xaj 1120 shp(1,3) = (s+0.5)/xaj 1121 1122 endif 1123 1124 end 1125 1126 subroutine bm2trn(s,cs,sn,nst,ndf,itype) 1127 1128! Itype: 1129 1130! 1 Transform matrix s(nst,nst) 1131! 2 Transform vector s(nst,1) 1132 1133 implicit none 1134 1135 integer nst,ndf,itype, nss 1136 integer i,j,n 1137 real*8 cs,sn,t,tol 1138 1139 real*8 s(nst,*) 1140 1141 save 1142 1143 data tol/ 1.d-12 / 1144 1145 if(cs.gt.1.0d0-tol) return 1146 1147 nss = nst - mod(nst,ndf) 1148 1149 if(itype.eq.1) then 1150 1151 do i = 1,nss,ndf 1152 j = i + 1 1153 do n = 1,nss 1154 t = s(n,i)*cs - s(n,j)*sn 1155 s(n,j) = s(n,i)*sn + s(n,j)*cs 1156 s(n,i) = t 1157 end do ! n 1158 end do ! i 1159 do i = 1,nss,ndf 1160 j = i + 1 1161 do n = 1,nss 1162 t = s(i,n)*cs - s(j,n)*sn 1163 s(j,n) = s(i,n)*sn + s(j,n)*cs 1164 s(i,n) = t 1165 end do ! n 1166 end do ! i 1167 1168 else 1169 1170 do i=1,nss,ndf 1171 j = i + 1 1172 t = s(i,1)*cs + s(j,1)*sn 1173 s(j,1) = -s(i,1)*sn + s(j,1)*cs 1174 s(i,1) = t 1175 end do ! i 1176 1177 endif 1178 1179 end 1180 1181 subroutine frcn2d(ix,p,ndf,dt,st,numnp) 1182 1183 implicit none 1184 1185 integer ndf,numnp 1186 integer i,j,ll 1187 1188 integer ix(*) 1189 real*8 dt(numnp),st(numnp,*),p(ndf,*) 1190 1191 save 1192 1193 do i = 1,2 1194 1195 ll = ix(i) 1196 if(ll.gt.0) then 1197 1198 dt(ll) = dt(ll) + 1.d0 1199 1200! Stress projections 1201 1202 do j = 1,3 1203 st(ll,j) = st(ll,j) + p(j,i) 1204 end do 1205 endif 1206 end do 1207 1208 end 1209 1210 subroutine frans2d(d,ul,xl,s,p,ndf,ndm,nst,isw) 1211 1212!-----[--+---------+---------+---------+---------+---------+---------+-] 1213! Two dimensional Euler-Bernoulli frame element 1214! N.B. ELASTIC ONLY 1215 1216! Control Information: 1217 1218! ndm - Spatial dimension of problem = 2 1219! ndf - Number degree-of-freedoms at node >= 3 1220! ( 1 = u_1 ; 2 = u_2 ; 3 = theta ) 1221! nen - Two node element (nel = 2) >= 2 1222!-----[--+---------+---------+---------+---------+---------+---------+-] 1223 implicit none 1224 1225 include 'bdata.h' 1226 include 'cdata.h' 1227 include 'eldata.h' 1228 include 'eltran.h' 1229 include 'iofile.h' 1230 include 'prld1.h' 1231 1232 1233 integer i,j,ii,jj,i1,j1,i2,j2,ndf,ndm,nst,isw 1234 real*8 cs,sn,le,xx,yy,xn,xm,vv,eps,chi,gam,EA,EI,RA 1235 real*8 b1,b2 1236 1237 real*8 d(*),ul(ndf,nen,*),xl(ndm,*),s(nst,*),p(*) 1238 real*8 sm(6,6),pm(6) 1239 1240 save 1241 1242! COMPUTE ELEMENT ARRAYS 1243 1244 cs = xl(1,2) - xl(1,1) 1245 sn = xl(2,2) - xl(2,1) 1246 le = sqrt(cs*cs+sn*sn) 1247 cs = cs/le 1248 sn = sn/le 1249 1250! Compute elment stiffness/residual arrays 1251 1252 if(isw.eq.3 .or. isw.eq.6) then 1253 1254! Set body loading factors 1255 1256 if(int(d(81)).gt.0) then 1257 b1 = 0.5*le*(d(11) + prldv(int(d(81)))*d(71)) 1258 else 1259 b1 = 0.5*le*d(11)*dm 1260 endif 1261 if(int(d(82)).gt.0) then 1262 b2 = 0.5*le*(d(12) + prldv(int(d(82)))*d(72)) 1263 else 1264 b2 = 0.5*le*d(12)*dm 1265 endif 1266 1267! Compute momentum residual and tangent matrix 1268 1269 call pzero(sm,36) 1270 EA = d(21)*d(32) 1271 EI = d(21)*d(33) 1272 RA = d(4)*d(32) 1273 call beam2d(s ,EA,EI,le,cs,sn,nst,ndf) 1274 call massf2(sm,pm,d(7),RA,le,cs,sn,6,3) 1275 1276! Stress and mass modification to residual and stiffness 1277 1278 i1 = 0 1279 i2 = 0 1280 do ii = 1,2 1281 p(1+i1) = p(1+i1) + b1 1282 p(2+i1) = p(2+i1) + b2 1283 do i = 1,ndf 1284 j1 = 0 1285 j2 = 0 1286 do jj = 1,2 1287 do j = 1,ndf 1288 1289! Residual 1290 1291 p(i+i1) = p(i+i1) - s(i+i1,j+j1)*ul(j,jj,1) 1292 & - sm(i+i2,j+j2)*ul(j,jj,5) 1293 1294! Tangent 1295 1296 s(i+i1,j+j1) = s(i+i1,j+j1)*ctan(1) 1297 & + sm(i+i2,j+j2)*ctan(3) 1298 1299 end do 1300 j1 = j1 + ndf 1301 j2 = j2 + 3 1302 end do 1303 1304! Lumped mass modification to residual and tangent 1305 1306 p(i+i1) = p(i+i1) - pm(i+i2)*ul(i,ii,5) 1307 1308 s(i+i1,i+i1) = s(i+i1,i+i1) + pm(i+i2)*ctan(3) 1309 1310 end do 1311 i1 = i1 + ndf 1312 i2 = i2 + 3 1313 end do 1314 1315! Compute element output quantities 1316 1317 elseif(isw.eq.4) then 1318 1319 xx = 0.5d0*(xl(1,1)+xl(1,2)) 1320 yy = 0.5d0*(xl(2,1)+xl(2,2)) 1321 eps = (cs*(ul(1,2,1)-ul(1,1,1)) 1322 & + sn*(ul(2,2,1)-ul(2,1,1)))/le 1323 chi = (ul(3,2,1)-ul(3,1,1))/le 1324 gam =-12.d0*(-sn*(ul(1,1,1)-ul(1,2,1)) 1325 & + cs*(ul(2,1,1)-ul(2,2,1)))/le**3 1326 & - 6.d0*(ul(3,1,1)+ul(3,2,1))/le/le 1327 xn = d(21)*d(32)*eps 1328 xm = d(21)*d(33)*chi 1329 vv = d(21)*d(33)*gam 1330 1331 mct = mct - 1 1332 if(mct.le.0) then 1333 write(iow,2000) o,head 1334 if(ior.lt.0) write(*,2000) o,head 1335 mct = 50 1336 endif 1337 1338 write(iow,2001) n,ma,xx,yy,xn,xm,vv,eps,chi,gam 1339 if(ior.lt.0) write(*,2001) n,ma,xx,yy,xn,xm,vv,eps,chi,gam 1340 1341! Compute element mass arrays 1342 1343 elseif(isw.eq.5) then 1344 RA = d(4)*d(32) 1345 call massf2(s,p,d(7),RA,le,cs,sn,nst,ndf) 1346 end if 1347 1348! Formats 1349 13502000 format(a1,20a4//5x,'Beam Element Stresses'// 1351 & ' Elmt Matl x-coord y-coord ', 1352 & ' Force Moment Shear'/ 1353 & 43x,'Strain Curvature Gamma') 1354 13552001 format(i8,i6,1p,5e12.4/38x,1p,3e12.4) 1356 1357 end 1358 1359 subroutine beam2d(s,ea,ei,le,cs,sn,nst,ndf) 1360 1361! Stiffness for frame element 1362 1363 implicit none 1364 1365 integer nst,ndf,i,j,k 1366 real*8 ea,ei,le,cs,sn,t 1367 1368 real*8 s(nst,nst) 1369 1370 i = ndf + 1 1371 j = ndf + 2 1372 k = ndf + 3 1373 1374 t = ea/le 1375 s(1,1) = t 1376 s(i,i) = t 1377 s(1,i) =-t 1378 s(i,1) =-t 1379 1380 t = 12.d0*ei/le**3 1381 s(2,2) = t 1382 s(j,j) = t 1383 s(2,j) =-t 1384 s(j,2) =-t 1385 t = (ei+ei)/le 1386 s(3,3) = t + t 1387 s(k,k) = t + t 1388 s(3,k) = t 1389 s(k,3) = t 1390 t = 6.d0*ei/le**2 1391 s(2,3) = t 1392 s(3,2) = t 1393 s(2,k) = t 1394 s(k,2) = t 1395 s(3,j) =-t 1396 s(j,3) =-t 1397 s(j,k) =-t 1398 s(k,j) =-t 1399 1400 call rotaf2(s,cs,sn,nst,ndf) 1401 1402 end 1403 1404 subroutine massf2(s,p,cfac,ra,le,cs,sn,nst,ndf) 1405 1406! Frame mass matrix 1407 1408 implicit none 1409 1410 integer nst,ndf,i,j,l,ii(4) 1411 real*8 cfac,lfac,ra,le,cs,sn,t,dv,s1,s2,s3 1412 real*8 p(nst),s(nst,nst),sg(4),ag(4),bb(4) 1413 1414 data ii /2,3,5,6/ 1415 1416 ii(3) = ndf + 2 1417 ii(4) = ndf + 3 1418 1419! Lumped mass matrix 1420 1421 t = 0.5d0*ra*le 1422 p(1) = t 1423 p(2) = t 1424 p(ndf+1) = t 1425 p(ndf+2) = t 1426 1427! Consistent mass matrix 1428 1429 s(1 ,1 ) = 0.6666666666666667d0*t 1430 s(1 ,ndf+1) = 0.3333333333333333d0*t 1431 s(ndf+1,ndf+1) = s(1,1) 1432 1433 call int1d(4,sg,ag) 1434 1435 do l = 1,4 1436 dv = t*ag(l) 1437 s1 = 0.5d0 + 0.5d0*sg(l) 1438 s2 = s1*s1 1439 s3 = s1*s2 1440 bb(3) = 3.d0*s2 - s3 - s3 1441 bb(4) = le*(s3 - s2) 1442 bb(1) = 1.d0 - bb(3) 1443 bb(2) = le*(s1 - s2) + bb(4) 1444 do i = 1,4 1445 s1 = bb(i)*dv 1446 do j = i,4 1447 s(ii(i),ii(j)) = s(ii(i),ii(j)) + s1*bb(j) 1448 end do 1449 end do 1450 end do 1451 1452 do i = 1,6 1453 do j = 1,i 1454 s(i,j) = s(j,i) 1455 end do 1456 end do 1457 1458 lfac = 1.d0 - cfac 1459 do i = 1,6 1460 do j = 1,6 1461 s(i,j) = cfac*s(i,j) 1462 end do 1463 s(i,i) = s(i,i) + lfac*p(i) 1464 end do 1465 1466 call rotaf2(s,cs,sn,nst,ndf) 1467 1468 end 1469 1470 subroutine rotaf2(s,cs,sn,nst,ndf) 1471 1472! Rotate arrays from local to global dof's 1473 1474 implicit none 1475 1476 integer nst,ndf,i,j,n 1477 real*8 cs,sn,t 1478 real*8 s(nst,nst) 1479 1480! Check angle (perform rotation if necessary) 1481 1482 if(cs.lt.0.99999999d0) then 1483 1484 do i = 1,nst,ndf 1485 j = i + 1 1486 do n = 1,nst 1487 t = s(n,i)*cs - s(n,j)*sn 1488 s(n,j) = s(n,i)*sn + s(n,j)*cs 1489 s(n,i) = t 1490 end do 1491 end do 1492 do i = 1,nst,ndf 1493 j = i + 1 1494 do n = 1,nst 1495 t = cs*s(i,n) - sn*s(j,n) 1496 s(j,n) = sn*s(i,n) + cs*s(j,n) 1497 s(i,n) = t 1498 end do 1499 end do 1500 1501 end if 1502 1503 end 1504 1505 subroutine franf2d(d,ul,xl,ix,s,r,ndf,ndm,nst,isw) 1506 1507!-----[--+---------+---------+---------+---------+---------+---------+-] 1508! Purpose: Two dimensional Euler-Bernoulli frame element: Second 1509! order theory. Includes one enhanced mode for axial and 1510! bending strain match. 1511!-----[--.----+----.----+----.-----------------------------------------] 1512 implicit none 1513 1514 include 'bdata.h' 1515 include 'bm2com.h' 1516 include 'cdata.h' 1517 include 'cdat1.h' 1518 include 'eldata.h' 1519 include 'eltran.h' 1520 include 'hdata.h' 1521 include 'iofile.h' 1522 include 'pmod2d.h' 1523 include 'prld1.h' 1524 include 'prstrs.h' 1525 include 'ptdat6.h' 1526 include 'tdata.h' 1527 1528 include 'comblk.h' 1529 1530 logical noconv 1531 integer ndf,ndm,nst,isw 1532 integer i,ii,itmax, j,jj, ll,lint, mm,nn 1533 real*8 cs,sn,b1,b2,dva,dvi,xjac, energy 1534 real*8 len, hle 1535 real*8 rhoa, rhoi, ctan1, ctan3 1536 1537 integer ix(*) 1538 1539 real*8 aa(4,4,2),shpw(4,2,4),shpt(4,2,4),shpu(2,2,4),dx(4) 1540 real*8 sg(4),wg(4), dudx(4), dwdx(4), eps(4), kap(4), gam(4) 1541 real*8 bmat(2,3,2),baii(4,2),forc(4,5),xx(2),b(2),nxi(3) 1542 real*8 xl(ndm,*),ul(ndf,nen,*) 1543 real*8 d(*),r(ndf,*),s(nst,nst), gen(3,2) 1544 real*8 duen, uen, ben, hen, EA, tol 1545 1546 save 1547 1548 data itmax / 20 / 1549 data tol / 1.d-10 / 1550 1551! Second order (nonlinear) deformation BEAM element. 1552 1553! d(1)*d(32) = EA 1554! d(1)*d(33) = EI 1555! d(4)*d(32) = rho*A 1556! d(4)*d(33) = rho*I 1557 1558! Compute element length and direction cosines 1559 1560 if(isw.ge.2) then 1561 cs = xl(1,nel) - xl(1,1) 1562 sn = xl(2,nel) - xl(2,1) 1563 len = sqrt(cs*cs + sn*sn) 1564 cs = cs/len 1565 sn = sn/len 1566 hle = 0.5d0*len 1567 1568 rhoa = d(4)*d(32) 1569 rhoi = d(4)*d(33) 1570 endif 1571 1572! Read data 1573 1574 if(isw.eq.1) then 1575 1576! Increment history storage if necessary 1577 1578 nh1 = nh1 + 1 ! Enhanced strain parameter 1579 1580! Compute mass array 1581 1582 elseif(isw.eq.5) then 1583 1584 lint = nel + 1 1585 call int1d(lint,sg,wg) 1586 1587! Compute lumped mass matrix 1588 1589 call bm2trn (xl,cs,sn,ndm*nel,ndm,2) 1590 do ll = 1,lint 1591 1592! Compute shape functions 1593 1594 call shp1db(sg(ll),xl,shpu,ndm,nel,xjac) 1595 1596 dva = wg(ll)*xjac*rhoa 1597 dvi = wg(ll)*xjac*rhoi*d(69) 1598 1599! For each node j compute db = rho*shape*dv 1600 1601 do j = 1,nel 1602 b1 = shpu(2,j,ll)*dva 1603 b2 = shpu(2,j,ll)*dvi 1604 1605! Compute a lumped mass 1606 1607 r(1,j) = r(1,j) + b1 1608 r(2,j) = r(2,j) + b1 1609 r(3,j) = r(3,j) + b2 1610 end do ! j 1611 end do ! ll 1612 1613! Place in consistent mass 1614 1615 jj = 0 1616 do j = 1,nel 1617 s(jj+1,jj+1) = r(1,j) 1618 s(jj+2,jj+2) = r(2,j) 1619 s(jj+3,jj+3) = r(3,j) 1620 jj = jj + ndf 1621 end do ! j 1622 1623 else 1624 1625! Quadrature terms 1626 1627 lint = nel + 1 1628 call int1d(lint,sg,wg) 1629 1630! Transform to local coordinates 1631 1632 call bm2trn (ul(1,1,1),cs,sn,ndf*nel,ndf,2) 1633 call bm2trn (ul(1,1,2),cs,sn,ndf*nel,ndf,2) 1634 call bm2trn (ul(1,1,4),cs,sn,ndf*nel,ndf,2) 1635 call bm2trn (ul(1,1,5),cs,sn,ndf*nel,ndf,2) 1636 call bm2trn (xl ,cs,sn,ndm*nel,ndm,2) 1637 dva = rhoa*hle 1638 dvi = rhoi*hle*d(69) 1639 1640 do ll = 1,lint 1641 1642! Shape functions 1643 1644 call shp1dh(sg(ll),len,shpw(1,1,ll),shpt(1,1,ll)) 1645 shpu(1,1,ll) = -1.d0/len 1646 shpu(1,2,ll) = -shpu(1,1,ll) 1647 shpu(2,1,ll) = 0.5d0 - 0.5d0*sg(ll) 1648 shpu(2,2,ll) = 0.5d0 + 0.5d0*sg(ll) 1649 1650 dx(ll) = wg(ll)*hle 1651 1652! Form displacement derivatives from nodal displacements 1653 1654 dwdx(ll) = shpw(1,1,ll)*ul(2,1,1) + shpw(1,2,ll)*ul(2,2,1) 1655 & + shpt(1,1,ll)*ul(3,1,1) + shpt(1,2,ll)*ul(3,2,1) 1656 1657 dudx(ll) = 0.5d0*(ul(1,2,1) - ul(1,1,1))/hle 1658 1659 eps(ll) = dudx(ll) + 0.5d0*(dudx(ll)*dudx(ll) 1660 & + dwdx(ll)*dwdx(ll)) 1661 1662 kap(ll) = shpw(2,1,ll)*ul(2,1,1) + shpw(2,2,ll)*ul(2,2,1) 1663 & + shpt(2,1,ll)*ul(3,1,1) + shpt(2,2,ll)*ul(3,2,1) 1664 1665 gam(ll) = shpw(3,1,ll)*ul(2,1,1) + shpw(3,2,ll)*ul(2,2,1) 1666 & + shpt(3,1,ll)*ul(3,1,1) + shpt(3,2,ll)*ul(3,2,1) 1667 end do ! ll 1668 1669! Enhanced strain computation 1670 1671 if(etype.eq.3) then 1672 uen = hr(nh1) 1673 ii = 0 1674 noconv = .true. 1675 do while(noconv) 1676 1677 ii = ii + 1 1678 1679! Zero enhanced terms 1680 1681 ben = 0.0d0 1682 hen = 0.0d0 1683 1684 do ll = 1,lint 1685 1686 EA = d(1)*d(32) 1687 hen = hen + sg(ll)*EA*sg(ll)*dx(ll)*ctan(1) 1688 ben = ben - sg(ll)*EA*(eps(ll) + sg(ll)*uen)*dx(ll) 1689 1690 end do ! ll 1691 1692 hen = 1.d0/ hen 1693 duen = ben * hen 1694 uen = uen + duen 1695 1696 if(abs(duen).le.tol*abs(uen) .or. ben.eq.0.0d0) then 1697 noconv = .false. 1698 elseif(ii.gt.itmax) then 1699 noconv = .false. 1700 write(*,*) 'WARNING -- No convergence in FRANF2D' 1701 write(*,*) duen,uen,hen 1702 endif 1703 1704 end do ! while 1705 1706! Save enhance mode parameter 1707 1708 hr(nh2) = uen 1709 1710 else 1711 1712 hen = 0.0d0 1713 1714 end if ! etype 1715 1716! Stiffness and residual computation 1717 1718 if(isw.eq.3 .or. isw.eq.6) then 1719 1720! Zero enhanced coupling array 1721 1722 do i = 1,3 1723 gen(i,1) = 0.0d0 1724 gen(i,2) = 0.0d0 1725 end do ! i 1726 1727! Final tangent form 1728 1729 do ll = 1,lint 1730 1731 forc(1,ll) = d(1)*d(32)*(eps(ll) + sg(ll)*uen)*dx(ll) 1732 forc(2,ll) = d(1)*d(33)*kap(ll)*dx(ll) 1733 aa(1,1,1) = d(1)*d(32) 1734 aa(2,2,1) = d(1)*d(33) 1735 aa(1,2,1) = 0.0d0 1736 aa(2,1,1) = 0.0d0 1737 1738 ctan1 = ctan(1)*dx(ll) 1739 do jj = 1,4 1740 do ii = 1,4 1741 aa(ii,jj,1) = aa(ii,jj,1)*ctan1 1742 end do ! ii 1743 end do ! jj 1744 1745! Compute strain-displacement matrices for two nodes 1746 1747 do i = 1,2 1748 bmat(1,1,i) = shpu(1,i,ll)*(1.d0 + dudx(ll)) 1749 bmat(2,1,i) = 0.0d0 1750 bmat(1,2,i) = shpw(1,i,ll)*dwdx(ll) 1751 bmat(2,2,i) = shpw(2,i,ll) 1752 bmat(1,3,i) = shpt(1,i,ll)*dwdx(ll) 1753 bmat(2,3,i) = shpt(2,i,ll) 1754 end do ! i 1755 1756! Mechanical tangent terms 1757 1758 mm = 0 1759 do ii = 1,nel 1760 1761 do i = 1,3 1762 1763! B^T * AA 1764 1765 baii(i,1) = (bmat(1,i,ii)*aa(1,1,1) + 1766 & bmat(2,i,ii)*aa(2,1,1)) 1767 baii(i,2) = (bmat(1,i,ii)*aa(1,2,1) + 1768 & bmat(2,i,ii)*aa(2,2,1)) 1769 1770! Residual 1771 1772 r(i,ii) = r(i,ii) - bmat(1,i,ii)*forc(1,ll) 1773 & - bmat(2,i,ii)*forc(2,ll) 1774 1775! Enhanced stiffness 1776 1777 gen(i,ii) = gen(i,ii) + baii(i,1)*sg(ll) 1778 1779 end do ! i 1780 1781! Tangent 1782 1783 if(isw.eq.3) then 1784 1785 nxi(1) = shpu(1,ii,ll)*forc(1,ll)*ctan(1) 1786 nxi(2) = shpw(1,ii,ll)*forc(1,ll)*ctan(1) 1787 nxi(3) = shpt(1,ii,ll)*forc(1,ll)*ctan(1) 1788 nn = 0 1789 do jj = 1,nel 1790 1791! Material part 1792 1793 do j = 1,3 1794 do i = 1,3 1795 s(mm+i,nn+j) = s(mm+i,nn+j) 1796 & + baii(i,1)*bmat(1,j,jj) 1797 & + baii(i,2)*bmat(2,j,jj) 1798 end do ! i 1799 end do ! j 1800 1801! Geometric part 1802 1803 s(mm+1,nn+1) = s(mm+1,nn+1) + nxi(1)* shpu(1,jj,ll) 1804 s(mm+2,nn+2) = s(mm+2,nn+2) + nxi(2)*shpw(1,jj,ll) 1805 s(mm+2,nn+3) = s(mm+2,nn+3) + nxi(2)*shpt(1,jj,ll) 1806 s(mm+3,nn+2) = s(mm+3,nn+2) + nxi(3)*shpw(1,jj,ll) 1807 s(mm+3,nn+3) = s(mm+3,nn+3) + nxi(3)*shpt(1,jj,ll) 1808 1809 nn = nn + ndf 1810 end do ! jj 1811 endif 1812 mm = mm + ndf 1813 end do ! ii 1814 end do ! ll 1815 1816! Lumped inertia contributions 1817 1818 ctan3 = ctan(3) + d(77)*ctan(2) 1819 jj = 0 1820 do ii = 1,nel 1821 r(1,ii) = r(1,ii) - dva*(ul(1,ii,5)+d(77)*ul(1,ii,4)) 1822 r(2,ii) = r(2,ii) - dva*(ul(2,ii,5)+d(77)*ul(2,ii,4)) 1823 r(3,ii) = r(3,ii) - dvi*(ul(3,ii,5)+d(77)*ul(3,ii,4)) 1824 s(jj+1,jj+1) = s(jj+1,jj+1) + dva*ctan3 1825 s(jj+2,jj+2) = s(jj+2,jj+2) + dva*ctan3 1826 s(jj+3,jj+3) = s(jj+3,jj+3) + dvi*ctan3 1827 jj = jj + ndf 1828 end do 1829 1830! Transform stiffness and residual to global coordinates 1831 1832 if(isw.eq.3) then 1833 1834! Static condensation 1835 1836 nn = 0 1837 do jj = 1,nel 1838 do j = 1,3 1839 duen = gen(j,jj)*hen 1840 mm = 0 1841 do ii = 1,nel 1842 do i = 1,3 1843 s(i+mm,j+nn) = s(i+mm,j+nn) - gen(i,ii)*duen 1844 end do ! i 1845 mm = mm + ndf 1846 end do ! ii 1847 end do ! j 1848 nn = nn + ndf 1849 end do ! jj 1850 1851 call bm2trn (s,cs,sn,nst,ndf,1) 1852 endif 1853 call bm2trn ( r,cs,-sn,nst,ndf,2) 1854 1855! Set body loading factors 1856 1857 if(int(d(74)).gt.0) then 1858 b(1) = hle*(d(11) + prldv(int(d(74)))*d(71)) 1859 else 1860 b(1) = hle*d(11)*dm 1861 endif 1862 if(int(d(75)).gt.0) then 1863 b(2) = hle*(d(12) + prldv(int(d(75)))*d(72)) 1864 else 1865 b(2) = hle*d(12)*dm 1866 endif 1867 1868! Add body force components 1869 1870 r(1,1) = r(1,1) + b(1) 1871 r(1,2) = r(1,2) + b(1) 1872 r(2,1) = r(2,1) + b(2) 1873 r(2,2) = r(2,2) + b(2) 1874 1875! Output forces 1876 1877 elseif(isw.eq.4 .or. isw.eq.8) then 1878 1879! Member forces 1880 1881 do ll = 1,lint 1882 1883 defa(1,1) = eps(ll) + sg(ll)*uen 1884 defa(2,1) = gam(ll) 1885 defa(3,1) = kap(ll) 1886 1887 forc(1,ll) = d(1)*d(32)*defa(1,1) 1888 forc(2,ll) = d(1)*d(33)*defa(2,1) 1889 forc(3,ll) = d(1)*d(33)*defa(3,1) 1890 1891! Stress output 1892 1893 if(isw.eq.4) then 1894 1895 do i = 1,ndm 1896 xx(i) = 0. 1897 do ii = 1,nel 1898 xx(i) = xx(i) + xl(i,ii)*shpu(2,ii,ll) 1899 end do ! ii 1900 end do ! i 1901 mct = mct - 3 1902 if (mct.le.0) then 1903 write(iow,2001) o,head,ttim 1904 if(ior.lt.0) write(*,2001) o,head,ttim 1905 mct = 50 1906 endif 1907 write(iow,2002) n,ma,(xx(i),i=1,2), 1908 & (forc(i,ll),i=1,3),(defa(i,1),i=1,3) 1909 if(ior.lt.0) then 1910 write(*,2002) n,ma,(xx(i),i=1,2), 1911 & (forc(i,ll),i=1,3),(defa(i,1),i=1,3) 1912 endif 1913 1914! Stress projections save 1915 1916 else 1917 1918! Compute strain-displacement matrices for two nodes 1919 1920 do i = 1,2 1921 bmat(1,1,i) = shpu(1,i,ll)*(1.d0 + dudx(ll)) 1922 bmat(2,1,i) = 0.0d0 1923 bmat(1,2,i) = shpw(1,i,ll)*dwdx(ll) 1924 bmat(2,2,i) = shpw(2,i,ll) 1925 bmat(1,3,i) = shpt(1,i,ll)*dwdx(ll) 1926 bmat(2,3,i) = shpt(2,i,ll) 1927 end do ! i 1928 1929! End forces 1930 1931 do ii = 1,nel 1932 1933 do i = 1,3 1934 r(i,ii) = r(i,ii) - (bmat(1,i,ii)*forc(1,ll) 1935 & + bmat(2,i,ii)*forc(3,ll))*dx(ll) 1936 end do ! i 1937 1938 end do ! ii 1939 end if 1940 end do ! ll 1941 1942! Projection on end notes (uses reactions) 1943 1944 if(isw.eq.8) then 1945 do i = 1,3 1946 r(i,2) = -r(i,2) 1947 end do 1948 call frcn2d(ix,r,ndf,hr(nph),hr(nph+numnp),numnp) 1949 endif 1950 1951! Compute energy 1952 1953 elseif(isw.eq.13) then 1954 1955 dva = hle*rhoa 1956 dvi =hle*rhoi*d(69) 1957 1958! Compute internal energy 1959 1960 do ll = 1,lint 1961 1962! Compute energy density from stress and deformation 1963 1964 call shp1db(sg(ll),xl,shpu,ndm,nel,xjac) 1965 dx(ll) = wg(ll)*xjac 1966 call b2mods (d,ul,forc,aa,energy,shpu,ndf,nen,isw) 1967 1968! Accumulate energy 1969 1970 epl(8) = epl(8) + 0.5d0*energy*dx(ll) 1971 1972 end do ! ll 1973 1974! Compute kinetic energy for lumped mass 1975 1976 epl(7) = epl(7) + 0.5d0*dva*(ul(1,1,4)**2 + ul(1,2,4)**2 1977 & + ul(2,1,4)**2 + ul(2,2,4)**2) 1978 & + 0.5d0*dvi*(ul(3,1,4)**2 + ul(3,2,4)**2) 1979 1980 endif 1981 1982 endif 1983 1984! Formats 1985 19862001 format(a1,20a4/5x,'time',e13.5,5x,' element forces '// 1987 & 43x,'********* FORCE / STRAIN *********'/ 1988 & 3x,'element material', 1989 & 3x,'1-coord',3x,'2-coord',6x,'n-dir',8x,'s-dir',8x,'m-dir'/) 1990 19912002 format(2i10,0p,2f10.3,1p,3e13.4/40x,1p,3e13.4) 1992 1993 end 1994