FrontISTR  5.7.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 !-------------------------------------------------------------------------------
6 module mviscoelastic
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), intent(in) :: 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  d(:,:)=0.d0
107 
108  ddt = dt
109 
110  ina(1) = temp
111  call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
112  if( ierr ) then
113  ee = matl%variables(m_youngs)
114  pp = matl%variables(m_poisson)
115  else
116  ee = outa(1)
117  pp = outa(2)
118  endif
119  if( matl%mtype>viscoelastic ) ddt=trs(temp,matl%mtype, matl%variables)*dt
120 
121  if( ddt<=0.d0 ) then
122  call calelasticmatrix( matl, d3, d, temp )
123  return
124  endif
125 
126  ! Set elastic parameters for G (mu) and lambda
127 
128  g = ee/(2.d0*(1.d0 + pp))
129  k = ee/(3.d0*(1.d0 - 2.d0*pp))
130 
131  ! Set properties for integrating the q_i terms
132 
133  mu_0 = 0.0d0
134  gfac = 0.0d0
135 
136  call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
137  if( ierr ) stop "Viscoelastic properties not defined"
138 
139  do n = 1,dicval%tbrow
140  mu_n = dicval%tbval(1,n)
141  dtau = ddt/dicval%tbval(2,n)
142  exp_n = dexp(-dtau)
143 
144  dq_n = mu_n * hvisc(dtau,exp_n)
145  gfac = gfac + dq_n
146  mu_0 = mu_0 + mu_n
147 
148  end do
149 
150  mu_0 = 1.d0 - mu_0
151  gfac = gfac + mu_0
152 
153  ! Set tangent parameters
154  gg = g*gfac
155  kg = k - 0.6666666666666666d0*gg
156  do j =1,3
157  do i = 1,3
158  d(i,j) = kg
159  end do
160  d(j,j) = d(j,j) + 2.d0*gg
161  end do
162 
163  do i = 4,6
164  d(i,i) = gg
165  end do
166 
167 
168  end subroutine
169 
170  !-------------------------------------------------------------------------------
172  subroutine updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
173  use m_elasticlinear
174  type( tmaterial ), intent(in) :: matl
175  integer, intent(in) :: sectType
176  real(kind=kreal), intent(in) :: eps(6)
177  real(kind=kreal), intent(out) :: sig(6)
178  real(kind=kreal), intent(inout) :: vsig(:)
179  real(kind=kreal), intent(in) :: dt
180  real(kind=kreal), intent(in) :: temp
181  real(kind=kreal), intent(in) :: tempn
182 
183  integer i,n
184  real(kind=kreal) :: g,k,kth, exp_n,mu_0,mu_n,dq_n,dtau, theta
185  real(kind=kreal) :: ina(1), outa(2), ee, pp, ddt
186 
187  real(kind=kreal) :: devstrain(6), en(6), d(6,6)
188  type(ttable), pointer :: dicval
189  logical :: ierr
190 
191  ddt = dt
192 
193  ina(1) = temp
194  call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
195  if( ierr ) then
196  ee = matl%variables(m_youngs)
197  pp = matl%variables(m_poisson)
198  else
199  ee = outa(1)
200  pp = outa(2)
201  endif
202  if( matl%mtype>viscoelastic ) ddt=trsinc(temp, tempn, matl%mtype, matl%variables)*dt
203 
204  if( ddt<=0.d0 ) then
205  call calelasticmatrix( matl, secttype, d, temp )
206  sig = matmul( d(:,:), eps(:) )
207  vsig = 0.d0
208  return
209  endif
210 
211  ! Set elastic parameters for G (mu) and lambda
212 
213  g = ee/(2.d0*(1.d0 + pp))
214  k = ee/(3.d0*(1.d0 - 2.d0*pp))
215 
216  ! Compute volumetric strain and deviatoric components
217 
218  theta = (eps(1) + eps(2) + eps(3))/3.d0
219  do i = 1,3
220  devstrain(i) = eps(i) - theta
221  end do
222  do i = 4,6
223  devstrain(i) = eps(i)*0.5d0
224  end do
225 
226  ! Set properties for integrating the q_i terms
227 
228  sig(:) = 0.0d0
229  mu_0 = 0.0d0
230 
231  call fetch_table( mc_viscoelastic, matl%dict, dicval, ierr )
232  if( ierr ) stop "Viscoelastic properties not defined"
233 
234  en(:) = vsig(12*dicval%tbrow+1:12*dicval%tbrow+6)
235 
236  do n = 1,dicval%tbrow
237  mu_n = dicval%tbval(1,n)
238  dtau = ddt/dicval%tbval(2,n)
239  exp_n = dexp(-dtau)
240 
241  dq_n = mu_n * hvisc(dtau,exp_n)
242  mu_0 = mu_0 + mu_n
243 
244  do i = 1,6
245  vsig((n-1)*12+i+6) = exp_n*vsig((n-1)*12+i) + dq_n*(devstrain(i) - en(i))
246  sig(i) = sig(i) + vsig((n-1)*12+i+6)
247  end do
248  end do
249 
250  ! Finish stress update
251 
252  mu_0 = 1.d0 - mu_0
253  do i = 1,6
254  sig(i) = 2.d0*g*(mu_0*devstrain(i) + sig(i))
255  end do
256 
257  ! Add elastic bulk term
258 
259  kth = k*theta*3.0d0
260  do i = 1,3
261  sig(i) = sig(i) + kth
262  end do
263 
264  end subroutine
265 
267  subroutine updateviscoelasticstate( gauss )
268  use mmechgauss
269  type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
270 
271  integer :: i, n, nrow
272  real(kind=kreal) :: thetan, vstrain(6)
273  nrow = fetch_tablerow( mc_viscoelastic, gauss%pMaterial%dict )
274  do n = 1,nrow
275  do i = 1,6
276  gauss%fstatus((n-1)*12+i) = gauss%fstatus((n-1)*12+i+6)
277  end do
278  end do
279  vstrain = gauss%strain
280  thetan = (vstrain(1) + vstrain(2) + vstrain(3))/3.d0
281  do i = 1,3
282  gauss%fstatus(nrow*12+i) = vstrain(i) - thetan
283  end do
284  do i = 4,6
285  gauss%fstatus(nrow*12+i) = vstrain(i)*0.5d0
286  end do
287  end subroutine
288 
289 end module
mmaterial::mc_viscoelastic
character(len=dict_key_length) mc_viscoelastic
Definition: material.f90:145
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:140
mviscoelastic
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
mmaterial::m_youngs
integer(kind=kint), parameter m_youngs
Definition: material.f90:92
mviscoelastic::calviscoelasticmatrix
subroutine calviscoelasticmatrix(matl, sectType, dt, D, temp)
This subroutine calculates tangent moduli for isotropic viscoelastic material.
Definition: Viscoelastic.f90:92
mviscoelastic::updateviscoelastic
subroutine updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
This subroutine provides to update stress for viscoelastic material.
Definition: Viscoelastic.f90:173
mmaterial::m_poisson
integer(kind=kint), parameter m_poisson
Definition: material.f90:93
mviscoelastic::updateviscoelasticstate
subroutine updateviscoelasticstate(gauss)
Update viscoplastic state.
Definition: Viscoelastic.f90:268
mmaterial::viscoelastic
integer(kind=kint), parameter viscoelastic
Definition: material.f90:77
m_fstr::eps
real(kind=kreal) eps
Definition: m_fstr.f90:142
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
m_fstr::dt
real(kind=kreal) dt
ANALYSIS CONTROL for NLGEOM and HEAT.
Definition: m_fstr.f90:139
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_elasticlinear
This module provides functions for elastic material.
Definition: ElasticLinear.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
mmaterial::tmaterial
Structure to manage all material related data.
Definition: material.f90:167
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
Definition: ElasticLinear.f90:16