1!$Id:$ 2 subroutine fld2d1(d,ul,xl,ix,s,r,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! Purpose: 2-D Finite Deformation Elasticity Routine 11! Remark: This is a standard displacement model 12 13! Inputs: 14! d(*) - Material set parameters 15! ul(ndf,*) - Nodal solution parameters for element 16! xl(ndm,*) - Nodal coordinates for element 17! ix(*) - Element nodal connection list 18! ndf - Number dof/node 19! ndm - Spatial dimension of mesh 20! nst - Dimension of element arrays 21! isw - Switch to control action 22 23! Outputs: 24! s(nst,*) - Element matrix 25! p(nst) - Element vector 26!-----[--.----+----.----+----.-----------------------------------------] 27 28 implicit none 29 30 include 'debugs.h' 31 include 'bdata.h' 32 include 'cdata.h' 33 include 'eldata.h' 34 include 'elengy.h' 35 include 'elplot.h' 36 include 'eltran.h' 37 include 'hdata.h' 38 include 'iofile.h' 39 include 'pmod2d.h' 40 include 'prstrs.h' 41 include 'ptdat6.h' 42 include 'qudshp.h' 43 include 'comblk.h' 44 45 integer ndf,ndm,nst,isw 46 integer i,i1, j,jj,j1, l, nn,nhv, istrt 47 48 real*8 bdb,bd3,dl,dc,di,dvol0,dvol,dmas0, jac0 49 real*8 cfac,lfac,xx1,xx2,xx3,yy1, tempi, ta 50 51 integer ix(*) 52 real*8 d(*),ul(ndf,nen,*),xl(ndm,*),s(nst,*) 53 real*8 r(ndf,*),r1(3,64),al(2),ac(2),vl(2) 54 real*8 bbd(4,2),bb(6),shpr(64) 55 real*8 df(3,3,64),f(9,2,64),detf(2,64) 56 real*8 ds(6,6,5),dd(6,6),sigv(9),sigl(9,64) 57 58 save 59 60 data lint / 0 / 61 data f /1152*0.0d0 / 62 63! TEMPORARY SET OF TEMPERATURE VALUE 64 65 data ta / 0.0d0/ 66 67! No action for isw = 1 68 69 if(isw.eq.1) return 70 71! Set quadrature and go to process isw 72 73 call quadr2d(d,.true.) 74 75 if(debug) call mprint(sg2,3,lint,3,'SG2') 76 77! COMPUTE TANGENT STIFFNESS AND RESIDUAL FORCE VECTOR 78 79! Compute shape functions and derivatives in reference configuration 80 81 do l = 1,lint 82 call interp2d(l, xl,ix,ndm,nel, .false.) 83 84 if(debug) call mprint(shp2(1,1,l),3,nel,3,'SHP') 85 end do 86 if(debug) call mprint(jac,1,lint,1,'JAC') 87 88! Compute deformation gradient and determinant; transform shape 89! functions to current configuration. 90 91 do l = 1,lint 92 call kine2d(shp2(1,1,l),xl,ul,f(1,1,l),df(1,1,l), 93 & detf(1,l),ndm,ndf,nel,nen) 94 end do 95 96! Integration order set to static 97 98 if(d(7).lt.0.0d0) then 99 cfac = 0.0d0 100 lfac = 0.0d0 101 else 102 cfac = d(7) 103 lfac = 1.0d0 - cfac 104 endif 105 106 nhv = nint(d(15)) 107 istrt = nint(d(84)) 108 109 if(mod(isw,4).eq.0) go to 4 110 111! LOOP OVER GAUSS POINTS 112 113 nn = 0 114 do l = 1,lint 115 116! Set reference coordinates 117 118 xx1 = 0.0d0 119 xx2 = 0.0d0 120 do i = 1,nel 121 xx1 = xx1 + shp2(3,i,l)*xl(1,i) 122 xx2 = xx2 + shp2(3,i,l)*xl(2,i) 123 end do 124 125! Check for axisymmetry 126 127 if(stype.eq.3) then 128 jac0 = jac(l) 129 dvol0 = jac0*xx1 130 do i = 1,nel 131 shpr(i) = shp2(3,i,l)/xx1 132 end do ! i 133 else 134 jac0 = 0.0d0 135 dvol0 = jac(l) 136 do i = 1,nel 137 shpr(i) = 0.0d0 138 end do ! i 139 end if 140 dvol = dvol0*detf(1,l) 141 dmas0 = dvol0*d(4) 142 143! Compute Cauchy stresses and spatial tangent tensor 144 145 call modlfd(l,d,f(1,1,l),df(1,1,l),detf(1,l),ta, 146 & hr(nn+nh1),hr(nn+nh2),nhv,istrt,ds,sigv,bb, 147 & .false.,isw) 148 149 if(isw.eq.13) then 150 151 epl(8) = epl(8) + estore*dvol0 152 153! Compute velocity at point 154 155 vl(1) = shp2(3,1,l)*ul(1,1,4) + shp2(3,2,l)*ul(1,2,4) 156 & + shp2(3,3,l)*ul(1,3,4) + shp2(3,4,l)*ul(1,4,4) 157 158 vl(2) = shp2(3,1,l)*ul(2,1,4) + shp2(3,2,l)*ul(2,2,4) 159 & + shp2(3,3,l)*ul(2,3,4) + shp2(3,4,l)*ul(2,4,4) 160 161 tempi = 0.0d0 162 do i = 1,nel 163 tempi = tempi 164 & + (ul(1,i,4)**2 + ul(2,i,4)**2)*shp2(3,i,l) 165 end do ! i 166 167! Accumulate kinetic energy 168 169 epl(7) = epl(7) + 0.5d0*(lfac*tempi 170 & + cfac*(vl(1)**2 + vl(2)**2))*dmas0 171 172 elseif(isw.ne.14) then 173 174! Store stress values for tplot 175 176 j1 = 6*(l-1) 177 do j = 1,4 178 tt(j+j1) = sigv(j) 179 end do 180 181! Multiply tangent moduli and stresses by volume element. 182 183 do i = 1,4 184 sigv(i) = sigv(i)*dvol 185 do j = 1,4 186 dd(i,j) = ds(i,j,1)*dvol*ctan(1) 187 end do 188 end do 189 190! Compute accelerations 191 192 al(1) = cfac*(shp2(3,1,l)*ul(1,1,5) + shp2(3,2,l)*ul(1,2,5) 193 & + shp2(3,3,l)*ul(1,3,5) + shp2(3,4,l)*ul(1,4,5)) 194 195 al(2) = cfac*(shp2(3,1,l)*ul(2,1,5) + shp2(3,2,l)*ul(2,2,5) 196 & + shp2(3,3,l)*ul(2,3,5) + shp2(3,4,l)*ul(2,4,5)) 197 198! COMPUTE STRESS DIVERGENCE AND INERTIA TERMS 199 200 xx1 = xx1*d(65)**2 201 xx2 = xx2*d(65)**2 202 203 do i = 1,nel 204 205! Compute inertial and body load effects 206 207 ac(1) = (al(1) + lfac*ul(1,i,5))*dmas0 208 & - d(11)*dvol0 - xx1*dmas0 209 ac(2) = (al(2) + lfac*ul(2,i,5))*dmas0 210 & - d(12)*dvol0 - xx2*dmas0 211 212! Stress divergence term (used in geometric stiffness) 213 214 r1(1,i) = shp2(1,i,l)*sigv(1) + shp2(2,i,l)*sigv(4) 215 r1(2,i) = shp2(1,i,l)*sigv(4) + shp2(2,i,l)*sigv(2) 216 217! Element residual 218 219 r(1,i) = r(1,i) - r1(1,i) - ac(1)*shp2(3,i,l) 220 & - shpr(i)*sigv(3) 221 r(2,i) = r(2,i) - r1(2,i) - ac(2)*shp2(3,i,l) 222 223 end do 224 225! COMPUTE K (s(nst,nst) = K) 226 227 if(isw.eq.3) then 228 229! PART 1. - Geometric and inertial part. 230 231 dc = cfac*ctan(3)*dmas0 232 dl = lfac*ctan(3)*dmas0 233 i1 = 0 234 do i = 1,nel 235 236 do jj = 1,2 237 s(i1+jj,i1+jj) = s(i1+jj,i1+jj) + shp2(3,i,l)*dl 238 end do 239 240 di = dc*shp2(3,i,l) 241 242! Include geometric stiffness 243 244 if(gflag) then 245 bd3 = shpr(i)*sigv(3)*ctan(1) 246 j1 = 0 247 do j = 1,i 248 bdb = (r1(1,i)*shp2(1,j,l) 249 & + r1(2,i)*shp2(2,j,l))*ctan(1) 250 & + di*shp2(3,j,l) 251 s(i1+1,j1+1) = s(i1+1,j1+1) + bdb + bd3*shpr(j) 252 s(i1+2,j1+2) = s(i1+2,j1+2) + bdb 253 j1 = j1 + ndf 254 end do 255 256! Include inertia only 257 258 else 259 j1 = 0 260 do j = 1,i 261 bdb = di*shp2(3,j,l) 262 s(i1+1,j1+1) = s(i1+1,j1+1) + bdb 263 s(i1+2,j1+2) = s(i1+2,j1+2) + bdb 264 j1 = j1 + ndf 265 end do 266 endif 267 268 i1 = i1 + ndf 269 end do 270 271! PART 2. - Tangent modulus part (based upon dd-array) 272 273 i1 = 0 274 do i = 1,nel 275 276! Compute bmat-t * dd * dvol 277 278 do jj = 1,4 279 280 bbd(jj,1) = shp2(1,i,l)*dd(1,jj) 281 & + shpr( i )*dd(3,jj) 282 & + shp2(2,i,l)*dd(4,jj) 283 284 bbd(jj,2) = shp2(1,i,l)*dd(4,jj) 285 & + shp2(2,i,l)*dd(2,jj) 286 287 end do ! jj 288 289! Compute tangent stiffness 290 291 j1 = 0 292 do j = 1,i 293 294 s(i1+1,j1+1) = s(i1+1,j1+1) + bbd(1,1)*shp2(1,j,l) 295 & + bbd(3,1)*shpr( j ) 296 & + bbd(4,1)*shp2(2,j,l) 297 298 s(i1+2,j1+1) = s(i1+2,j1+1) + bbd(1,2)*shp2(1,j,l) 299 & + bbd(3,2)*shpr( j ) 300 & + bbd(4,2)*shp2(2,j,l) 301 302 s(i1+1,j1+2) = s(i1+1,j1+2) + bbd(4,1)*shp2(1,j,l) 303 & + bbd(2,1)*shp2(2,j,l) 304 305 s(i1+2,j1+2) = s(i1+2,j1+2) + bbd(4,2)*shp2(1,j,l) 306 & + bbd(2,2)*shp2(2,j,l) 307 308 j1 = j1 + ndf 309 end do 310 311 i1 = i1 + ndf 312 end do 313 314 endif ! end of tangent 315 316 endif ! end of isw options 317 318 nn = nn + nhv 319 320 end do 321 322 if(isw.eq.3) then 323 324! Compute upper part by symmetry 325 326 do j = 1,nst 327 do i = 1,j 328 s(i,j) = s(j,i) 329 end do 330 end do 331 332 if(debug) then 333 call mprint(r,ndf,nel,ndf,'Resid') 334 call mprint(s,nst,nst,nst,'Stiff') 335 endif 336 337 endif 338 339 return 340 341! OUTPUT STRESSES 342 343 4 xx1 = 0.d0 344 xx2 = 0.d0 345 xx3 = 0.d0 346 do i = 1,6 347 sigv(i) = 0.0d0 348 end do ! i 349 350! LOOP OVER GAUSS POINTS 351 352 nn = 0 353 354 do l = 1,lint 355 356! Compute Cauchy stresses and spatial tangent tensor at t-n+1 357 358 call modlfd(l,d,f(1,1,l),df(1,1,l),detf(1,l),ta, 359 & hr(nn+nh1),hr(nn+nh2),nhv,istrt,ds,sigl(1,l),bb, 360 & .false.,isw) 361 362 yy1 = 1.d0/dble(nel) 363 do i=1,nel 364 tempi = yy1 * shp2(3,i,l) 365 xx1 = xx1 + tempi*xl(1,i) 366 xx2 = xx2 + tempi*xl(2,i) 367 end do 368 369! Compute average stresses and jacobian for printing 370 371 do i = 1,4 372 sigv(i) = sigv(i) + yy1 * sigl(i,l) 373 end do 374 375 nn = nn + nhv 376 377 end do 378 379! OUTPUT STRESSES 380 381 if(isw.eq.4) then 382 mct = mct - 2 383 if(mct.le.0) then 384 write(iow,2001) o,head 385 if(ior.lt.0) write(*,2001) o,head 386 mct = 50 387 endif 388 389 write(iow,2002) n,ma,(sigv(jj),jj=1,6),xx1,xx2,xx3 390 if(ior.lt.0) then 391 write(*,2002) n,ma,(sigv(jj),jj=1,6),xx1,xx2,xx3 392 end if 393 394 elseif(isw.eq.8) then 395 396! Project stress values to nodes 397 398 call stcn2d(ix,sigl,shp2,jac,hr(nph),hr(nph+numnp),hr(ner), 399 & erav,lint,nel,9,numnp) 400 401 end if 402 403! Format statements 404 4052001 format(a1,20a4//5x,'Element Stresses'//' Elmt Matl', 406 1 ' 11-stress 22-stress 33-stress 12-stress', 407 2 ' 23-stress 13-stress'/16x,'1-coord 2-coord 3-coord ') 408 4092002 format(2i6,1p,6e11.3/12x,1p,6e11.3/12x,0p,3f11.5) 410 411 end 412