FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
Viscoelastic.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 private :: hvisc, trs, trsinc
12
13contains
14
16 real(kind=kreal) function hvisc(x,expx)
17 real(kind=kreal), intent(in) :: x
18 real(kind=kreal), intent(in) :: expx
19
20 if(x < 1.d-04) then
21
22 hvisc = 1.d0 - 0.5d0*x*(1.d0 - x/3.d0*(1.d0 &
23 - 0.25d0*x*(1.d0 - 0.2d0*x)))
24
25 else
26
27 hvisc = (1.d0 - expx)/x
28
29 endif
30
31 end function
32
34 real(kind=kreal) function trsinc(tn1,tn,mtype, mvar)
35 real(kind=kreal), intent(in) :: tn1
36 real(kind=kreal), intent(in) :: tn
37 integer, intent(in) :: mtype
38 real(kind=kreal), intent(in) :: mvar(:)
39
40 real(kind=kreal) :: hsn1, hsn, asn1, asn
41
42 if(mtype==viscoelastic+2) then ! Arrhenius
43
44 hsn = -mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
45 hsn1 = -mvar(2)*( 1.d0/(tn1-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
46
47 else ! WLF
48
49 if( mvar(3)+tn1-mvar(1)<=0.d0 ) then
50 trsinc= -1.d0; return
51 endif
52 if( mvar(3)+tn-mvar(1)<=0.d0 ) then
53 trsinc= -1.d0; return
54 endif
55 hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
56 hsn1 = mvar(2)*( tn1-mvar(1) )/ (mvar(3)+tn1-mvar(1) )*dlog(10.d0)
57
58 endif
59
60 if( dabs(hsn1-hsn)<1.d-5 ) then
61 trsinc = dexp((hsn1+hsn)*0.5d0)
62 else
63 asn1 = dexp(hsn1)
64 asn = dexp(hsn)
65 trsinc = (asn1-asn)/(hsn1-hsn)
66 endif
67
68 end function
69
71 real(kind=kreal) function trs(tn,mtype, mvar)
72 real(kind=kreal), intent(in) :: tn
73 integer, intent(in) :: mtype
74 real(kind=kreal), intent(in) :: mvar(:)
75
76 real(kind=kreal) :: hsn
77
78 if(mtype==viscoelastic+2) then ! Arrhenius
79 hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
80 else ! WLF
81 hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
82 endif
83
84 trs = dexp(hsn)
85
86 end function
87
88 !-------------------------------------------------------------------------------
90 !-------------------------------------------------------------------------------
91 subroutine calviscoelasticmatrix(matl, sectType, dt, D, temp)
93 type( tmaterial ), intent(in) :: matl
94 integer, intent(in) :: secttype
95 real(kind=kreal), intent(in) :: dt
96 real(kind=kreal), intent(out) :: d(:,:)
97 real(kind=kreal), optional :: temp
98
99 integer i,j, n
100 real(kind=kreal) :: g,gg,k,kg, gfac,exp_n,mu_0,mu_n,dq_n,dtau
101 real(kind=kreal) :: ina(1), outa(2), ee, pp, ddt
102
103 type(ttable), pointer :: dicval
104 logical :: ierr
105
106
107 d(:,:)=0.d0
108 if( dt==0.d0 ) then
109 if( present( temp ) ) then
110 call calelasticmatrix( matl, d3, d, temp )
111 else
112 call calelasticmatrix( matl, d3, d )
113 endif
114 return
115 endif
116
117 ddt = dt
118 if( present(temp) ) then
119 ina(1) = temp
120 call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
121 if( ierr ) then
122 ee = matl%variables(m_youngs)
123 pp = matl%variables(m_poisson)
124 else
125 ee = outa(1)
126 pp = outa(2)
127 endif
128 if( matl%mtype>viscoelastic ) ddt=trs(temp,matl%mtype, matl%variables)*dt
129 else
130 ee = matl%variables(m_youngs)
131 pp = matl%variables(m_poisson)
132 endif
133
134 ! Set elastic parameters for G (mu) and lambda
135
136 g = ee/(2.d0*(1.d0 + pp))
137 k = ee/(3.d0*(1.d0 - 2.d0*pp))
138
139
140
141 ! Set properties for integrating the q_i terms
142
143 mu_0 = 0.0d0
144 gfac = 0.0d0
145
146 call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
147 if( ierr ) stop "Viscoelastic properties not defined"
148
149 do n = 1,dicval%tbrow
150 mu_n = dicval%tbval(1,n)
151 dtau = ddt/dicval%tbval(2,n)
152 exp_n = dexp(-dtau)
153
154 dq_n = mu_n * hvisc(dtau,exp_n)
155 gfac = gfac + dq_n
156 mu_0 = mu_0 + mu_n
157
158 end do
159
160 mu_0 = 1.d0 - mu_0
161 gfac = gfac + mu_0
162
163 ! Set tangent parameters
164 gg = g*gfac
165 kg = k - 0.6666666666666666d0*gg
166 do j =1,3
167 do i = 1,3
168 d(i,j) = kg
169 end do
170 d(j,j) = d(j,j) + 2.d0*gg
171 end do
172
173 do i = 4,6
174 d(i,i) = gg
175 end do
176
177
178 end subroutine
179
180 !-------------------------------------------------------------------------------
182 subroutine updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
184 type( tmaterial ), intent(in) :: matl
185 integer, intent(in) :: sectType
186 real(kind=kreal), intent(in) :: eps(6)
187 real(kind=kreal), intent(out) :: sig(6)
188 real(kind=kreal), intent(inout) :: vsig(:)
189 real(kind=kreal), intent(in) :: dt
190 real(kind=kreal), optional :: temp
191 real(kind=kreal), optional :: tempn
192
193 integer i,n
194 real(kind=kreal) :: g,k,kth, exp_n,mu_0,mu_n,dq_n,dtau, theta
195 real(kind=kreal) :: ina(1), outa(2), ee, pp, ddt
196
197 real(kind=kreal) :: devstrain(6), en(6), d(6,6)
198 type(ttable), pointer :: dicval
199 logical :: ierr
200
201 ddt = dt
202 if( present(temp) ) then
203 ina(1) = temp
204 call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
205 if( ierr ) then
206 ee = matl%variables(m_youngs)
207 pp = matl%variables(m_poisson)
208 else
209 ee = outa(1)
210 pp = outa(2)
211 endif
212 if( matl%mtype>viscoelastic ) ddt=trsinc(temp,tempn, matl%mtype, matl%variables)*dt
213 if( ddt<=0.d0 ) then
214 call calelasticmatrix( matl, secttype, d, temp )
215 sig = matmul( d(:,:), eps(:) )
216 vsig = 0.d0
217 return
218 endif
219 else
220 ee = matl%variables(m_youngs)
221 pp = matl%variables(m_poisson)
222 endif
223
224 ! Set elastic parameters for G (mu) and lambda
225
226 g = ee/(2.d0*(1.d0 + pp))
227 k = ee/(3.d0*(1.d0 - 2.d0*pp))
228
229 ! Compute volumetric strain and deviatoric components
230
231 theta = (eps(1) + eps(2) + eps(3))/3.d0
232 do i = 1,3
233 devstrain(i) = eps(i) - theta
234 end do
235 do i = 4,6
236 devstrain(i) = eps(i)*0.5d0
237 end do
238
239 ! Set properties for integrating the q_i terms
240
241 sig(:) = 0.0d0
242 mu_0 = 0.0d0
243
244 call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
245 if( ierr ) stop "Viscoelastic properties not defined"
246
247 en(:) = vsig(12*dicval%tbrow+1:12*dicval%tbrow+6)
248
249 do n = 1,dicval%tbrow
250 mu_n = dicval%tbval(1,n)
251 dtau = ddt/dicval%tbval(2,n)
252 exp_n = dexp(-dtau)
253
254 dq_n = mu_n * hvisc(dtau,exp_n)
255 mu_0 = mu_0 + mu_n
256
257 do i = 1,6
258 vsig((n-1)*12+i+6) = exp_n*vsig((n-1)*12+i) + dq_n*(devstrain(i) - en(i))
259 sig(i) = sig(i) + vsig((n-1)*12+i+6)
260 end do
261 end do
262
263 ! Finish stress update
264
265 mu_0 = 1.d0 - mu_0
266 do i = 1,6
267 sig(i) = 2.d0*g*(mu_0*devstrain(i) + sig(i))
268 end do
269
270 ! Add elastic bulk term
271
272 kth = k*theta*3.0d0
273 do i = 1,3
274 sig(i) = sig(i) + kth
275 end do
276
277 end subroutine
278
280 subroutine updateviscoelasticstate( gauss )
281 use mmechgauss
282 type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
283
284 integer :: i, n, nrow
285 real(kind=kreal) :: thetan, vstrain(6)
286 nrow = fetch_tablerow( mc_viscoelastic, gauss%pMaterial%dict )
287 do n = 1,nrow
288 do i = 1,6
289 gauss%fstatus((n-1)*12+i) = gauss%fstatus((n-1)*12+i+6)
290 end do
291 end do
292 vstrain = gauss%strain
293 thetan = (vstrain(1) + vstrain(2) + vstrain(3))/3.d0
294 do i = 1,3
295 gauss%fstatus(nrow*12+i) = vstrain(i) - thetan
296 end do
297 do i = 4,6
298 gauss%fstatus(nrow*12+i) = vstrain(i)*0.5d0
299 end do
300 end subroutine
301
302end module
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
This module provides functions for elastic material.
subroutine calelasticmatrix(matl, secttype, d, temp)
Calculate isotropic elastic matrix.
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter m_youngs
Definition: material.f90:84
character(len=dict_key_length) mc_viscoelastic
Definition: material.f90:124
integer(kind=kint), parameter viscoelastic
Definition: material.f90:70
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter m_poisson
Definition: material.f90:85
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:119
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.
subroutine calviscoelasticmatrix(matl, secttype, dt, d, temp)
This subroutine calculates tangent moduli for isotropic viscoelastic material.
subroutine updateviscoelastic(matl, secttype, eps, sig, vsig, dt, temp, tempn)
This subroutine provides to update stress for viscoelastic material.
Stucture to management all material relates data.
Definition: material.f90:144
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13