1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine e_c3d_u1(co,kon,lakonl,p1,p2,omx,bodyfx,nbody,s,sm, 20 & ff,nelem,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 21 & ielmat,ielorien,norien,orab,ntmat_, 22 & t0,t1,ithermal,vold,iperturb,nelemload, 23 & sideload,xload,nload,idist,sti,stx,iexpl,plicon, 24 & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime, 25 & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme, 26 & ttime,time,istep,iinc,coriolis,xloadold,reltime, 27 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 28 & ne0,ipkon,thicke, 29 & integerglob,doubleglob,tieset,istartset,iendset,ialset,ntie, 30 & nasym,ielprop,prop) 31! 32! computation of the element matrix and rhs for the user element 33! of type u1 34! 35! This is a beam type element. Reference: 36! Yunhua Luo, An Efficient 3D Timoshenko Beam Element with 37! Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1, 38! 2008, no. 3, 95-106 39! 40! special case for which the beam axis goes through the 41! center of gravity of the cross section and the 1-direction 42! corresponds with a principal axis 43! 44! 45! INPUT: 46! 47! co(1..3,i) coordinates of node i 48! kon(*) contains the topology of all elements. The 49! topology of element i starts at kon(ipkon(i)+1) 50! and continues until all nodes are covered. The 51! number of nodes depends on the element label 52! lakonl label of current element nelem (character*8) 53! p1(1..3) coordinates of a point on the rotation axis 54! (if applicable) 55! p2(1..3) unit vector on rotation axis (if applicable) 56! omx rotational speed square 57! bodyfx(1..3) acceleration vector 58! nbody number of body loads 59! nelem element number 60! nmethod procedure: 61! 1: static analysis 62! 2: frequency analysis 63! 3: buckling analysis 64! 4: (linear or nonlinear) dynamic analysis 65! 5: steady state dynamics analysis 66! 6: Coriolis frequency calculation 67! 7: flutter frequency calculation 68! 8: magnetostatics 69! 9: magnetodynamics 70! 10: electromagnetic eigenvalue problems 71! 11: superelement creation or Green function 72! calculation 73! 12: sensitivity analysis 74! elcon elastic constants (cf. List of variables and 75! their meaning in the User's Manual) 76! nelcon integers describing the elastic constant fields 77! (cf. User's Manual) 78! rhcon density constants (cf. User's Manual) 79! nrhcon integers describing the density constant fields 80! (cf. User's Manual) 81! alcon thermal expansion constants (cf. User's Manual) 82! nalcon integers describing the thermal expansion 83! constants (cf. User's Manual) 84! alzero thermal expansion reference values (cf. User's Manual) 85! ielmat(i) material for element i 86! ielorien(i) orientation for element i 87! norien number of orientations 88! orab(7,*) description of all local coordinate systems. 89! (cf. List of variables and their meaning in the 90! User's manual) 91! ntmat_ maximum number of material temperature points 92! t0(i) temperature in node i at start of calculation 93! t1(i) temperature in node i at the end of the current 94! increment 95! ithermal cf. ithermal(1) in the List of variables and 96! their meaning in the User's Manual 97! vold(0..mi(2),i) value of the variables in node i at the end 98! of the previous iteration 99! iperturb(*) describes the kind of nonlinearity of the 100! calculation, cf. User's Manual 101! nelemload field describing facial distributed load (cf. 102! User's Manual) 103! sideload field describing facial distributed load (cf. 104! User's Manual) 105! xload facial distributed load values (cf. 106! User's Manual) 107! nload number of facial distributed loads 108! idist 0: no distributed forces in the model 109! 1: distributed forces are present in the model 110! (body forces or thermal loads or residual 111! stresses or distributed facial loads) 112! sti(i,j,k) stress component i in integration point j 113! of element k at the end of the previous step 114! stx(i,j,k) stress component i in integration point j 115! of element k at the end of a static step 116! describing a reference buckling load 117! iexpl 0: structure: implicit dynamics, 118! fluid: incompressible 119! iexpl 1: structure: implicit dynamics, 120! fluid: compressible 121! iexpl 2: structure: explicit dynamics, 122! fluid: incompressible 123! iexpl 3: structure: explicit dynamics, 124! fluid: compressible 125! plicon,nplicon fields describing isotropic hardening of 126! a plastic material or spring constants of 127! a nonlinear spring (cf. User's Manual) 128! plkcon,nplkcon fields describing kinematic hardening of 129! a plastic material or gap conductance 130! constants (cf. User's Manual) 131! xstiff(i,j,k) stiffness (i=1...21) and conductivity 132! (i=22..27) constants in integration point j 133! of element k 134! npmat_ maximum number of plastic constants 135! dtime length of present time increment 136! matname(i) name of material i 137! mi(1) max # of integration points per element (max 138! over all elements) 139! mi(2) max degree of freedom per node (max over all 140! nodes) in fields like v(0:mi(2))... 141! mi(3) max number of layers in the structure 142! ncmat_ max number of elastic constants 143! mass(1) if 1: mass matrix needed, else not 144! mass(2) if 1: heat capacity matrix needed, else not 145! stiffness if 1: stiffness matrix needed, else not 146! buckling if 1: linear buckling calculation, else not 147! rhsi if 1: right hand side needed, else not 148! intscheme 1: integration point scheme for the calculation 149! of the initial acceleration in a dynamic 150! step 151! 0: any other case 152! ttime total time at the start of the present increment 153! time step time at the end of the present increment 154! istep current step number 155! iinc current increment number within the actual step 156! coriolis if 1: Coriolis forces are taken into account, 157! else not 158! xloadold load at the end of the previous increment 159! (cf. User's Manual) 160! reltime relative step time (between 0 and 1) 161! ipompc(i) pointer into field nodempc and coefmpc for 162! MPC i 163! nodempc field containing the nodes and dof's of all 164! MPC's (cf. User's Manual) 165! coefmpc field containing the coefficients of all 166! MPC's (cf. User's Manual) 167! nmpc number of MPC's in the model 168! ikmpc(i) field containing an integer code number 169! describing the dependent dof of some MPC: 170! ikmpc(i)=8*(node-1)+dir, where the dependent 171! dof is applied in direction dir of node "node" 172! ikmpc is sorted in ascending order 173! ilmpc(i) MPC number corresponding to ikmpc(i) (cf. 174! User's manual) 175! veold(j,i) time rate of variable j in node i at the end 176! of the previous iteration 177! ne0 element number before the first contact element; 178! contact elements are appended after ne0 179! ipkon(i) points to the location in field kon preceding 180! the topology of element i 181! thicke(j,i) thickness of layer j in node i 182! integerglob integer field needed for determining the 183! interface loading of a submodel 184! doubleglob real field needed for determining the 185! interface loading of a submodel 186! tieset(1..3,*) tie character information for tie i (cf. 187! User's Manual) 188! istartset integer describing set information (cf. 189! User's Manual) 190! iendset integer describing set information (cf. 191! User's Manual) 192! ialset integer describing set information (cf. 193! User's Manual) 194! ntie number of ties in the model 195! nasym 1: asymmetric contributions, else purely 196! symmetric 197! ielprop(i) points to the location in field prop preceding 198! the properties of element i 199! prop(*) contains the properties and some beam 200! elements (cf. User's Manual) 201! 202! 203! OUTPUT: 204! 205! s element stiffness matrix 206! sm element mass matrix 207! ff element load vector 208! xload facial distributed load values (cf. 209! User's Manual) 210! nmethod may be changed by the user, e.g. a 211! change to 0 for negative Jacobians leads 212! to a program exit in the calling program 213! for the other value: cf. above 214! 215 implicit none 216! 217 integer mass,stiffness,buckling,rhsi,coriolis 218! 219 character*8 lakonl 220 character*20 sideload(*) 221 character*80 matname(*),amat 222 character*81 tieset(3,*) 223! 224 integer konl(26),nelemload(2,*),nbody,nelem,mi(*),kon(*), 225 & ielprop(*),index,mattyp,ithermal(*),iperturb(*),nload,idist, 226 & i,j,k,i1,nmethod,kk,nelcon(2,*),nrhcon(*),nalcon(2,*), 227 & ielmat(mi(3),*),ielorien(mi(3),*),ipkon(*),indexe, 228 & ntmat_,nope,norien,ihyper,iexpl,kode,imat,iorien,istiff, 229 & ncmat_,intscheme,istep,iinc,ipompc(*),nodempc(3,*), 230 & nmpc,ikmpc(*),ilmpc(*),ne0,ndof,istartset(*),iendset(*), 231 & ialset(*),ntie,integerglob(*),nasym,nplicon(0:ntmat_,*), 232 & nplkcon(0:ntmat_,*),npmat_ 233! 234 real*8 co(3,*),xl(3,26),veold(0:mi(2),*),rho,s(60,60),bodyfx(3), 235 & ff(60),elconloc(21),coords(3),p1(3),elcon(0:ncmat_,ntmat_,*), 236 & p2(3),eth(6),rhcon(0:1,ntmat_,*),reltime,prop(*),tm(3,3), 237 & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),t0(*),t1(*), 238 & xloadold(2,*),vold(0:mi(2),*),xload(2,*),omx,e,un,um,tt, 239 & sm(60,60),sti(6,mi(1),*),stx(6,mi(1),*),t0l,t1l,coefmpc(*), 240 & elas(21),thicke(mi(3),*),doubleglob(*),dl,e2(3),e3(3), 241 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 242 & xstiff(27,mi(1),*),plconloc(802),dtime,ttime,time,tmg(12,12), 243 & a,xi11,xi12,xi22,xk,e1(3),offset1,offset2,y1,y2,y3,z1,z2,z3, 244 & sg(12,12) 245! 246! 247! 248 indexe=ipkon(nelem) 249! 250! material and orientation 251! 252 imat=ielmat(1,nelem) 253 amat=matname(imat) 254 if(norien.gt.0) then 255 iorien=max(0,ielorien(1,nelem)) 256 else 257 iorien=0 258 endif 259! 260 if(nelcon(1,imat).lt.0) then 261 ihyper=1 262 else 263 ihyper=0 264 endif 265! 266! properties of the cross section 267! 268 index=ielprop(nelem) 269 a=prop(index+1) 270 xi11=prop(index+2) 271 xi12=prop(index+3) 272 xi22=prop(index+4) 273 xk=prop(index+5) 274 e2(1)=prop(index+6) 275 e2(2)=prop(index+7) 276 e2(3)=prop(index+8) 277 offset1=prop(index+9) 278 offset2=prop(index+10) 279! 280 nope=2 281 ndof=6 282! 283 if(dabs(xi12).gt.0.d0) then 284 write(*,*) '*ERROR in e_c3d_u1: no nonzero cross moment' 285 write(*,*) ' of inertia for this type of element' 286 call exit(201) 287 endif 288! 289 if(dabs(offset1).gt.0.d0) then 290 write(*,*) '*ERROR in e_c3d_u1: no offset in 1-direction' 291 write(*,*) ' for this type of element' 292 call exit(201) 293 endif 294! 295 if(dabs(offset2).gt.0.d0) then 296 write(*,*) '*ERROR in e_c3d_u1: no offset in 2-direction' 297 write(*,*) ' for this type of element' 298 call exit(201) 299 endif 300! 301! computation of the coordinates of the local nodes 302! 303 do i=1,nope 304 konl(i)=kon(indexe+i) 305 do j=1,3 306 xl(j,i)=co(j,konl(i)) 307 enddo 308 enddo 309! 310! local axes e1-e2-e3 (e2 is given by the user (in the input deck 311! this is called e1), e1 is parallel to the beam axis) 312! 313 do j=1,3 314 e1(j)=xl(j,2)-xl(j,1) 315 enddo 316! 317 dl=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3)) 318! 319 do j=1,3 320 e1(j)=e1(j)/dl 321 enddo 322! 323! e3 = e1 x e2 324! 325 e3(1)=e1(2)*e2(3)-e1(3)*e2(2) 326 e3(2)=e1(3)*e2(1)-e1(1)*e2(3) 327 e3(3)=e1(1)*e2(2)-e1(2)*e2(1) 328! 329! transformation matrix from the global into the local system 330! 331 do j=1,3 332 tm(1,j)=e1(j) 333 tm(2,j)=e2(j) 334 tm(3,j)=e3(j) 335 enddo 336! 337! initialisation for distributed forces 338! 339 if(rhsi.eq.1) then 340 if(idist.ne.0) then 341 do i=1,ndof*nope 342 ff(i)=0.d0 343 enddo 344 endif 345 endif 346! 347! displacements for 2nd order static and modal theory 348! 349 if(((iperturb(1).eq.1).or.(iperturb(2).eq.1)).and. 350 & (stiffness.eq.1).and.(buckling.eq.0)) then 351 write(*,*) '*ERROR in e_c3d_u1: no second order' 352 write(*,*) ' calculation for this type of element' 353 call exit(201) 354 endif 355! 356! initialisation of sm 357! 358 if((mass.eq.1).or.(buckling.eq.1).or.(coriolis.eq.1)) then 359 write(*,*) '*ERROR in e_c3d_u1: no dynamic or buckling' 360 write(*,*) ' calculation for this type of element' 361 call exit(201) 362 endif 363! 364! initialisation of s 365! 366 do i=1,ndof*nope 367 do j=1,ndof*nope 368 s(i,j)=0.d0 369 enddo 370 enddo 371! 372! computation of the matrix 373! 374! calculating the temperature in the integration 375! point 376! 377 t0l=0.d0 378 t1l=0.d0 379 if(ithermal(1).eq.1) then 380 do i1=1,nope 381 t0l=t0l+t0(konl(i1))/2.d0 382 t1l=t1l+t1(konl(i1))/2.d0 383 enddo 384 elseif(ithermal(1).ge.2) then 385 write(*,*) '*ERROR in e_c3d_u1: no thermal' 386 write(*,*) ' calculation for this type of element' 387 call exit(201) 388 endif 389 tt=t1l-t0l 390! 391! calculating the coordinates of the integration point 392! for material orientation purposes (for cylindrical 393! coordinate systems) 394! 395 if(iorien.gt.0) then 396 write(*,*) '*ERROR in e_c3d_u1: no orientation' 397 write(*,*) ' calculation for this type of element' 398 call exit(201) 399 endif 400! 401! for deformation plasticity: calculating the Jacobian 402! and the inverse of the deformation gradient 403! needed to convert the stiffness matrix in the spatial 404! frame of reference to the material frame 405! 406 kode=nelcon(1,imat) 407! 408 if(kode.eq.2) then 409 mattyp=1 410 else 411 mattyp=0 412 endif 413! 414! material data and local stiffness matrix 415! 416 istiff=0 417 call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 418 & imat,amat,iorien,coords,orab,ntmat_,elas,rho, 419 & nelem,ithermal,alzero,mattyp,t0l,t1l, 420 & ihyper,istiff,elconloc,eth,kode,plicon, 421 & nplicon,plkcon,nplkcon,npmat_, 422 & plconloc,mi(1),dtime,kk, 423 & xstiff,ncmat_) 424! 425 if(mattyp.eq.1) then 426 e=elconloc(1) 427 un=elconloc(2) 428 um=e/(1.d0+un) 429 um=um/2.d0 430 elseif(mattyp.eq.2) then 431 write(*,*) '*ERROR in e_c3d_u1: no orthotropic material' 432 write(*,*) ' calculation for this type of element' 433 call exit(201) 434 else 435 write(*,*) '*ERROR in e_c3d_u1: no anisotropic material' 436 write(*,*) ' calculation for this type of element' 437 call exit(201) 438 endif 439! 440! initialization for the body forces 441! 442 if(rhsi.eq.1) then 443 if(nbody.ne.0) then 444 write(*,*) '*ERROR in e_c3d_u1: no body forces' 445 write(*,*) ' for this type of element' 446 call exit(201) 447 endif 448 endif 449! 450 if(buckling.eq.1) then 451 write(*,*) '*ERROR in e_c3d_u1: no buckling ' 452 write(*,*) ' calculation for this type of element' 453 call exit(201) 454 endif 455! 456! determination of the stiffness, and/or mass and/or 457! buckling matrix 458! 459 if((stiffness.eq.1).or.(mass.eq.1).or.(buckling.eq.1).or. 460 & (coriolis.eq.1)) then 461! 462 if(((iperturb(1).ne.1).and.(iperturb(2).ne.1)).or. 463 & (buckling.eq.1)) then 464! 465! stiffness matrix 466! 467 y1=xk*um*a*e*xi11*(12.d0*e*xi11+xk*um*a*dl*dl) 468 y2=(12.d0*e*xi11-xk*um*a*dl*dl)**2 469 y3=4.d0*e*xi11*((xk*um*a)**2*dl**4+ 470 & 3.d0*xk*um*a*dl*dl*e*xi11+ 471 & 36.d0*(e*xi11)**2) 472 z1=xk*um*a*e*xi22*(12.d0*e*xi22+xk*um*a*dl*dl) 473 z2=(12.d0*e*xi22-xk*um*a*dl*dl)**2 474 z3=4.d0*e*xi22*((xk*um*a)**2*dl**4+ 475 & 3.d0*xk*um*a*dl*dl*e*xi22+ 476 & 36.d0*(e*xi22)**2) 477! 478! stiffness matrix S' in local coordinates 479! 480 s(1,1)=e*a/dl 481 s(1,7)=-s(1,1) 482 s(2,2)=12.d0*y1/(dl*y2) 483 s(2,6)=6.d0*y1/y2 484 s(2,8)=-s(2,2) 485 s(2,12)=s(2,6) 486 s(3,3)=12.d0*z1/(dl*z2) 487 s(3,5)=-6.d0*z1/z2 488 s(3,9)=s(3,3) 489 s(3,11)=s(3,5) 490 s(4,4)=um*(xi11+xi22)/dl 491 s(4,10)=-s(4,4) 492 s(5,5)=z3/(dl*z2) 493 s(5,9)=6.d0*z1/z2 494 s(5,11)=-2.d0*e*xi22*(72.d0*(e*xi22)**2- 495 & (xk*um*a)**2*dl**4- 496 & 30.d0*xk*um*a*dl*dl*e*xi22)/(dl*z2) 497 s(6,6)=y3/(dl*y2) 498 s(6,8)=-6.d0*y1/y2 499 s(6,12)=-2.d0*e*xi11*(-(xk*um*a)**2*dl**4- 500 & 30.d0*xk*um*a*dl*dl*e*xi11+ 501 & 72.d0*(e*xi11)**2)/(dl*y2) 502 s(1,7)=-s(1,1) 503 s(7,7)=s(1,1) 504 s(8,8)=12.d0*y1/(dl*y2) 505 s(8,12)=-6.d0*y1/y2 506 s(9,9)=12.d0*z1/(dl*z2) 507 s(9,11)=6.d0*z1/z2 508 s(10,10)=s(4,4) 509 s(11,11)=z3/(dl*z2) 510 s(12,12)=y3/(dl*y2) 511! 512! completing the symmetric part 513! 514 do i=1,12 515 do j=1,i 516 s(i,j)=s(j,i) 517 enddo 518 enddo 519! 520! 12 x 12 transformation matrix 521! 522 do i=1,12 523 do j=1,12 524 tmg(i,j)=0.d0 525 enddo 526 enddo 527 do i=1,3 528 do j=1,3 529 tmg(i,j)=tm(i,j) 530 tmg(i+3,j+3)=tm(i,j) 531 tmg(i+6,j+6)=tm(i,j) 532 tmg(i+9,j+9)=tm(i,j) 533 enddo 534 enddo 535! 536! stiffness matrix in global coordinates: S = T^T*S'*T 537! 538 do i=1,12 539 do j=1,12 540 sg(i,j)=0.d0 541 do k=1,12 542 sg(i,j)=sg(i,j)+s(i,k)*tmg(k,j) 543 enddo 544 enddo 545 enddo 546! 547! only upper triangular matrix 548! 549 do i=1,12 550 do j=i,12 551 s(i,j)=0.d0 552 do k=1,12 553 s(i,j)=s(i,j)+tmg(k,i)*sg(k,j) 554 enddo 555 enddo 556 enddo 557! 558 endif 559! 560 endif 561! 562 return 563 end 564 565