FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
creep.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 mcreep
7  use hecmw_util
8  use mmaterial
10 
11  implicit none
12 
13 contains
14 
17  subroutine iso_creep(matl, sectType, stress, strain, extval,plstrain, &
18  dtime,ttime,stiffness, temp)
19  type( tmaterial ), intent(in) :: matl
20  integer, intent(in) :: sectType
21  real(kind=kreal), intent(in) :: stress(6)
22  real(kind=kreal), intent(in) :: strain(6)
23  real(kind=kreal), intent(in) :: extval(:)
24  real(kind=kreal), intent(in) :: plstrain
25  real(kind=kreal), intent(in) :: ttime
26  real(kind=kreal), intent(in) :: dtime
27  real(kind=kreal), intent(out) :: stiffness(6,6)
28  real(kind=kreal), optional :: temp
29 
30  integer :: i, j
31  logical :: ierr
32  real(kind=kreal) :: ina(1), outa(3)
33  real(kind=kreal) :: xxn, aa
34 
35  real(kind=kreal) :: c3,e,un,g,stri(6),p,dstri,c4,c5,f,df, eqvs
36 
37  aa = 0.0d0
38  xxn = 0.0d0
39  !
40  ! elastic
41  !
42  call calelasticmatrix( matl, secttype, stiffness )
43  if( dtime==0.d0 .or. all(stress==0.d0) ) return
44  !
45  ! elastic constants
46  !
47  if( present(temp) ) then
48  ina(1) = temp
49  call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr, ina )
50  else
51  call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr )
52  endif
53  if( ierr ) then
54  stop "error in isotropic elasticity definition"
55  else
56  e=outa(1)
57  un=outa(2)
58  endif
59 
60  ! Norton
61  if( matl%mtype==norton ) then ! those with no yield surface
62  if( present( temp ) ) then
63  ina(1) = temp
64  call fetch_tabledata( mc_norton, matl%dict, outa, ierr, ina )
65  else
66  call fetch_tabledata( mc_norton, matl%dict, outa, ierr )
67  endif
68  xxn=outa(2)
69  aa=outa(1)*((ttime+dtime)**(outa(3)+1.d0)-ttime**(outa(3)+1.d0))/(outa(3)+1.d0)
70  endif
71 
72  g=0.5d0*e/(1.d0+un)
73  !
74  ! creep
75  !
76  stri(:)=stress(:)
77  p=-(stri(1)+stri(2)+stri(3))/3.d0
78  do i=1,3
79  stri(i)=stri(i)+p
80  enddo
81  !
82  dstri=dsqrt(1.5d0*(stri(1)*stri(1)+stri(2)*stri(2)+stri(3)*stri(3)+ &
83  2.d0*(stri(4)*stri(4)+stri(5)*stri(5)+stri(6)*stri(6))) )
84  !
85  ! unit trial vector
86  !
87  stri(:)=stri(:)/dstri
88 
89  eqvs = dstri
90  if( eqvs<1.d-10 ) eqvs=1.d-10
91  f=aa*eqvs**xxn
92  df=xxn*f/eqvs
93  !
94  ! stiffness matrix
95  !
96  c3=6.d0*g*g
97  c4=c3*plstrain/(dstri+3.d0*g*plstrain)
98  c3=c4-c3*df/(3.d0*g*df+1.d0)
99  c5=c4/3.d0
100 
101  do i=1,6
102  do j=1,6
103  stiffness(i,j) = stiffness(i,j) +c3*stri(i)*stri(j)
104  enddo
105  enddo
106  do i=1,3
107  stiffness(i,i) = stiffness(i,i) - c4
108  do j=1,3
109  stiffness(i,j) = stiffness(i,j) +c5
110  enddo
111  enddo
112  do i=4,6
113  stiffness(i,i) = stiffness(i,i) - c4/2.d0
114  enddo
115 
116  end subroutine
117 
120  subroutine update_iso_creep(matl, sectType, strain, stress, extval,plstrain, &
121  dtime,ttime,temp)
122  type( tmaterial ), intent(in) :: matl
123  integer, intent(in) :: secttype
124  real(kind=kreal), intent(in) :: strain(6)
125  real(kind=kreal), intent(inout) :: stress(6)
126  real(kind=kreal), intent(inout) :: extval(:)
127  real(kind=kreal), intent(out) :: plstrain
128  real(kind=kreal), intent(in) :: ttime
129  real(kind=kreal), intent(in) :: dtime
130  real(kind=kreal), optional :: temp
131 
132  integer :: i
133  logical :: ierr
134  real(kind=kreal) :: ina(1), outa(3)
135  real(kind=kreal) :: xxn, aa
136 
137  real(kind=kreal) :: e,un,g,dg,ddg,stri(6),p,dstri,f,df, eqvs
138 
139  aa = 0.0d0
140  xxn = 0.0d0
141 
142  if( dtime==0.d0 ) return
143  !
144  ! elastic constants
145  !
146  if( present(temp) ) then
147  ina(1) = temp
148  call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr, ina )
149  else
150  call fetch_tabledata( mc_isoelastic, matl%dict, outa(1:2), ierr )
151  endif
152  if( ierr ) then
153  stop "error in isotropic elasticity definition"
154  else
155  e=outa(1)
156  un=outa(2)
157  endif
158 
159  ! Norton
160  if( matl%mtype==norton ) then ! those with no yield surface
161  if( present( temp ) ) then
162  ina(1) = temp
163  call fetch_tabledata( mc_norton, matl%dict, outa, ierr, ina )
164  else
165  call fetch_tabledata( mc_norton, matl%dict, outa, ierr )
166  endif
167  if( ierr ) then
168  stop "error in isotropic elasticity definition"
169  else
170  xxn=outa(2)
171  aa=outa(1)*((ttime+dtime)**(outa(3)+1.d0)-ttime**(outa(3)+1.d0))/(outa(3)+1.d0)
172  endif
173  endif
174 
175  g=0.5d0*e/(1.d0+un)
176 
177  !
178  ! creep
179  !
180  stri(:)=stress(:)
181  p=-(stri(1)+stri(2)+stri(3))/3.d0
182  do i=1,3
183  stri(i)=stri(i)+p
184  enddo
185  !
186  dstri=dsqrt(1.5d0*(stri(1)*stri(1)+stri(2)*stri(2)+stri(3)*stri(3)+ &
187  2.d0*(stri(4)*stri(4)+stri(5)*stri(5)+stri(6)*stri(6))) )
188  !
189  ! determination of the consistency parameter
190  !
191  dg=0.d0
192  do
193  if( matl%mtype==norton ) then
194  eqvs = dstri-3.d0*g*dg
195  f=aa*eqvs**xxn
196  df=xxn*f/eqvs
197  ddg = (f-dg)/(3.d0*g*df+1.d0)
198  dg = dg+ddg
199  if((ddg<dg*1.d-6).or.(ddg<1.d-12)) exit
200  endif
201  enddo
202 
203  stri(:) = stri(:)-3.d0*g*dg*stri(:)/dstri
204  stress(1:3) = stri(1:3)-p
205  stress(4:6) = stri(4:6)
206 
207  !
208  ! state variables
209  !
210  plstrain= dg
211  extval(1)=eqvs
212 
213  end subroutine
214 
216  subroutine updateviscostate( gauss )
217  use mmechgauss
218  type(tgaussstatus), intent(inout) :: gauss ! status of curr gauss point
219 
220  gauss%fstatus(2) = gauss%fstatus(2)+gauss%plstrain
221  end subroutine
222 
223 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 provides functions for creep calculation.
Definition: creep.f90:6
subroutine updateviscostate(gauss)
Update viscoplastic state.
Definition: creep.f90:217
subroutine iso_creep(matl, sectType, stress, strain, extval, plstrain, dtime, ttime, stiffness, temp)
This subroutine calculates stiffness for elastically isotropic materials with isotropic creep.
Definition: creep.f90:19
subroutine update_iso_creep(matl, sectType, strain, stress, extval, plstrain, dtime, ttime, temp)
This subroutine calculates stresses and creep status for an elastically isotropic material with isotr...
Definition: creep.f90:122
This module summarizes all infomation of material properties.
Definition: material.f90:6
character(len=dict_key_length) mc_norton
Definition: material.f90:125
integer(kind=kint), parameter norton
Definition: material.f90:71
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
Stucture to management all material relates data.
Definition: material.f90:144
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13