CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
rubber.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine rubber (elconloc, elas, emec, kode, didc, d2idc2, dibdc, d2ibdc2, dudc, d2udc2, dldc, d2ldc2, dlbdc, d2lbdc2, ithermal, icmd, beta, stre)
 

Function/Subroutine Documentation

◆ rubber()

subroutine rubber ( real*8, dimension(*)  elconloc,
real*8, dimension(*)  elas,
real*8, dimension(*)  emec,
integer  kode,
real*8, dimension(3,3,3)  didc,
real*8, dimension(3,3,3,3,3)  d2idc2,
real*8, dimension(3,3,3)  dibdc,
real*8, dimension(3,3,3,3,3)  d2ibdc2,
real*8, dimension(3,3)  dudc,
real*8, dimension(3,3,3,3)  d2udc2,
real*8, dimension(3,3,3)  dldc,
real*8, dimension(3,3,3,3,3)  d2ldc2,
real*8, dimension(3,3,3)  dlbdc,
real*8, dimension(3,3,3,3,3)  d2lbdc2,
integer  ithermal,
integer  icmd,
real*8, dimension(*)  beta,
real*8, dimension(*)  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 !
72 c 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.5)+
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
v1
static double * v1
Definition: mafillsmmain_se.c:40
Hosted by OpenAircraft.com, (Michigan UAV, LLC)