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 resultsmech_us3(co,kon,ipkon,lakon,ne,v, 20 & stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero, 21 & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr, 22 & iprestr,eme,iperturb,fn,iout,qa,vold,nmethod, 23 & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon, 24 & xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd, 25 & ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc, 26 & reltime,calcul_fn,calcul_qa,calcul_cauchy,nener, 27 & ikin,nal,ne0,thicke,emeini,i,ielprop,prop,t0g,t1g) 28! 29! calculates nal,qa,fn,xstiff,ener,eme,eei,stx for user element 1 30! 31! This is a beam type element. Reference: 32! Yunhua Luo, An Efficient 3D Timoshenko Beam Element with 33! Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1, 34! 2008, no. 3, 95-106 35! 36! 37! special case for which the beam axis goes through the 38! center of gravity of the cross section and the 1-direction 39! corresponds with a principal axis 40! 41! note that the strain components are Lagrange strain components, 42! the stress components are Piola-Kirchhoff components of the 43! second kind. For linear geometric calculations the strains 44! reduce to the infinitesimal strains 45! 46! 47! INPUT: 48! 49! co(1..3,i) coordinates of node i 50! kon(*) contains the topology of all elements. The 51! topology of element i starts at kon(ipkon(i)+1) 52! and continues until all nodes are covered. The 53! number of nodes depends on the element label 54! ipkon(i) points to the location in field kon preceding 55! the topology of element i 56! lakon(i) label of element i (character*8) 57! ne highest element number in the mesh 58! v(0..mi(2),i) value of the variables in node i at the end 59! of the present iteration 60! elcon elastic constants (cf. List of variables and 61! their meaning in the User's Manual) 62! nelcon integers describing the elastic constant fields 63! (cf. User's Manual) 64! rhcon density constants (cf. User's Manual) 65! nrhcon integers describing the density constant fields 66! (cf. User's Manual) 67! alcon thermal expansion constants (cf. User's Manual) 68! nalcon integers describing the thermal expansion 69! constants (cf. User's Manual) 70! alzero thermal expansion reference values (cf. User's Manual) 71! ielmat(i) material for element i 72! ielorien(i) orientation for element i 73! norien number of orientations 74! orab(7,*) description of all local coordinate systems. 75! (cf. List of variables and their meaning in the 76! User's manual) 77! ntmat_ maximum number of material temperature points 78! t0(i) temperature in node i at start of calculation 79! t1(i) temperature in node i at the end of the current 80! increment 81! ithermal(1..2) cf. List of variables and 82! their meaning in the User's Manual 83! prestr(i,j,k) residual stress component i in integration point j 84! of element k 85! iprestr if 0: no residual stresses 86! else: residual stresses 87! iperturb(*) describes the kind of nonlinearity of the 88! calculation, cf. User's Manual 89! iout if -2: v is assumed to be known and is used to 90! calculate strains, stresses..., no result output 91! corresponds to iout=-1 with in addition the 92! calculation of the internal energy density 93! if -1: v is assumed to be known and is used to 94! calculate strains, stresses..., no result 95! output; is used to take changes in SPC's and 96! MPC's at the start of a new increment or 97! iteration into account 98! if 0: v is calculated from the system solution 99! and strains, stresses.. are calculated, 100! no result output 101! if 1: v is calculated from the system solution and 102! strains, stresses.. are calculated, requested 103! results output 104! if 2: v is assumed to be known and is used to 105! calculate strains, stresses..., requested 106! results output 107! vold(0..mi(2),i) value of the variables in node i at the end 108! of the previous iteration 109! nmethod procedure: 110! 1: static analysis 111! 2: frequency analysis 112! 3: buckling analysis 113! 4: (linear or nonlinear) dynamic analysis 114! 5: steady state dynamics analysis 115! 6: Coriolis frequency calculation 116! 7: flutter frequency calculation 117! 8: magnetostatics 118! 9: magnetodynamics 119! 10: electromagnetic eigenvalue problems 120! 11: superelement creation or Green function 121! calculation 122! 12: sensitivity analysis 123! veold(j,i) time rate of variable j in node i at the end 124! of the previous iteration 125! dtime length of present time increment 126! time step time at the end of the present increment 127! ttime total time at the start of the present increment 128! plicon,nplicon fields describing isotropic hardening of 129! a plastic material or spring constants of 130! a nonlinear spring (cf. User's Manual) 131! plkcon,nplkcon fields describing kinematic hardening of 132! a plastic material or gap conductance 133! constants (cf. User's Manual) 134! xstateini(i,j,k) state variable component i in integration point j 135! of element k at the start of the present increment 136! xstate(i,j,k) state variable component i in integration point j 137! of element k at the end of the present increment 138! npmat_ maximum number of plastic constants 139! matname(i) name of material i 140! mi(1) max # of integration points per element (max 141! over all elements) 142! mi(2) max degree of freedom per node (max over all 143! nodes) in fields like v(0:mi(2))... 144! mi(3) max number of layers in the structure 145! ielas 0: no elastic iteration: irreversible effects 146! are allowed 147! 1: elastic iteration, i.e. no irreversible 148! deformation allowed 149! icmd not equal to 3: calculate stress and stiffness 150! 3: calculate only stress 151! ncmat_ max number of elastic constants 152! nstate_ max number of state variables in any integration 153! point 154! stiini(i,j,k) stress component i in integration point j 155! of element k at the start of the present 156! increment (= end of last increment) 157! vini(0..mi(2),i) value of the variables in node i at the start 158! of the present increment 159! enerini(j,k) internal energy density at integration point j 160! of element k at the start of the present increment 161! istep current step number 162! iinc current increment number within the actual step 163! reltime relative step time (between 0 and 1) 164! calcul_fn if 0: no nodal forces have to be calculated 165! else: nodal forces are required on output 166! calcul_qa if 0: no mean forces have to be calculated 167! else: mean forces are required on output 168! calcul_cauchy if 0: no Cauchy stresses are required 169! else: Cauchy stresses are required on output: have 170! to be calculated from the PK2 stresses 171! nener if 0: internal energy calculation is not required 172! else: internal energy is required on output 173! ikin if 0: kinetic energy calculation is not requred 174! else: kinetic energy is required on output 175! ne0 largest element number without contact elements (are 176! stored after all other elements) 177! thicke(j,i) layer thickness for layer j in element i 178! emeini(i,j,k) mechanical strain component i in integration point j 179! of element k at the start of the present increment 180! i actual element at stake 181! ielprop(i) points to the location in field prop preceding 182! the properties of element i 183! prop(*) contains the properties and some beam 184! elements (cf. User's Manual) 185! 186! 187! OUTPUT: 188! 189! stx(i,j,k) stress component i in integration point j 190! of element k at the end of the present 191! iteration 192! eme(i,j,k) mechanical strain component i in integration point j 193! of element k at the end of the present iteration 194! fn(j,i) nodal force component j in node i 195! qa(1..4) qa(1): average force 196! qa(2): average flux 197! ... cf. User's Manual 198! xstiff(i,j,k) stiffness (i=1...21) and conductivity 199! (i=22..27) constants in integration point j 200! of element k 201! ener(j,k) internal energy density at integration point j 202! of element k at the end of the present increment 203! eei(i,j,k) total strain component i in integration point j 204! of element k at the end of the present iteration 205! nal number of nodal force contributions 206! 207 implicit none 208 ! 209 character*8 lakon(*) 210 character*80 amat,matname(*) 211 ! 212 integer kon(*),konl(26),mi(*), 213 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*), 214 & ielorien(mi(3),*),ntmat_,ipkon(*),ne0,nelem, 215 & istep,iinc,ne,mattyp,ithermal(*),iprestr,i,j,k,m1,m2,jj, 216 & i1,kk,nener,indexe,nope,norien,iperturb(*),iout, 217 & nal,icmd,ihyper,nmethod,kode,imat,iorien,ielas, 218 & istiff,ncmat_,nstate_,ikin,ielprop(*), 219 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn, 220 & calcul_cauchy,calcul_qa,index,node,jjj,j1,j2 221 ! 222 real*8 co(3,*),v(0:mi(2),*),stiini(6,mi(1),*), 223 & stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*), 224 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*), 225 & alcon(0:6,ntmat_,*),vini(0:mi(2),*),q1, 226 & alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*), 227 & q(0:mi(2),26),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*), 228 & vold(0:mi(2),*),eloc(9),elconloc(21),eth(6),coords(3), 229 & ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*), 230 & veold(0:mi(2),*),e,un,um,tt,dl,qa(*),t0l,t1l,dtime,time,ttime, 231 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 232 & xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802), 233 & xstateini(nstate_,mi(1),*),tm(3,3),a,reltime,kap, 234 & thicke(mi(3),*),emeini(6,mi(1),*),aly,alz,bey,bez,xi(2), 235 & vlp(6,2),aa,Bs(2,18),Bb(3,18),Bm(3,18),tau(2),Smf(3), 236 & h,xg(3,3),x(3,3),ushell(18),fintg(18),ueg(18),pgauss(3), 237 & Ds(2,2),Qs(2,2),Qin(3,3),Dm(3,3),Db(3,3),temp_grad0, 238 & Kp(18,18),Km(18,18),tmg(18,18),Kshell(18,18), 239 & alp(3),temp_grad,thick_gp(2),ftherm(18),beta(6), 240 & d3,Ae,tmt(3,3),Ys(2),Emf(3),gpthick(3,2), 241 & Ethm(3),Ethf(3),Em(3),Ef(3),alph6(6),t0g(2,*),t1g(2,*), 242 & Dsi(2,2),Dmi(3,3),Dbi(3,3),di,dett,dettt 243 ! 244 nope = 3 245 ! 246 Bs(:,:) = 0.d0 247 Bb(:,:) = 0.d0 248 Bm(:,:) = 0.d0 249 beta(:) = 0.d0 250 ! 251 gpthick(1,1) = +1.d0 252 gpthick(2,1) = 0.d0 253 gpthick(3,1) = -1.d0 254 gpthick(1,2) = 2.d0/6.d0 255 gpthick(2,2) = 8.d0/6.d0 256 gpthick(3,2) = 2.d0/6.d0 257 ! 258 259 ! 260 di = 1.d0/3.d0 261 ! 262 indexe=ipkon(i) 263 ! 264 ! properties of the shell section 265 ! 266 index = ielprop(i) 267 h = prop(index+1) 268 dett = h/2.d0 269 dettt = h**3/8.d0 270 ! 271 ! material and orientation 272 ! 273 imat = ielmat(1,i) 274 amat = matname(imat) 275 ! 276 if(norien.gt.0) then 277 iorien = max(0,ielorien(1,i)) 278 else 279 iorien = 0 280 endif 281 ! 282 if(nelcon(1,imat).lt.0) then 283 ihyper = 1 284 else 285 ihyper = 0 286 endif 287 ! 288 kode = nelcon(1,imat) 289 ! 290 if(kode.eq.2) then 291 mattyp=1 292 else 293 mattyp=0 294 endif 295 ! 296 ! computation of the coordinates and displacements of the local nodes 297 ! 298 ! vl: u,v,w,rx,ry,rz 299 do k = 1,nope 300 konl(k) = kon(indexe+k) 301 do j = 1,3 302 xg(k,j) = co(j,konl(k)) ! global coordinates element nodes 303 vl(j,k) = v(j,konl(k)) ! uvw 304 vl(j+3,k) = v(j+3,konl(k)) ! rxyz 305 enddo 306 enddo 307! 308! gp temp gradient based nodal gradients 309 if(ithermal(1).ne.0) then 310 temp_grad = (t1g(1,konl(1))+t1g(1,konl(2))+t1g(1,konl(3)))*di 311 temp_grad0 = (t0g(1,konl(1))+t0g(1,konl(2))+t0g(1,konl(3)))*di 312 endif 313 ! 314 ! e0-frame base tm (3x3) -> tmg (18x18) 315 call us3_csys(xg,tm,tmg) ! | call us3_csys_cr(xg,tm,tmg) 316 call us3_xu(x,ushell,ueg,xg,tm,vl(1:6,1:3)) ! coords,displ in e0-frame & global displ 317 ! 318 ! stress and strains 319 ! 320 call us3_Bs(x,Bs) ! constant 321 call us3_Bb(x(:,1),x(:,2),Bb) ! constant 322 ! 323 Dm(:,:) = 0.d0 324 Db(:,:) = 0.d0 325 Ds(:,:) = 0.d0 326 ! 327 do jjj = 1,3 ! ---- start loop gauß points ---- 328 ! 329 ! natural coordinates z-dirc. -1,0,1 330 ! 331 if(jjj.eq.1) then 332 aa = -0.5d0*h 333 elseif(jjj.eq.2) then 334 aa = 0.d0 335 else 336 aa = 0.5d0*h 337 endif 338 ! 339 if(ithermal(1).ne.0) then 340 t0l = (t0(konl(1))+t0(konl(2))+t0(konl(3)))/3.d0 341 & +temp_grad0*aa 342 t1l = (t1(konl(1))+t1(konl(2))+t1(konl(3)))/3.d0 343 & +temp_grad*aa 344 endif 345 ! 346 istiff = 0 347 call us3_materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon, 348 & nalcon,imat,amat,iorien,pgauss,orab,ntmat_, 349 & elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper, 350 & istiff,elconloc,eth,kode,plicon,nplicon, 351 & plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,jjj, 352 & xstiff,ncmat_) 353 ! 354 if(mattyp.eq.1) then 355 e = elas(1) 356 un = elas(2) 357 rho = rhcon(1,1,imat) 358 alp(1) = eth(1) !alcon(1,1,imat) 359 alp(2) = eth(1) !alcon(1,1,imat) 360 alp(3) = 0.d0 361 elseif(mattyp.eq.2) then 362 write(*,*) '*ERROR in e_c3d_US3: no orthotropic material' 363 write(*,*) ' calculation for this type of element' 364 call exit(201) 365 else 366 write(*,*) '*ERROR in e_c3d_US3: no anisotropic material' 367 write(*,*) ' calculation for this type of element' 368 call exit(201) 369 endif 370 ! 371 call us3_linel_Qi(e,un,Qin,Qs) 372 Dm = Dm + Qin*gpthick(jjj,2)*dett 373 Db = Db + Qin*((gpthick(jjj,1))**2)*gpthick(jjj,2)*dettt 374 Ds = Ds + Qs*gpthick(jjj,2)*dett 375 ! 376 ! save elastic constants in xstff -> e_c3d_us3 377 ! 378 ! material data; for linear elastic materials 379 ! this includes the calculation of the stiffness 380 ! matrix 381 kode = 2 382 call linel(kode,mattyp,beta,eme,stre,elas,elconloc, 383 & iorien,orab,pgauss) 384 do m1=1,21 385 xstiff(m1,jjj,i) = elas(m1) ! elas for each gp saved in xstiff 386 enddo 387 ! 388 call bm_ANDES(x,Bm,Qin,h,di,di) 389 ! 390 ! Total strains 391 ! 392 Ys = matmul(Bs,ushell) 393 Em = matmul(Bm,ushell) 394 Ef = aa*(matmul(Bb,ushell)) 395 Emf = Em + Ef 396 ! 397 ! -> eei total strains out 398 ! 399 if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100)) then 400 eei(1,jjj,i) = Emf(1) 401 eei(2,jjj,i) = Emf(2) 402 eei(3,jjj,i) = 0.d0 403 eei(4,jjj,i) = Emf(3) 404 eei(5,jjj,i) = Ys(2) 405 eei(6,jjj,i) = Ys(1) 406 endif 407 ! 408 if(ithermal(1).ne.0) then 409 t0l = (t0(konl(1))+t0(konl(2))+t0(konl(3)))/3.d0+temp_grad0*aa 410 t1l = (t1(konl(1))+t1(konl(2))+t1(konl(3)))/3.d0+temp_grad*aa 411 ! 412 Ethm = alp*(t1l-t0l) ! mem 413 Ethf = alp*(temp_grad-temp_grad0)*aa ! bend. 414 Emf = Emf-(Ethf+Ethm) 415 endif 416 ! 417 ! mechanical strains 418 ! 419 eme(1,jjj,i) = Emf(1) 420 eme(2,jjj,i) = Emf(2) 421 eme(3,jjj,i) = 0.d0 422 eme(4,jjj,i) = Emf(3) 423 eme(5,jjj,i) = Ys(2) 424 eme(6,jjj,i) = Ys(1) 425 ! 426 Smf = matmul(Qin,Emf) 427 tau = matmul(Qs,Ys) 428 ! 429 stx(1,jjj,i) = Smf(1) ! sxx 430 stx(2,jjj,i) = Smf(2) ! syy 431 stx(3,jjj,i) = 0.d0 ! szz 432 stx(4,jjj,i) = Smf(3) ! txy 433 stx(5,jjj,i) = tau(1) ! tyz 434 stx(6,jjj,i) = tau(2) ! txz 435 ! 436 stre(1) = Smf(1) ! sxx 437 stre(2) = Smf(2) ! syy 438 stre(3) = 0.d0 ! szz 439 stre(4) = Smf(3) ! txy 440 stre(5) = tau(1) ! tyz 441 stre(6) = tau(2) ! txz 442 ! 443 enddo ! ---- end loop gauß points ---- 444 ! 445 ! Thermal loads 446 ! 447 if(ithermal(1).ne.0) then 448 ! 449 call us3_ae(x,Ae) ! element area 450 ftherm(:) = 0.d0 451 ! 452 do k = 1,3 453 ! 454 t0l = ((t0(1)+t0(2)+t0(3))/3.d0)!+temp_grad0*0.5*h*gpthick(k,1) 455 t1l = ((t1(1)+t1(2)+t1(3))/3.d0)!+temp_grad*0.5*h*gpthick(k,1) 456 ! 457 istiff=1 458 call us3_materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon, 459 & nalcon,imat,amat,iorien,coords,orab,ntmat_,elas,rho, 460 & i,ithermal,alzero,mattyp,t0l,t1l, 461 & ihyper,istiff,elconloc,eth,kode,plicon, 462 & nplicon,plkcon,nplkcon,npmat_, 463 & plconloc,mi(1),dtime,k, 464 & xstiff,alcon) 465 e = elas(1) 466 un = elas(2) 467 rho = rhcon(1,1,imat) 468 alp(1) = eth(1)!alcon(1,1,imat) 469 alp(2) = eth(2)!alcon(1,1,imat) 470 alp(3) = 0.d0 471 ! 472 call us3_linel_Qi(e,un,Qin,Qs) 473 ! 474 ftherm=ftherm-matmul(matmul(transpose(Bm),Qin),alp) 475 & *(t1l-t0l)*Ae*gpthick(k,2)*dett 476 ftherm=ftherm-matmul(matmul(transpose(Bb),Qin) ,alp)* 477 & (temp_grad-temp_grad0)*Ae*(gpthick(k,1)**2)*dettt*gpthick(k,2) 478 enddo 479! t0l = ((t0(1)-t0(2)-t0(3))/3.d0)!+temp_grad0*0.5*h*g_p(k,1) 480! t1l = ((t1(1)+t1(2)+t1(3))/3.d0)!+temp_grad*0.5*h*g_p(k,1) 481! ftherm=ftherm-matmul(matmul(transpose(Bm),Dm),alp) 482! & *(t1l-t0l)*Ae 483! ftherm=ftherm-matmul(matmul(transpose(Bb),Db) 484! & ,alp)*(temp_grad-temp_grad0)*Ae*h 485 do k = 1,nope 486 konl(k) = kon(indexe+k) 487 do j = 1,3 488 fn(j,konl(k)) = fn(j,konl(k)) +ftherm(j+(k-1)*6) ! f_th 489 fn(j+3,konl(k)) = fn(j+3,konl(k))+ftherm(j+3+(k-1)*6) ! r_th 490 enddo 491 enddo 492 endif 493 ! 494 ! Internal forces based on displacments 495 ! 496 if(calcul_fn.eq.1)then 497 ! 498 ! stiffness matirx (6 dofs 3 nodes -> 18x18) 499 ! 500 call us3_Kp(x,Db,Ds,Kp) ! plate part (CS-DSG) - in e0-frame 501 call us3_Km(x,Km,Dm/h,h) ! membrane part (ANDES) - in e0-frame 502 ! 503 Kshell = matmul(matmul(transpose(tmg),(Km+Kp)),tmg) 504 fintg = matmul(Kshell,ueg) 505 do k = 1,nope 506 konl(k) = kon(indexe+k) 507 do j = 1,3 508 fn(j,konl(k)) = fn(j,konl(k)) +fintg(j+(k-1)*6) ! fi 509 fn(j+3,konl(k)) = fn(j+3,konl(k))+fintg(j+3+(k-1)*6) ! ri 510 enddo 511 enddo 512 endif 513 ! 514 ! 515 ! q contains the contributions to the nodal force in the nodes 516 ! belonging to the element at stake from other elements (elements 517 ! already treated). These contributions have to be 518 ! subtracted to get the contributions attributable to the element 519 ! at stake only 520 ! 521 if(calcul_qa.eq.1) then 522 do m1=1,nope 523 do m2=1,3 524 qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1)) 525 enddo 526 enddo 527 nal=nal+3*nope 528 endif 529 ! 530 return 531 end 532