-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathNeoHookeanNumeric.f90
183 lines (167 loc) · 7.02 KB
/
NeoHookeanNumeric.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
module umatutils
!!! Utility module to be used for UMATs
implicit none
private
public dp, delta, m31tensorprod, m33det, m33tensorprod, mapnotation
integer, parameter :: dp=kind(0.d0)
real(dp), parameter :: delta(3, 3)=reshape([1, 0, 0, 0, 1, 0, 0, 0, 1],&
[3, 3])
contains
subroutine mapnotation(sigma, ccj, ntens, stress, ddsdde)
!! Map the matrix notation to vector notation
real(dp), intent(in) :: sigma(3, 3), ccj(3, 3, 3, 3)
integer, intent(in) :: ntens
real(dp), intent(out) :: stress(ntens), ddsdde(ntens, ntens)
integer, parameter :: notationmap(6, 2)=reshape([1, 2, 3, 1, 1, 2, 1,&
2, 3, 2, 3, 3], [6, 2])
integer :: k1, k2
do k1 = 1, 6
stress(k1) = sigma(notationmap(k1, 1), notationmap(k1, 2))
do k2 = 1, 6
ddsdde(k1, k2) = ccj(notationmap(k1, 1), notationmap(k1, 2),&
notationmap(k2, 1), notationmap(k2, 2))
end do
end do
end subroutine mapnotation
function m31tensorprod(a, b)
!! Return the tensor product from two 3x1 vectors
real(dp), dimension(3), intent(in) :: a, b
real(dp) :: m31tensorprod(3, 3)
integer :: k1, k2
forall(k1=1:3, k2=1:3)
m31tensorprod(k1, k2) = a(k1) * b(k2)
end forall
return
end function m31tensorprod
function m33tensorprod(a, b)
!! Return the tensor product from two 3x3 matrices
real(dp), dimension(3, 3), intent(in) :: a, b
real(dp) :: m33tensorprod(3, 3, 3, 3)
integer :: k1, k2, k3, k4
forall(k1=1:3, k2=1:3, k3=1:3, k4=1:3)
m33tensorprod(k1, k2, k3, k4) = a(k1, k2) * b(k3, k4)
end forall
return
end function m33tensorprod
real(dp) function m33det(m33)
!! Return determinant of a 3 by 3 matrix
real(dp), intent(in) :: m33(3, 3)
m33det = m33(1, 1) * m33(2, 2) * m33(3, 3)&
- m33(1, 1) * m33(2, 3) * m33(3, 2)&
- m33(1, 2) * m33(2, 1) * m33(3, 3)&
+ m33(1, 2) * m33(2, 3) * m33(3, 1)&
+ m33(1, 3) * m33(2, 1) * m33(3, 2)&
- m33(1, 3) * m33(2, 2) * m33(3, 1)
return
end function m33det
end module umatutils
module numerichyper
!!! Numeric method for hyperelastic materials
!!! Module to be embedded in the UMAT
!!! Paramters
!!! ---------
!!! det : Jacobian
!!! dfgrd : deformation gradient
!!! fbar : modified deformation gradient
!!! rcg : right Cauchy-Green tensor
!!! rcgbar : modified right Cauchy-Green tensor
!!! ibar1 : first invariant of rcgbar
!!! sigma : Cauchy stress
!!! ccj : tangent modulus for Cauchy stress in Jaumann rate
use umatutils, only: dp, delta, m33det, m31tensorprod, mapnotation
implicit none
private
public update_stress_ddsdde
contains
function getpsi(rcg, props) result(psi)
!! Return strain energy density given C and material properties
real(dp), intent(in) :: rcg(3, 3), props(:)
real(dp) :: rcgbar(3, 3), det, ibar1, psi
real(dp) :: c10, d1 ! Specific to Neo-Hookean model
integer :: k1
c10 = props(1)
d1 = props(2)
det = sqrt(m33det(rcg))
rcgbar = det**(-2._dp/3) * rcg
ibar1 = sum([(rcgbar(k1, k1), k1=1, 3)])
psi = c10 * (ibar1 - 3) + (det - 1)**2 / d1
return
end function getpsi
subroutine update_stress_ddsdde(props, dfgrd, ntens, stress, ddsdde)
!! Update the stress and ddsdde for Neo-Hookean
real(dp), intent(in) :: props(:), dfgrd(3, 3)
integer, intent(in) :: ntens
real(dp), intent(out) :: stress(ntens), ddsdde(ntens, ntens)
real(dp), parameter :: eps_s=1d-4, eps_c=1d-6
real(dp) :: sigma(3, 3), ccj(3, 3, 3, 3)
call get_sigma_ccj(dfgrd, props, eps_c, eps_s, sigma, ccj)
call mapnotation(sigma, ccj, ntens, stress, ddsdde)
end subroutine update_stress_ddsdde
function gettau(dfgrd, props, eps_s) result(tau)
!! Return the Cauchy stress given deformation gradient
real(dp), intent(in) :: dfgrd(3, 3), props(:), eps_s
real(dp) :: tau(3, 3), rcg(3, 3), psi, rcgptb(3, 3), psiptb, pk2(3, 3)
integer :: k1, k2
rcg = matmul(transpose(dfgrd), dfgrd)
psi = getpsi(rcg, props)
! Outer two layers are indices for the 2nd PK stress
do k1 = 1, 3
do k2 = 1, 3
rcgptb = rcg + eps_s * (&
m31tensorprod(delta(:, k1), delta(:, k2)) +&
m31tensorprod(delta(:, k2), delta(:, k1)))
psiptb = getpsi(rcgptb, props)
pk2(k1, k2) = (psiptb - psi) / eps_s
end do
end do
! Calculate Kirchoff stress from 2nd PK stress
tau = matmul(matmul(dfgrd, pk2), transpose(dfgrd))
return
end function gettau
subroutine get_sigma_ccj(dfgrd, props, eps_c, eps_s, sigma, ccj)
real(dp), intent(in) :: dfgrd(3, 3), props(:), eps_c, eps_s
real(dp), intent(out) :: sigma(3, 3), ccj(3, 3, 3, 3)
real(dp) :: tau(3, 3), tauptb(3, 3), fptb(3, 3), det
integer :: k3, k4
det = m33det(dfgrd)
tau = gettau(dfgrd, props, eps_s)
sigma = tau / det ! Pass out sigma
! Use k3 & k4 rather than k1 & k2 to denote that it's the k, l
! component of the elasticity tensor being calculated
do k3 = 1, 3
do k4 = 1, 3
fptb = dfgrd + eps_c / 2 * (&
matmul(m31tensorprod(delta(:, k3), delta(:, k4)), dfgrd) +&
matmul(m31tensorprod(delta(:, k4), delta(:, k3)), dfgrd))
tauptb = gettau(fptb, props, eps_s)
ccj(:, :, k3, k4) = (tauptb - tau) / eps_c / det
end do
end do
end subroutine get_sigma_ccj
end module numerichyper
subroutine umat(stress,statev,ddsdde,sse,spd,scd,&
rpl,ddsddt,drplde,drpldt,&
stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,&
ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,&
celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)
use numerichyper, only: update_stress_ddsdde
use umatutils, only: dp
implicit none
! This is a hack. The content of the 'aba_param.inc' is simply the
! Following two lines. I commented the first one, and modified the
! second one.
! include 'aba_param.inc'
! implicit real*8(a-h,o-z)
! parameter (nprecd=2)
integer, parameter :: nprecd=2
character*80, intent(in) :: cmname
real(dp), intent(in) :: stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),&
props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3), dtime, temp, &
dtemp, celent
integer, intent(in) :: ndi, nshr, ntens, nstatv, nprops, noel, npt, layer,&
kspt, kstep, kinc
real(dp), intent(inout) :: stress(ntens), statev(nstatv), sse, spd, scd,&
rpl, ddsdde(ntens, ntens), ddsddt(ntens), drplde(ntens), drpldt, pnewdt
real(dp) :: c10, d1
call update_stress_ddsdde(props, dfgrd1, ntens, stress, ddsdde)
end subroutine umat