FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
11contains
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
424
425end 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 calupdateelasticarrudaboyce(matl, secttype, dstrain, dstress)
This subroutine provides to update stress and strain for Arrude-Royce 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 calelasticmooneyrivlin(matl, secttype, cijkl, strain)
This subroutine provides elastic tangent coefficient for Mooney-Rivlin hyperelastic material.
subroutine calelasticarrudaboyce(matl, secttype, cijkl, strain)
This subroutine provides elastic tangent coefficient for Arruda-Boyce hyperelastic material.
subroutine calelasticmooneyrivlinaniso(matl, secttype, cijkl, strain, cdsys)
This subroutine provides elastic tangent coefficient for anisotropic Mooney-Rivlin hyperelastic mater...
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 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.
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