FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
Hyperelastic.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
7  use hecmw_util
8  use mmaterial
9  implicit none
10 
11 contains
12 
14  subroutine cderiv( matl, sectType, ctn, itn, &
15  inv1b, inv2b, inv3b, dibdc, d2ibdc2, strain, direction, inv4b, dibdc_ani, d2ibdc2_ani )
16  type( tmaterial ), intent(in) :: matl
17  integer, intent(in) :: sectType
18  real(kind=kreal), intent(out) :: inv1b
19  real(kind=kreal), intent(out) :: inv2b
20  real(kind=kreal), intent(out) :: inv3b
21  real(kind=kreal), intent(out) :: dibdc(3,3,3)
22  real(kind=kreal), intent(out) :: d2ibdc2(3,3,3,3,3)
23  real(kind=kreal), intent(in) :: strain(6)
24  real(kind=kreal), intent(out) :: ctn(3,3)
25  real(kind=kreal), intent(out) :: itn(3,3)
26  real(kind=kreal), intent(in), optional :: direction(3)
27  real(kind=kreal), intent(out), optional :: inv4b
28  real(kind=kreal), intent(out), optional :: dibdc_ani(3,3)
29  real(kind=kreal), intent(out), optional :: d2ibdc2_ani(3,3,3,3)
30 
31  integer :: i, j, k, l, m, n
32 
33  real(kind=kreal) :: inv1, inv2, inv3, inv33
34  real(kind=kreal) :: delta(3,3)
35  real(kind=kreal) :: didc(3,3,3), ctninv(3,3)
36  real(kind=kreal) :: d2idc2(3,3,3,3,3)
37 
38  ! Kronecker's delta
39  delta(:,:) = 0.d0
40  delta(1,1) = 1.d0
41  delta(2,2) = 1.d0
42  delta(3,3) = 1.d0
43 
44  ! Green-Lagrange strain to right Cauchy-Green deformation tensor
45  ctn(1,1)=strain(1)*2.d0+1.d0
46  ctn(2,2)=strain(2)*2.d0+1.d0
47  ctn(3,3)=strain(3)*2.d0+1.d0
48  ctn(1,2)=strain(4); ctn(2,1)=ctn(1,2)
49  ctn(2,3)=strain(5); ctn(3,2)=ctn(2,3)
50  ctn(3,1)=strain(6); ctn(1,3)=ctn(3,1)
51 
52  itn(:,:) = delta(:,:)
53 
54  ! ----- calculate the invariant of C
55  inv1 = ctn(1,1)+ctn(2,2)+ctn(3,3)
56  inv2 = ctn(2,2)*ctn(3,3)+ctn(1,1)*ctn(3,3)+ctn(1,1)*ctn(2,2) &
57  -ctn(2,3)*ctn(2,3)-ctn(1,3)*ctn(1,3)-ctn(1,2)*ctn(1,2)
58  inv3 = ctn(1,1)*ctn(2,2)*ctn(3,3) &
59  +ctn(2,1)*ctn(3,2)*ctn(1,3) &
60  +ctn(3,1)*ctn(1,2)*ctn(2,3) &
61  -ctn(3,1)*ctn(2,2)*ctn(1,3) &
62  -ctn(2,1)*ctn(1,2)*ctn(3,3) &
63  -ctn(1,1)*ctn(3,2)*ctn(2,3)
64  inv33 = inv3**(-1.d0/3.d0)
65 
66  ! ----- calculate the inverse of C
67  ctninv(1,1)=(ctn(2,2)*ctn(3,3)-ctn(2,3)*ctn(2,3))/inv3
68  ctninv(2,2)=(ctn(1,1)*ctn(3,3)-ctn(1,3)*ctn(1,3))/inv3
69  ctninv(3,3)=(ctn(1,1)*ctn(2,2)-ctn(1,2)*ctn(1,2))/inv3
70  ctninv(1,2)=(ctn(1,3)*ctn(2,3)-ctn(1,2)*ctn(3,3))/inv3
71  ctninv(1,3)=(ctn(1,2)*ctn(2,3)-ctn(2,2)*ctn(1,3))/inv3
72  ctninv(2,3)=(ctn(1,2)*ctn(1,3)-ctn(1,1)*ctn(2,3))/inv3
73  ctninv(2,1)=ctninv(1,2)
74  ctninv(3,1)=ctninv(1,3)
75  ctninv(3,2)=ctninv(2,3)
76 
77  ! ----- calculate the derivative of the C-invariant with respect to c(i,j)
78  do i=1,3
79  do j=1,3
80  didc(i,j,1) = delta(i,j)
81  didc(i,j,2) = inv1*delta(i,j)-ctn(i,j)
82  didc(i,j,3) = inv3*ctninv(i,j)
83  enddo
84  enddo
85 
86  ! ----- calculate the secont derivative of the C-invariant
87  do k=1,3
88  do l=1,3
89  do m=1,3
90  do n=1,3
91  d2idc2(k,l,m,n,1)=0.d0
92  d2idc2(k,l,m,n,2)=delta(k,l)*delta(m,n)- &
93  (delta(k,m)*delta(l,n)+delta(k,n)*delta(l,m))/2.d0
94  d2idc2(k,l,m,n,3)=inv3*(ctninv(m,n)*ctninv(k,l)- &
95  (ctninv(k,m)*ctninv(n,l)+ctninv(k,n)*ctninv(m,l))/2.d0)
96  enddo
97  enddo
98  enddo
99  enddo
100 
101  ! ----- derivatives for the reduced invariants
102  inv1b = inv1*inv33
103  inv2b = inv2*inv33*inv33
104  inv3b = dsqrt(inv3)
105 
106  ! --- first derivative the reduced c-invarians w.r.t. c(i,j)
107  do i=1,3
108  do j=1,3
109  dibdc(i,j,1) = -inv33**4*inv1*didc(i,j,3)/3.d0 + inv33*didc(i,j,1)
110  dibdc(i,j,2) = -2.d0*inv33**5*inv2*didc(i,j,3)/3.d0 + inv33**2*didc(i,j,2)
111  dibdc(i,j,3) = didc(i,j,3)/(2.d0*dsqrt(inv3))
112  enddo
113  enddo
114 
115  ! --- second derivative of the reduced c-invariants w.r.t. c(i,j)
116  do i=1,3
117  do j=1,3
118  do k=1,3
119  do l=1,3
120  d2ibdc2(i,j,k,l,1) = 4.d0/9.d0*inv33**7*inv1*didc(i,j,3)*didc(k,l,3) &
121  -inv33**4/3.d0*(didc(k,l,1)*didc(i,j,3)+didc(i,j,1)*didc(k,l,3)) &
122  -inv33**4/3.d0*inv1*d2idc2(i,j,k,l,3) &
123  +inv33*d2idc2(i,j,k,l,1)
124  d2ibdc2(i,j,k,l,2) = 10.d0/9.d0*inv33**8*inv2*didc(i,j,3)*didc(k,l,3) &
125  -2.d0/3.d0*inv33**5*(didc(k,l,2)*didc(i,j,3)+didc(i,j,2)*didc(k,l,3)) &
126  -2.d0/3.d0*inv33**5*inv2*d2idc2(i,j,k,l,3) &
127  +inv33**2*d2idc2(i,j,k,l,2)
128  d2ibdc2(i,j,k,l,3) = -didc(i,j,3)*didc(k,l,3)/(4.d0*inv3**1.5d0) &
129  +d2idc2(i,j,k,l,3)/(2.d0*dsqrt(inv3))
130  enddo
131  enddo
132  enddo
133  enddo
134 
135  if( present(direction) ) call cderiv_aniso( ctn, inv3, didc(:,:,3), d2idc2(:,:,:,:,3), &
136  direction, inv4b, dibdc_ani, d2ibdc2_ani )
137 
138  end subroutine cderiv
139 
142  subroutine cderiv_aniso( ctn, inv3, didc_3, d2idc2_3, direction, inv4b, dibdc_ani, d2ibdc2_ani )
143  real(kind=kreal), intent(in) :: ctn(3,3)
144  real(kind=kreal), intent(in) :: inv3
145  real(kind=kreal), intent(in) :: didc_3(3,3)
146  real(kind=kreal), intent(in) :: d2idc2_3(3,3,3,3)
147  real(kind=kreal), intent(in) :: direction(3)
148  real(kind=kreal), intent(out) :: inv4b
149  real(kind=kreal), intent(out) :: dibdc_ani(3,3)
150  real(kind=kreal), intent(out) :: d2ibdc2_ani(3,3,3,3)
151 
152  integer :: i, j, k, l, m, n
153  real(kind=kreal) :: inv4, inv33, inv3m43, inv4d3
154  real(kind=kreal) :: d2ibdc24
155 
156  inv33 = inv3**(-1.d0/3.d0) ! = I_3^{-1/3}
157  inv3m43 = inv33/inv3 ! = I_3^{-4/3}
158 
159  inv4 = 0.d0
160  do i=1,3
161  do j=1,3
162  inv4 = inv4 + direction(j)*ctn(j,i)*direction(i)
163  enddo
164  enddo
165  inv4b = inv33*inv4 ! I4
166  inv4d3 = inv4/inv3 ! I4*I_3^{-1}
167 
168  ! --- first derivative the reduced 4th c-invarians w.r.t. c(i,j)
169  do i=1,3
170  do j=1,3
171  dibdc_ani(i,j) = inv33*( (-1.d0/3.d0)*didc_3(i,j)*inv4d3+direction(i)*direction(j) )
172  enddo
173  enddo
174 
175  ! --- second derivative of the reduced c-invariants w.r.t. c(i,j)
176  do k=1,3
177  do l=1,3
178  do m=1,3
179  do n=1,3
180  d2ibdc24 = (2.d0/3.d0)*didc_3(k,l)*didc_3(m,n)*inv4d3
181  d2ibdc24 = d2ibdc24 - d2idc2_3(k,l,m,n)*inv4
182  d2ibdc24 = d2ibdc24 - direction(k)*direction(l)*didc_3(m,n)
183  d2ibdc24 = d2ibdc24 - didc_3(k,l)*direction(m)*direction(n)
184  d2ibdc2_ani(k,l,m,n) = (1.d0/3.d0)*inv3m43*d2ibdc24
185  enddo
186  enddo
187  enddo
188  enddo
189 
190  end subroutine
191 
192 
193  !-------------------------------------------------------------------------------
195  !
196  !-------------------------------------------------------------------------------
197  subroutine calelasticarrudaboyce( matl, sectType, cijkl, strain )
198  type( tmaterial ), intent(in) :: matl
199  integer, intent(in) :: secttype
200  real(kind=kreal), intent(out) :: cijkl(3,3,3,3)
201  real(kind=kreal), intent(in) :: strain(6)
202 
203  integer :: i, j, k, l
204  real(kind=kreal) :: ctn(3,3), itn(3,3)
205  real(kind=kreal) :: inv1b, inv2b, inv3b
206  real(kind=kreal) :: dibdc(3,3,3)
207  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
208  real(kind=kreal) :: constant(3), coef
209 
210  call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
211  dibdc, d2ibdc2, strain )
212  constant(1:3)=matl%variables(m_plconst1:m_plconst3)
213  coef = constant(2)
214 
215  forall( i=1:3, j=1:3, k=1:3, l=1:3 )
216  cijkl(i,j,k,l) = constant(1)*(1.d0/(10.d0*coef**2) &
217  +66.d0*inv1b/(1050.d0*coef**4) &
218  +228.d0*inv1b**2/(7000.d0*coef**6) &
219  +10380.d0*inv1b**3/(673750.d0*coef**8)) &
220  *dibdc(i,j,1)*dibdc(k,l,1) &
221  +constant(1)*(0.5d0+inv1b/(10.d0*coef**2) &
222  +33.d0*inv1b**2/(1050.d0*coef**4) &
223  +76.d0*inv1b**3/(7000.d0*coef**6) &
224  +2595.d0*inv1b**4/(673750.d0*coef**8)) &
225  *d2ibdc2(i,j,k,l,1) &
226  +(1.d0+1.d0/inv3b**2)*dibdc(i,j,3)*dibdc(k,l,3)/constant(3) &
227  +(inv3b-1.d0/inv3b)*d2ibdc2(i,j,k,l,3)/constant(3)
228  end forall
229  cijkl(:,:,:,:) = 4.d0*cijkl(:,:,:,:)
230 
231  end subroutine calelasticarrudaboyce
232 
233  !-------------------------------------------------------------------------------
235  subroutine calupdateelasticarrudaboyce( matl, sectType, dstrain, dstress )
236  type( tmaterial ), intent(in) :: matl
237  integer, intent(in) :: secttype
238  real(kind=kreal), intent(out) :: dstress(6)
239  real(kind=kreal), intent(in) :: dstrain(6)
240 
241  integer :: i, j
242  real(kind=kreal) :: ctn(3,3), itn(3,3)
243  real(kind=kreal) :: inv1b, inv2b, inv3b
244  real(kind=kreal) :: dibdc(3,3,3)
245  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
246  real(kind=kreal) :: constant(3), coef
247  real(kind=kreal) :: pkstress(3,3)
248 
249 
250  constant(1:3)=matl%variables(m_plconst1:m_plconst3)
251  coef = constant(2)
252  call cderiv( matl, secttype, ctn, itn, &
253  inv1b, inv2b, inv3b, &
254  dibdc, d2ibdc2, dstrain )
255 
256 
257  ! ----- calculate the stress
258  do i=1,3
259  do j=1,3
260  pkstress(i,j) = constant(1)*( 0.5d0+inv1b/(10.d0*coef**2) &
261  +33.d0*inv1b*inv1b/(1050.d0*coef**4) &
262  +76.d0*inv1b**3/(7000.d0*coef**6) &
263  +2595.d0*inv1b**4/(673750.d0*coef**8))*dibdc(i,j,1) &
264  +(inv3b-1.d0/inv3b)*dibdc(i,j,3)/constant(3)
265  enddo
266  enddo
267 
268  dstress(1) = 2.d0*pkstress(1,1)
269  dstress(2) = 2.d0*pkstress(2,2)
270  dstress(3) = 2.d0*pkstress(3,3)
271  dstress(4) = pkstress(1,2) + pkstress(2,1)
272  dstress(5) = pkstress(2,3) + pkstress(3,2)
273  dstress(6) = pkstress(1,3) + pkstress(3,1)
274 
275  end subroutine calupdateelasticarrudaboyce
276 
277  !-------------------------------------------------------------------------------
279  !-------------------------------------------------------------------------------
280  subroutine calelasticmooneyrivlin( matl, sectType, cijkl, strain )
281  type( tmaterial ), intent(in) :: matl
282  integer, intent(in) :: secttype
283  real(kind=kreal), intent(out) :: cijkl(3,3,3,3)
284  real(kind=kreal), intent(in) :: strain(6)
285 
286  integer :: k, l, m, n
287  real(kind=kreal) :: ctn(3,3), itn(3,3)
288  real(kind=kreal) :: inv1b, inv2b, inv3b
289  real(kind=kreal) :: dibdc(3,3,3)
290  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
291  real(kind=kreal) :: constant(3)
292 
293 
294  constant(1:3)=matl%variables(m_plconst1:m_plconst3)
295  call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
296  dibdc, d2ibdc2, strain )
297 
298  forall( k=1:3, l=1:3, m=1:3, n=1:3 )
299  cijkl(k,l,m,n) = d2ibdc2(k,l,m,n,1)*constant(1) + &
300  d2ibdc2(k,l,m,n,2)*constant(2) + &
301  2.d0*(dibdc(k,l,3)*dibdc(m,n,3)+ &
302  (inv3b-1.d0)*d2ibdc2(k,l,m,n,3))/constant(3)
303  end forall
304  cijkl(:,:,:,:)=4.d0*cijkl(:,:,:,:)
305 
306  end subroutine calelasticmooneyrivlin
307 
308  !-------------------------------------------------------------------------------
310  !-------------------------------------------------------------------------------
311  subroutine calupdateelasticmooneyrivlin( matl, sectType, strain, stress )
312  type( tmaterial ), intent(in) :: matl
313  integer, intent(in) :: secttype
314  real(kind=kreal), intent(out) :: stress(6)
315  real(kind=kreal), intent(in) :: strain(6)
316 
317  integer :: k, l
318  real(kind=kreal) :: ctn(3,3), itn(3,3)
319  real(kind=kreal) :: inv1b, inv2b, inv3b
320  real(kind=kreal) :: dibdc(3,3,3)
321  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
322  real(kind=kreal) :: constant(3)
323  real(kind=kreal) :: dudc(3,3)
324 
325  constant(1:3)=matl%variables(m_plconst1:m_plconst3)
326  call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
327  dibdc, d2ibdc2, strain )
328 
329 
330  ! ----- stress
331  do l=1,3
332  do k=1,3
333  dudc(k,l) = dibdc(k,l,1)*constant(1)+dibdc(k,l,2)*constant(2) &
334  +2.d0*(inv3b-1.d0)*dibdc(k,l,3)/constant(3)
335  enddo
336  enddo
337 
338  stress(1)=2.d0*dudc(1,1)
339  stress(2)=2.d0*dudc(2,2)
340  stress(3)=2.d0*dudc(3,3)
341  stress(4)=2.d0*dudc(1,2)
342  stress(5)=2.d0*dudc(2,3)
343  stress(6)=2.d0*dudc(1,3)
344 
345  end subroutine calupdateelasticmooneyrivlin
346 
347  !-------------------------------------------------------------------------------
349  !-------------------------------------------------------------------------------
350  subroutine calelasticmooneyrivlinaniso( matl, sectType, cijkl, strain, cdsys )
351  type( tmaterial ), intent(in) :: matl
352  integer, intent(in) :: secttype
353  real(kind=kreal), intent(out) :: cijkl(3,3,3,3)
354  real(kind=kreal), intent(in) :: strain(6)
355  real(kind=kreal), intent(in) :: cdsys(3,3)
356 
357  integer :: i, j, k, l, m, n, jj
358  real(kind=kreal) :: ctn(3,3), itn(3,3)
359  real(kind=kreal) :: inv1b, inv2b, inv3b, inv4b
360  real(kind=kreal) :: dibdc(3,3,3)
361  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
362  real(kind=kreal) :: constant(5)
363  real(kind=kreal) :: dibdc_ani(3,3)
364  real(kind=kreal) :: d2ibdc2_ani(3,3,3,3)
365 
366 
367  constant(1:5)=matl%variables(m_plconst1:m_plconst5)
368  call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
369  dibdc, d2ibdc2, strain, cdsys(1,1:3), inv4b, dibdc_ani, d2ibdc2_ani )
370 
371  forall( k=1:3, l=1:3, m=1:3, n=1:3 )
372  cijkl(k,l,m,n) = d2ibdc2(k,l,m,n,1)*constant(1) + &
373  d2ibdc2(k,l,m,n,2)*constant(2) + &
374  2.d0*(dibdc(k,l,3)*dibdc(m,n,3)+ &
375  (inv3b-1.d0)*d2ibdc2(k,l,m,n,3))/constant(3)+ &
376  (2.d0*constant(4)+6.d0*(inv4b-1.d0)*constant(5))*dibdc_ani(k,l)*dibdc_ani(m,n)+ &
377  (inv4b-1.d0)*(2.d0*constant(4)+3.d0*(inv4b-1.d0)*constant(5))*d2ibdc2_ani(k,l,m,n)
378  end forall
379  cijkl(:,:,:,:)=4.d0*cijkl(:,:,:,:)
380 
381  end subroutine calelasticmooneyrivlinaniso
382 
383  !-------------------------------------------------------------------------------
385  !-------------------------------------------------------------------------------
386  subroutine calupdateelasticmooneyrivlinaniso( matl, sectType, strain, stress, cdsys )
387  type( tmaterial ), intent(in) :: matl
388  integer, intent(in) :: secttype
389  real(kind=kreal), intent(out) :: stress(6)
390  real(kind=kreal), intent(in) :: strain(6)
391  real(kind=kreal), intent(in) :: cdsys(3,3)
392 
393  integer :: k, l
394  real(kind=kreal) :: ctn(3,3), itn(3,3)
395  real(kind=kreal) :: inv1b, inv2b, inv3b, inv4b
396  real(kind=kreal) :: dibdc(3,3,3)
397  real(kind=kreal) :: d2ibdc2(3,3,3,3,3)
398  real(kind=kreal) :: constant(5)
399  real(kind=kreal) :: dudc(3,3)
400  real(kind=kreal) :: dibdc_ani(3,3)
401  real(kind=kreal) :: d2ibdc2_ani(3,3,3,3)
402 
403  constant(1:5)=matl%variables(m_plconst1:m_plconst5)
404  call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
405  dibdc, d2ibdc2, strain, cdsys(1,1:3), inv4b, dibdc_ani, d2ibdc2_ani )
406 
407  ! ----- stress
408  do l=1,3
409  do k=1,3
410  dudc(k,l) = dibdc(k,l,1)*constant(1)+dibdc(k,l,2)*constant(2) &
411  +2.d0*(inv3b-1.d0)*dibdc(k,l,3)/constant(3) &
412  +(inv4b-1.d0)*(2.d0*constant(4)+3.d0*(inv4b-1.d0)*constant(5))*dibdc_ani(k,l)
413  enddo
414  enddo
415 
416  stress(1)=2.d0*dudc(1,1)
417  stress(2)=2.d0*dudc(2,2)
418  stress(3)=2.d0*dudc(3,3)
419  stress(4)=2.d0*dudc(1,2)
420  stress(5)=2.d0*dudc(2,3)
421  stress(6)=2.d0*dudc(1,3)
422 
423  end subroutine calupdateelasticmooneyrivlinaniso
424 
425 end module
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
subroutine calelasticmooneyrivlin(matl, sectType, cijkl, strain)
This subroutine provides elastic tangent coefficient for Mooney-Rivlin hyperelastic material.
subroutine cderiv_aniso(ctn, inv3, didc_3, d2idc2_3, direction, inv4b, dibdc_ani, d2ibdc2_ani)
This subroutine calculates derivative of the 4th reduced invariant I4 with respect to Cauchy-Green te...
subroutine calupdateelasticmooneyrivlinaniso(matl, sectType, strain, stress, cdsys)
This subroutine provides to update stress and strain for anisotropic Mooney-Rivlin material.
subroutine calupdateelasticmooneyrivlin(matl, sectType, strain, stress)
This subroutine provides to update stress and strain for Mooney-Rivlin material.
subroutine calupdateelasticarrudaboyce(matl, sectType, dstrain, dstress)
This subroutine provides to update stress and strain for Arrude-Royce material.
subroutine calelasticarrudaboyce(matl, sectType, cijkl, strain)
This subroutine provides elastic tangent coefficient for Arruda-Boyce hyperelastic material.
subroutine cderiv(matl, sectType, ctn, itn, inv1b, inv2b, inv3b, dibdc, d2ibdc2, strain, direction, inv4b, dibdc_ani, d2ibdc2_ani)
This subroutine calculates derivative of the invariant with respect to Cauchy-Green tensor.
subroutine calelasticmooneyrivlinaniso(matl, sectType, cijkl, strain, cdsys)
This subroutine provides elastic tangent coefficient for anisotropic Mooney-Rivlin hyperelastic mater...
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter m_plconst5
Definition: material.f90:94
integer(kind=kint), parameter m_plconst1
Definition: material.f90:90
integer(kind=kint), parameter m_plconst3
Definition: material.f90:92
Stucture to management all material relates data.
Definition: material.f90:144