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 rubber(elconloc,elas,emec,kode,didc, 20 & d2idc2,dibdc,d2ibdc2,dudc,d2udc2,dldc,d2ldc2,dlbdc,d2lbdc2, 21 & ithermal,icmd,beta,stre) 22! 23! calculates stiffness and stresses for rubber and elastomeric 24! foam materials 25! 26! icmd=3: stress at mechanical strain 27! else: stress and stiffness matrix at mechanical strain 28! 29 implicit none 30! 31 integer ogden,hyperfoam,taylor 32! 33 integer nelconst,kode,kk(84),i,j,k,l,m,nt,icmd,istart,iend, 34 & nc,n,ithermal(*),ii,jj,mm,neig 35! 36 real*8 elconloc(*),elas(*),emec(*),didc(3,3,3), 37 & d2idc2(3,3,3,3,3),dibdc(3,3,3),d2ibdc2(3,3,3,3,3),dudc(3,3), 38 & d2udc2(3,3,3,3),dldc(3,3,3),d2ldc2(3,3,3,3,3),dlbdc(3,3,3), 39 & d2lbdc2(3,3,3,3,3),v1,v2,v3,c(3,3),cinv(3,3),d(3,3),djth, 40 & coef,bb,cc,cm,cn,tt,pi,dd,al(3),v1b,v2b,v3b, 41 & alb(3),beta(*),v33,v36,all(3),term,stre(*),total,coefa, 42 & coefb,coefd,coefm,constant(21) 43! 44 kk=(/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3, 45 & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,3,3,1,3, 46 & 1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,1,2,2,3,1,3,2,3, 47 & 2,3,2,3/) 48! 49! copy the elastic constants into a new field, such that 50! they can be mixed without influencing the field in the 51! calling program 52! 53 do i=1,21 54 constant(i)=elconloc(i) 55 enddo 56! 57! type of hyperelastic law; taylor stands for everything 58! which involves parts of a taylor expansion in terms of the 59! reduced Green deformation invariants 60! 61 ogden=0 62 hyperfoam=0 63 taylor=0 64 if((kode.lt.-3).and.(kode.gt.-7)) then 65 ogden=1 66 elseif((kode.lt.-14).and.(kode.gt.-18)) then 67 hyperfoam=1 68 else 69 taylor=1 70 endif 71! 72c if(icmd.eq.1) then 73 istart=1 74 iend=1 75! 76! reclassifying some classes of hyperelastic materials as 77! subclasses of the polynomial model 78! 79 if(((kode.lt.-1).and.(kode.gt.-4)).or. 80 & ((kode.lt.-6).and.(kode.gt.-13)).or. 81 & (kode.eq.-14)) then 82 if(kode.eq.-2) then 83 kode=-7 84 nelconst=1 85 elseif((kode.eq.-3).or.(kode.eq.-10)) then 86 constant(3)=constant(2) 87 constant(2)=0.d0 88 kode=-7 89 nelconst=1 90 elseif(kode.eq.-11) then 91 constant(7)=constant(4) 92 constant(6)=constant(3) 93 constant(5)=0.d0 94 constant(4)=0.d0 95 constant(3)=constant(2) 96 constant(2)=0.d0 97 kode=-8 98 nelconst=2 99 elseif((kode.eq.-12).or.(kode.eq.-14)) then 100 constant(12)=constant(6) 101 constant(11)=constant(5) 102 constant(10)=constant(4) 103 constant(9)=0.d0 104 constant(8)=0.d0 105 constant(7)=0.d0 106 constant(6)=constant(3) 107 constant(5)=0.d0 108 constant(4)=0.d0 109 constant(3)=constant(2) 110 constant(2)=0.d0 111 kode=-9 112 nelconst=3 113 elseif(kode.eq.-7) then 114 nelconst=1 115 elseif(kode.eq.-8) then 116 nelconst=2 117 elseif(kode.eq.-9) then 118 nelconst=3 119 endif 120 endif 121! 122! major loop 123! 124 do ii=istart,iend 125! 126! calculation of the Green deformation tensor for the total 127! strain and the thermal strain 128! 129 do i=1,3 130 c(i,i)=emec(i)*2.d0+1.d0 131 enddo 132 c(1,2)=2.d0*emec(4) 133 c(1,3)=2.d0*emec(5) 134 c(2,3)=2.d0*emec(6) 135! 136! calculation of the invariants of c 137! 138 v1=c(1,1)+c(2,2)+c(3,3) 139 v2=c(2,2)*c(3,3)+c(1,1)*c(3,3)+c(1,1)*c(2,2)- 140 & (c(2,3)*c(2,3)+c(1,3)*c(1,3)+c(1,2)*c(1,2)) 141 v3=c(1,1)*(c(2,2)*c(3,3)-c(2,3)*c(2,3)) 142 & -c(1,2)*(c(1,2)*c(3,3)-c(1,3)*c(2,3)) 143 & +c(1,3)*(c(1,2)*c(2,3)-c(1,3)*c(2,2)) 144 v33=v3**(-1.d0/3.d0) 145 v36=v3**(-1.d0/6.d0) 146! 147! calculation of the thermal strain jacobian 148! (not really needed) 149! 150 djth=1.d0 151! 152! inversion of c 153! 154 cinv(1,1)=(c(2,2)*c(3,3)-c(2,3)*c(2,3))/v3 155 cinv(2,2)=(c(1,1)*c(3,3)-c(1,3)*c(1,3))/v3 156 cinv(3,3)=(c(1,1)*c(2,2)-c(1,2)*c(1,2))/v3 157 cinv(1,2)=(c(1,3)*c(2,3)-c(1,2)*c(3,3))/v3 158 cinv(1,3)=(c(1,2)*c(2,3)-c(2,2)*c(1,3))/v3 159 cinv(2,3)=(c(1,2)*c(1,3)-c(1,1)*c(2,3))/v3 160 cinv(2,1)=cinv(1,2) 161 cinv(3,1)=cinv(1,3) 162 cinv(3,2)=cinv(2,3) 163! 164! creation of the delta Dirac matrix d 165! 166 do j=1,3 167 do i=1,3 168 d(i,j)=0.d0 169 enddo 170 enddo 171 do i=1,3 172 d(i,i)=1.d0 173 enddo 174! 175! derivative of the c-invariants with respect to c(k,l) 176! 177 do l=1,3 178 do k=1,l 179 didc(k,l,1)=d(k,l) 180 didc(k,l,2)=v1*d(k,l)-c(k,l) 181 didc(k,l,3)=v3*cinv(k,l) 182 enddo 183 enddo 184! 185! second derivative of the c-invariants w.r.t. c(k,l) 186! and c(m,n) 187! 188 if(icmd.ne.3) then 189 nt=0 190 do i=1,21 191 k=kk(nt+1) 192 l=kk(nt+2) 193 m=kk(nt+3) 194 n=kk(nt+4) 195 nt=nt+4 196 d2idc2(k,l,m,n,1)=0.d0 197 d2idc2(k,l,m,n,2)=d(k,l)*d(m,n)- 198 & (d(k,m)*d(l,n)+d(k,n)*d(l,m))/2.d0 199 d2idc2(k,l,m,n,3)=v3*(cinv(m,n)*cinv(k,l)- 200 & (cinv(k,m)*cinv(n,l)+cinv(k,n)*cinv(m,l))/2.d0) 201 enddo 202 endif 203! 204! derivatives for the reduced invariants used in rubber materials 205! 206 v1b=v1*v33 207 v2b=v2*v33*v33 208 v3b=dsqrt(v3)/djth 209! 210! first derivative of the reduced c-invariants w.r.t. c(k,l) 211! 212 do l=1,3 213 do k=1,l 214 if(taylor.eq.1) then 215 dibdc(k,l,1)=-v33**4*v1*didc(k,l,3)/3.d0 216 & +v33*didc(k,l,1) 217 dibdc(k,l,2)=-2.d0*v33**5*v2*didc(k,l,3)/3.d0 218 & +v33**2*didc(k,l,2) 219 endif 220 dibdc(k,l,3)=didc(k,l,3)/(2.d0*dsqrt(v3)*djth) 221 enddo 222 enddo 223! 224! second derivative of the reduced c-invariants w.r.t. c(k,l) 225! and c(m,n) 226! 227 if(icmd.ne.3) then 228 nt=0 229 do i=1,21 230 k=kk(nt+1) 231 l=kk(nt+2) 232 m=kk(nt+3) 233 n=kk(nt+4) 234 nt=nt+4 235 if(taylor.eq.1) then 236 d2ibdc2(k,l,m,n,1)=4.d0/9.d0*v33**7*v1*didc(k,l,3) 237 & *didc(m,n,3)-v33**4/3.d0*(didc(m,n,1)*didc(k,l,3) 238 & +didc(k,l,1)*didc(m,n,3))-v33**4/3.d0*v1* 239 & d2idc2(k,l,m,n,3)+v33*d2idc2(k,l,m,n,1) 240 d2ibdc2(k,l,m,n,2)=10.d0*v33**8/9.d0*v2*didc(k,l,3) 241 & *didc(m,n,3)-2.d0*v33**5/3.d0*(didc(m,n,2) 242 & *didc(k,l,3) 243 & +didc(k,l,2)*didc(m,n,3))-2.d0*v33**5/3.d0*v2* 244 & d2idc2(k,l,m,n,3)+v33**2*d2idc2(k,l,m,n,2) 245 endif 246 d2ibdc2(k,l,m,n,3)=-didc(k,l,3)*didc(m,n,3)/ 247 & (4.d0*djth*v3**1.5d0)+d2idc2(k,l,m,n,3)/ 248 & (2.d0*dsqrt(v3)*djth) 249 enddo 250 endif 251! 252! calculation of the principal stretches for the Ogden model and 253! hyperfoam materials 254! 255 if((ogden.eq.1).or.(hyperfoam.eq.1)) then 256! 257! taking the thermal jacobian into account 258! 259 if((kode.lt.-14).and.(kode.gt.-18)) then 260 dd=djth**(1.d0/3.d0) 261 else 262 dd=1.d0 263 endif 264! 265 pi=4.d0*datan(1.d0) 266! 267! determining the eigenvalues of c (Simo & Hughes) and taking 268! the square root to obtain the principal stretches 269! 270! neig is the number of different eigenvalues 271! 272 neig=3 273! 274 bb=v2-v1*v1/3.d0 275 cc=-2.d0*v1**3/27.d0+v1*v2/3.d0-v3 276 if(dabs(bb).le.1.d-10) then 277 if(dabs(cc).gt.1.d-10) then 278 al(1)=-cc**(1.d0/3.d0) 279 else 280 al(1)=0.d0 281 endif 282 al(2)=al(1) 283 al(3)=al(1) 284 neig=1 285 else 286 cm=2.d0*dsqrt(-bb/3.d0) 287 cn=3.d0*cc/(cm*bb) 288 if(dabs(cn).gt.1.d0) then 289 if(cn.gt.1.d0) then 290 cn=1.d0 291 else 292 cn=-1.d0 293 endif 294 endif 295 tt=datan2(dsqrt(1.d0-cn*cn),cn)/3.d0 296 al(1)=dcos(tt) 297 al(2)=dcos(tt+2.d0*pi/3.d0) 298 al(3)=dcos(tt+4.d0*pi/3.d0) 299! 300! check for two equal eigenvalues 301! 302 if((dabs(al(1)-al(2)).lt.1.d-5).or. 303 & (dabs(al(1)-al(3)).lt.1.d-5).or. 304 & (dabs(al(2)-al(3)).lt.1.d-5)) neig=2 305 al(1)=cm*al(1) 306 al(2)=cm*al(2) 307 al(3)=cm*al(3) 308 endif 309 do i=1,3 310 al(i)=dsqrt(al(i)+v1/3.d0) 311 all(i)=(6.d0*al(i)**5-4.d0*v1*al(i)**3+2.d0*al(i)*v2)*dd 312 enddo 313! 314! first derivative of the principal stretches w.r.t. c(k,l) 315! 316 if(neig.eq.3) then 317! 318! three different principal stretches 319! 320 do i=1,3 321 do l=1,3 322 do k=1,l 323 dldc(k,l,i)=(al(i)**4*didc(k,l,1) 324 & -al(i)**2*didc(k,l,2)+didc(k,l,3))/all(i) 325 enddo 326 enddo 327 enddo 328 elseif(neig.eq.1) then 329! 330! three equal principal stretches 331! 332 do i=1,3 333 do l=1,3 334 do k=1,l 335 dldc(k,l,i)=didc(k,l,1)/(6.d0*al(i)) 336 enddo 337 enddo 338 enddo 339 else 340! 341! two equal principal stretches 342! 343 do i=1,3 344 do l=1,3 345 do k=1,l 346 dldc(k,l,i)=(dcos(tt+(i-1)*2.d0*pi/3.d0)* 347 & (2.d0*v1*didc(k,l,1)-3.d0*didc(k,l,2))/ 348 & (3.d0*dsqrt(v1*v1-3.d0*v2))+didc(k,l,1)/3.d0)/ 349 & (2.d0*al(i)) 350 enddo 351 enddo 352 enddo 353 endif 354! 355! second derivative of the principal stretches w.r.t. c(k,l) 356! and c(m,n) 357! 358 if(icmd.ne.3) then 359 if(neig.eq.3) then 360! 361! three different principal stretches 362! 363 do i=1,3 364 nt=0 365 do j=1,21 366 k=kk(nt+1) 367 l=kk(nt+2) 368 m=kk(nt+3) 369 n=kk(nt+4) 370 nt=nt+4 371 d2ldc2(k,l,m,n,i)=(-30.d0*al(i)**4 372 & *dldc(k,l,i)*dldc(m,n,i)+al(i)**4 373 & *d2idc2(k,l,m,n,1) 374 & +4.d0*al(i)**3*(didc(k,l,1)*dldc(m,n,i) 375 & +didc(m,n,1) 376 & *dldc(k,l,i))+12.d0*v1*al(i)**2*dldc(k,l,i)* 377 & dldc(m,n,i)-d2idc2(k,l,m,n,2)*al(i)**2-2.d0 378 & *al(i)* 379 & didc(k,l,2)*dldc(m,n,i)-2.d0*v2*dldc(k,l,i)* 380 & dldc(m,n,i)-2.d0*al(i)*didc(m,n,2)*dldc(k,l,i) 381 & +d2idc2(k,l,m,n,3))/all(i) 382 enddo 383 enddo 384 elseif(neig.eq.1) then 385! 386! three equal principal stretches 387! 388 do i=1,3 389 nt=0 390 do j=1,21 391 k=kk(nt+1) 392 l=kk(nt+2) 393 m=kk(nt+3) 394 n=kk(nt+4) 395 nt=nt+4 396 d2ldc2(k,l,m,n,i)=(d2idc2(k,l,m,n,1)/6.d0 397 & -dldc(k,l,i)*dldc(m,n,i))/al(i) 398 enddo 399 enddo 400 else 401! 402! two equal principal stretches 403! 404 do i=1,3 405 nt=0 406 do j=1,21 407 k=kk(nt+1) 408 l=kk(nt+2) 409 m=kk(nt+3) 410 n=kk(nt+4) 411 nt=nt+4 412 d2ldc2(k,l,m,n,i)=(dcos(tt+(i-1)*2.d0*pi/3.d0)* 413 & (-(2.d0*v1*didc(k,l,1)-3.d0*didc(k,l,2))* 414 & (2.d0*v1*didc(m,n,1)-3.d0*didc(m,n,2))/ 415 & (6.d0*(v1*v1-3.d0*v2)**1.5d0)+ 416 & (2.d0*didc(k,l,1)*didc(m,n,1)+2.d0*v1* 417 & d2idc2(k,l,m,n,1)-3.d0*d2idc2(k,l,m,n,2))/ 418 & (3.d0*dsqrt(v1*v1-3.d0*v2))) 419 & +d2idc2(k,l,m,n,1)/3.d0)/(2.d0*al(i))- 420 & dldc(k,l,i)*dldc(m,n,i)/al(i) 421 enddo 422 enddo 423 endif 424 endif 425! 426! reduced principal stretches (Ogden model) 427! 428 if(ogden.eq.1) then 429! 430! calculation of the reduced principal stretches 431! 432 do i=1,3 433 alb(i)=al(i)*v36 434 enddo 435! 436! first derivative of the reduced principal stretches 437! w.r.t. c(k,l) 438! 439 do i=1,3 440 do l=1,3 441 do k=1,l 442 dlbdc(k,l,i)=-v36**7*al(i)*didc(k,l,3)/6.d0 443 & +v36*dldc(k,l,i) 444 enddo 445 enddo 446 enddo 447! 448! second derivative of the reduced principal stretches w.r.t. 449! c(k,l) and c(m,n) 450! 451 if(icmd.ne.3) then 452 do i=1,3 453 nt=0 454 do j=1,21 455 k=kk(nt+1) 456 l=kk(nt+2) 457 m=kk(nt+3) 458 n=kk(nt+4) 459 nt=nt+4 460 d2lbdc2(k,l,m,n,i)=7.d0*v36**13*al(i) 461 & *didc(k,l,3)*didc(m,n,3)/36.d0-v36**7/6.d0 462 & *(dldc(m,n,i)*didc(k,l,3)+al(i) 463 & *d2idc2(k,l,m,n,3)+dldc(k,l,i)*didc(m,n,3)) 464 & +v36*d2ldc2(k,l,m,n,i) 465 enddo 466 enddo 467 endif 468! 469 endif 470 endif 471! 472! calculation of the local stiffness matrix, and, if appropriate, 473! the stresses 474! 475! Polynomial model 476! 477 if((kode.lt.-6).and.(kode.gt.-10)) then 478! 479! first derivative of U w.r.t. c(k,l) 480! 481 do l=1,3 482 do k=1,l 483 dudc(k,l)=0.d0 484 enddo 485 enddo 486! 487 nc=0 488 do m=1,nelconst 489 do j=0,m 490 i=m-j 491 nc=nc+1 492 coef=constant(nc) 493 if(dabs(coef).lt.1.d-20) cycle 494 do l=1,3 495 do k=1,l 496 total=0.d0 497 if(i.gt.0) then 498 term=dibdc(k,l,1) 499 if(i.gt.1) term=i*term*(v1b-3.d0)**(i-1) 500 if(j.gt.0) term=term*(v2b-3.d0)**j 501 total=total+term 502 endif 503 if(j.gt.0) then 504 term=dibdc(k,l,2) 505 if(i.gt.0) term=term*(v1b-3.d0)**i 506 if(j.gt.1) term=j*term*(v2b-3.d0)**(j-1) 507 total=total+term 508 endif 509 dudc(k,l)=dudc(k,l)+total*coef 510 enddo 511 enddo 512 enddo 513 enddo 514 do m=1,nelconst 515 nc=nc+1 516 coef=constant(nc) 517 do l=1,3 518 do k=1,l 519 dudc(k,l)=dudc(k,l)+2.d0*m*(v3b-1.d0)** 520 & (2*m-1)*dibdc(k,l,3)/coef 521 enddo 522 enddo 523 enddo 524! 525! tangent stiffness matrix 526! second derivative of U w.r.t. c(k,l) and c(m,n) 527! 528 if(icmd.ne.3) then 529 nt=0 530 do i=1,21 531 k=kk(nt+1) 532 l=kk(nt+2) 533 m=kk(nt+3) 534 n=kk(nt+4) 535 nt=nt+4 536 d2udc2(k,l,m,n)=0.d0 537 enddo 538 nc=0 539 do mm=1,nelconst 540 do j=0,mm 541 i=mm-j 542 nc=nc+1 543 coef=constant(nc) 544 if(dabs(coef).lt.1.d-20) cycle 545 nt=0 546 do jj=1,21 547 k=kk(nt+1) 548 l=kk(nt+2) 549 m=kk(nt+3) 550 n=kk(nt+4) 551 nt=nt+4 552 total=0.d0 553 if(i.gt.1) then 554 term=dibdc(k,l,1)*dibdc(m,n,1)*i*(i-1) 555 if(i.gt.2) term=term*(v1b-3.d0)**(i-2) 556 if(j.gt.0) term=term*(v2b-3.d0)**j 557 total=total+term 558 endif 559 if((i.gt.0).and.(j.gt.0)) then 560 term=dibdc(k,l,1)*dibdc(m,n,2)+ 561 & dibdc(m,n,1)*dibdc(k,l,2) 562 if(i.gt.1) term=i*term*(v1b-3.d0)**(i-1) 563 if(j.gt.1) term=j*term*(v2b-3.d0)**(j-1) 564 total=total+term 565 endif 566 if(i.gt.0) then 567 term=d2ibdc2(k,l,m,n,1) 568 if(i.gt.1) term=i*term*(v1b-3.d0)**(i-1) 569 if(j.gt.0) term=term*(v2b-3.d0)**j 570 total=total+term 571 endif 572 if(j.gt.1) then 573 term=dibdc(k,l,2)*dibdc(m,n,2)*j*(j-1) 574 if(i.gt.0) term=term*(v1b-3.d0)**i 575 if(j.gt.2) term=term*(v2b-3.d0)**(j-2) 576 total=total+term 577 endif 578 if(j.gt.0) then 579 term=d2ibdc2(k,l,m,n,2) 580 if(i.gt.0) term=term*(v1b-3.d0)**i 581 if(j.gt.1) term=j*term*(v2b-3.d0)**(j-1) 582 total=total+term 583 endif 584 d2udc2(k,l,m,n)=d2udc2(k,l,m,n)+total*coef 585 enddo 586 enddo 587 enddo 588! 589 do mm=1,nelconst 590 nc=nc+1 591 coef=constant(nc) 592 nt=0 593 do i=1,21 594 k=kk(nt+1) 595 l=kk(nt+2) 596 m=kk(nt+3) 597 n=kk(nt+4) 598 nt=nt+4 599 if(mm.eq.1) then 600 term=(2.d0*dibdc(k,l,3)*dibdc(m,n,3)+ 601 & 2.d0*(v3b-1.d0)*d2ibdc2(k,l,m,n,3))/coef 602 else 603 term= 604 & 2.d0*mm*(v3b-1.d0)**(2*mm-2)/coef* 605 & ((2*mm-1)*dibdc(k,l,3)*dibdc(m,n,3) 606 & +(v3b-1.d0)*d2ibdc2(k,l,m,n,3)) 607 endif 608 d2udc2(k,l,m,n)=d2udc2(k,l,m,n)+term 609 enddo 610 enddo 611 endif 612 endif 613! 614! Ogden form 615! 616 if((kode.lt.-3).and.(kode.gt.-7)) then 617 if(kode.eq.-4) then 618 nelconst=1 619 elseif(kode.eq.-5) then 620 nelconst=2 621 elseif(kode.eq.-6) then 622 nelconst=3 623 endif 624! 625! first derivative of U w.r.t. c(k,l) 626! 627 do l=1,3 628 do k=1,l 629 dudc(k,l)=0.d0 630 enddo 631 enddo 632! 633 do m=1,nelconst 634 coefa=constant(2*m) 635 coefd=constant(2*nelconst+m) 636 coefm=constant(2*m-1) 637 do l=1,3 638 do k=1,l 639 term=0.d0 640 do i=1,3 641 term=term+alb(i)**(coefa-1.d0)*dlbdc(k,l,i) 642 enddo 643 dudc(k,l)=dudc(k,l)+2.d0*coefm/coefa 644 & *term+2.d0*m/coefd* 645 & (v3b-1.d0)**(2*m-1)*dibdc(k,l,3) 646 enddo 647 enddo 648 enddo 649! 650! tangent stiffness matrix 651! second derivative of U w.r.t. c(k,l) and c(m,n) 652! 653 if(icmd.ne.3) then 654 nt=0 655 do i=1,21 656 k=kk(nt+1) 657 l=kk(nt+2) 658 m=kk(nt+3) 659 n=kk(nt+4) 660 nt=nt+4 661 d2udc2(k,l,m,n)=0.d0 662 enddo 663 do mm=1,nelconst 664 coefa=constant(2*mm) 665 coefd=constant(2*nelconst+mm) 666 coefm=constant(2*mm-1) 667 nt=0 668 do jj=1,21 669 k=kk(nt+1) 670 l=kk(nt+2) 671 m=kk(nt+3) 672 n=kk(nt+4) 673 nt=nt+4 674 term=0.d0 675 do i=1,3 676 term=term+alb(i)**(coefa-2.d0)*dlbdc(k,l,i)* 677 & dlbdc(m,n,i) 678 enddo 679 term=term*(coefa-1.d0) 680 do i=1,3 681 term=term+alb(i)**(coefa-1.d0) 682 & *d2lbdc2(k,l,m,n,i) 683 enddo 684 term=term*2.d0*coefm/coefa 685 d2udc2(k,l,m,n)=d2udc2(k,l,m,n)+term+(2*mm)* 686 & (2*mm-1)/coefd*(v3b-1.d0)**(2*mm-2)* 687 & dibdc(k,l,3)*dibdc(m,n,3)+2*mm/coefd 688 & *(v3b-1.d0)**(2*mm-1)*d2ibdc2(k,l,m,n,3) 689 enddo 690 enddo 691 endif 692 endif 693! 694! Arruda-Boyce model 695! 696 if(kode.eq.-1) then 697 coef=constant(2) 698! 699! first derivative of U w.r.t. c(k,l) 700! 701 do l=1,3 702 do k=1,l 703 dudc(k,l)=constant(1)*(0.5d0+v1b/(10.d0* 704 & coef**2)+33.d0*v1b*v1b/(1050.d0* 705 & coef**4)+76.d0*v1b**3/(7000.d0* 706 & coef**6)+2595.d0*v1b**4/(673750.d0* 707 & coef**8))*dibdc(k,l,1)+(v3b-1.d0/v3b) 708 & *dibdc(k,l,3)/constant(3) 709 enddo 710 enddo 711! 712! tangent stiffness matrix 713! second derivative of U w.r.t. c(k,l) and c(m,n) 714! 715 if(icmd.ne.3) then 716 nt=0 717 do jj=1,21 718 k=kk(nt+1) 719 l=kk(nt+2) 720 m=kk(nt+3) 721 n=kk(nt+4) 722 nt=nt+4 723 d2udc2(k,l,m,n)=constant(1)*(1.d0/(10.d0* 724 & coef**2)+66.d0*v1b/(1050.d0*coef**4)+228.d0 725 & *v1b**2/(7000.d0*coef**6)+10380.d0*v1b**3/ 726 & (673750.d0*coef**8))*dibdc(k,l,1)*dibdc(m,n,1) 727 & +constant(1)*(0.5d0+v1b/(10.d0*coef**2) 728 & +33.d0*v1b**2/ 729 & (1050.d0*coef**4)+76.d0*v1b**3/(7000.d0*coef**6)+ 730 & 2595.d0*v1b**4/(673750.d0*coef**8)) 731 & *d2ibdc2(k,l,m,n,1) 732 & +(1.d0+1.d0/v3b**2)*dibdc(k,l,3)*dibdc(m,n,3)/ 733 & constant(3)+(v3b-1.d0/v3b)*d2ibdc2(k,l,m,n,3) 734 & /constant(3) 735 enddo 736 endif 737 endif 738! 739! elastomeric foam behavior 740! 741 if((kode.lt.-15).and.(kode.gt.-18)) then 742 if(kode.eq.-15) then 743 nelconst=1 744 elseif(kode.eq.-16) then 745 nelconst=2 746 elseif(kode.eq.-17) then 747 nelconst=3 748 endif 749! 750! first derivative of U w.r.t. c(k,l) 751! 752 do l=1,3 753 do k=1,l 754 dudc(k,l)=0.d0 755 enddo 756 enddo 757! 758 do m=1,nelconst 759 coefa=constant(2*m) 760 coefb=constant(2*nelconst+m)/(1.d0-2.d0 761 & *constant(2*nelconst+m)) 762 coefm=constant(2*m-1) 763 do l=1,3 764 do k=1,l 765 term=0.d0 766 do i=1,3 767 term=term+al(i)**(coefa-1.d0)*dldc(k,l,i) 768 enddo 769 dudc(k,l)=dudc(k,l)+2.d0*coefm/coefa 770 & *(term-v3b**(-coefa*coefb-1.d0)* 771 & dibdc(k,l,3)) 772 enddo 773 enddo 774 enddo 775! 776! tangent stiffness matrix 777! second derivative of U w.r.t. c(k,l) and c(m,n) 778! 779 if(icmd.ne.3) then 780 nt=0 781 do i=1,21 782 k=kk(nt+1) 783 l=kk(nt+2) 784 m=kk(nt+3) 785 n=kk(nt+4) 786 nt=nt+4 787 d2udc2(k,l,m,n)=0.d0 788 enddo 789 do mm=1,nelconst 790 coefa=constant(2*mm) 791 coefb=constant(2*nelconst+mm)/(1.d0-2.d0 792 & *constant(2*nelconst+mm)) 793 coefm=constant(2*mm-1) 794 nt=0 795 do jj=1,21 796 k=kk(nt+1) 797 l=kk(nt+2) 798 m=kk(nt+3) 799 n=kk(nt+4) 800 nt=nt+4 801 term=0.d0 802 do i=1,3 803 term=term+(coefa-1.d0)*al(i)**(coefa-2.d0) 804 & *dldc(k,l,i)*dldc(m,n,i) 805 & +al(i)**(coefa-1.d0)*d2ldc2(k,l,m,n,i) 806 enddo 807 d2udc2(k,l,m,n)=d2udc2(k,l,m,n) 808 & +2.d0*coefm/ 809 & coefa*(term+(coefa*coefb+1.d0)*v3b 810 & **(-coefa*coefb-2.d0)*dibdc(k,l,3) 811 & *dibdc(m,n,3)-v3b**(-coefa*coefb-1.d0) 812 & *d2ibdc2(k,l,m,n,3)) 813 enddo 814 enddo 815 endif 816 endif 817! 818! storing the stiffness matrix and/or the stress 819! 820 if(icmd.ne.3) then 821! 822! storing the stiffness matrix 823! 824 nt=0 825 do i=1,21 826 k=kk(nt+1) 827 l=kk(nt+2) 828 m=kk(nt+3) 829 n=kk(nt+4) 830 nt=nt+4 831 elas(i)=4.d0*d2udc2(k,l,m,n) 832 enddo 833 endif 834! 835! store the stress at mechanical strain 836! 837 stre(1)=2.d0*dudc(1,1) 838 stre(2)=2.d0*dudc(2,2) 839 stre(3)=2.d0*dudc(3,3) 840 stre(4)=2.d0*dudc(1,2) 841 stre(5)=2.d0*dudc(1,3) 842 stre(6)=2.d0*dudc(2,3) 843! 844 enddo 845! 846 return 847 end 848