FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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!-------------------------------------------------------------------------------
6module mcreep
7 use hecmw_util
8 use mmaterial
10
11 implicit none
12
13contains
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
223end 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 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
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 updateviscostate(gauss)
Update viscoplastic state.
Definition: creep.f90:217
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