14 subroutine cderiv( matl, sectType, ctn, itn, &
15 inv1b, inv2b, inv3b, dibdc, d2ibdc2, strain, direction, inv4b, dibdc_ani, d2ibdc2_ani )
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)
31 integer :: i, j, k, l, m, n
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)
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)
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)
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)
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)
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)
103 inv2b = inv2*inv33*inv33
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))
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))
135 if(
present(direction) )
call cderiv_aniso( ctn, inv3, didc(:,:,3), d2idc2(:,:,:,:,3), &
136 direction, inv4b, dibdc_ani, d2ibdc2_ani )
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)
152 integer :: i, j, k, l, m, n
153 real(kind=
kreal) :: inv4, inv33, inv3m43, inv4d3
154 real(kind=
kreal) :: d2ibdc24
156 inv33 = inv3**(-1.d0/3.d0)
162 inv4 = inv4 + direction(j)*ctn(j,i)*direction(i)
171 dibdc_ani(i,j) = inv33*( (-1.d0/3.d0)*didc_3(i,j)*inv4d3+direction(i)*direction(j) )
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
199 integer,
intent(in) :: secttype
200 real(kind=
kreal),
intent(out) :: cijkl(3,3,3,3)
201 real(kind=
kreal),
intent(in) :: strain(6)
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
210 call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
211 dibdc, d2ibdc2, strain )
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)
229 cijkl(:,:,:,:) = 4.d0*cijkl(:,:,:,:)
237 integer,
intent(in) :: secttype
238 real(kind=
kreal),
intent(out) :: dstress(6)
239 real(kind=
kreal),
intent(in) :: dstrain(6)
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)
252 call cderiv( matl, secttype, ctn, itn, &
253 inv1b, inv2b, inv3b, &
254 dibdc, d2ibdc2, dstrain )
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)
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)
282 integer,
intent(in) :: secttype
283 real(kind=
kreal),
intent(out) :: cijkl(3,3,3,3)
284 real(kind=
kreal),
intent(in) :: strain(6)
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)
295 call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
296 dibdc, d2ibdc2, strain )
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)
304 cijkl(:,:,:,:)=4.d0*cijkl(:,:,:,:)
313 integer,
intent(in) :: secttype
314 real(kind=
kreal),
intent(out) :: stress(6)
315 real(kind=
kreal),
intent(in) :: strain(6)
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)
326 call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
327 dibdc, d2ibdc2, strain )
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)
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)
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)
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)
368 call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
369 dibdc, d2ibdc2, strain, cdsys(1,1:3), inv4b, dibdc_ani, d2ibdc2_ani )
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)
379 cijkl(:,:,:,:)=4.d0*cijkl(:,:,:,:)
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)
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)
404 call cderiv( matl, secttype, ctn, itn, inv1b, inv2b, inv3b, &
405 dibdc, d2ibdc2, strain, cdsys(1,1:3), inv4b, dibdc_ani, d2ibdc2_ani )
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)
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)
integer(kind=4), parameter kreal
This module provides functions for hyperelastic calculation.
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.
integer(kind=kint), parameter m_plconst5
integer(kind=kint), parameter m_plconst1
integer(kind=kint), parameter m_plconst3
Stucture to management all material relates data.