FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 
13 contains
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)
92  use m_elasticlinear
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)
183  use m_elasticlinear
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 
302 end 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 updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
This subroutine provides to update stress for viscoelastic material.
subroutine calviscoelasticmatrix(matl, sectType, dt, D, temp)
This subroutine calculates tangent moduli for isotropic 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