FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
calMatMatrix.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 use mmechgauss
14 use mcreep
15 use muelastic
16 use mumat
17
18 implicit none
19
20contains
21
23 integer function getnlgeomflag( gauss )
24 type( tgaussstatus ), intent(in) :: gauss
25 getnlgeomflag = gauss%pMaterial%nlgeom_flag
26 end function
27
29 subroutine matlmatrix( gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp )
30 type( tgaussstatus ), intent(in) :: gauss
31 integer, intent(in) :: sectType
32 real(kind=kreal), intent(out) :: matrix(:,:)
33 real(kind=kreal), intent(in) :: time
34 real(kind=kreal), intent(in) :: dtime
35 real(kind=kreal), intent(in) :: cdsys(3,3)
36 real(kind=kreal), intent(in), optional :: temperature
37 integer(kind=kint), intent(in), optional :: isEp
38
39 integer :: i
40 integer :: flag
41 real(kind=kreal) :: cijkl(3,3,3,3)
42 type( tmaterial ), pointer :: matl
43 matl=>gauss%pMaterial
44
45 flag = 0
46 if( present(isep) )then
47 if( isep == 1 )flag = 1
48 endif
49
50 if( matl%mtype==userelastic ) then
51 call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
52 elseif( isviscoelastic(matl%mtype) ) then
53 if( present(temperature) ) then
54 call calviscoelasticmatrix( matl, secttype, dtime, matrix, temperature )
55 else
56 call calviscoelasticmatrix( matl, secttype, dtime, matrix )
57 endif
58 elseif( iselastic(matl%mtype) .or. flag==1 ) then
59 if(flag==1)then
61 else
62 i = getelastictype(gauss%pMaterial%mtype)
63 endif
64
65 if( i==0 ) then
66 if( present(temperature) ) then
67 call calelasticmatrix( matl, secttype, matrix, temperature )
68 else
69 call calelasticmatrix( matl, secttype, matrix )
70 endif
71 elseif( i==1 ) then
72 if( present(temperature) ) then
73 call calelasticmatrix_ortho( gauss%pMaterial, secttype, cdsys, matrix, temperature )
74 else
75 call calelasticmatrix_ortho( gauss%pMaterial, secttype, cdsys, matrix )
76 endif
77 else
78 print *, "Elasticity type", matl%mtype, "not supported"
79 stop
80 endif
81
82 elseif( matl%mtype==neohooke .or. matl%mtype==mooneyrivlin ) then
83 call calelasticmooneyrivlin( matl, secttype, cijkl, gauss%strain )
84 call mat_c2d( cijkl, matrix, secttype )
85 elseif( matl%mtype==arrudaboyce ) then
86 call calelasticarrudaboyce( matl, secttype, cijkl, gauss%strain )
87 call mat_c2d( cijkl, matrix, secttype )
88 elseif( matl%mtype==mooneyrivlin_aniso ) then
89 call calelasticmooneyrivlinaniso( matl, secttype, cijkl, gauss%strain, cdsys )
90 call mat_c2d( cijkl, matrix, secttype )
91 elseif( matl%mtype==userhyperelastic ) then
92 call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
93 elseif( iselastoplastic(matl%mtype) ) then
94 if( present( temperature ) ) then
95 call calelastoplasticmatrix( matl, secttype, gauss%stress, &
96 gauss%istatus(1), gauss%fstatus, gauss%plstrain, matrix, temperature )
97 else
98 call calelastoplasticmatrix( matl, secttype, gauss%stress, &
99 gauss%istatus(1), gauss%fstatus, gauss%plstrain, matrix )
100 endif
101 elseif( matl%mtype==usermaterial ) then
102 call umatlmatrix( matl%name, matl%variables(101:), gauss%strain, &
103 gauss%stress, gauss%fstatus, matrix, dtime, time )
104 elseif( matl%mtype==norton ) then
105 if( present( temperature ) ) then
106 call iso_creep( matl, secttype, gauss%stress, gauss%strain, gauss%fstatus, &
107 gauss%plstrain, dtime, time, matrix, temperature )
108 else
109 call iso_creep( matl, secttype, gauss%stress, gauss%strain, gauss%fstatus, &
110 gauss%plstrain, dtime, time, matrix )
111 endif
112 else
113 stop "Material type not supported!"
114 endif
115
116 end subroutine
117
118 !
120 subroutine stressupdate( gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn )
121 type( tgaussstatus ), intent(inout) :: gauss
122 integer, intent(in) :: sectType
123 real(kind=kreal), intent(in) :: strain(6)
124 real(kind=kreal), intent(out) :: stress(6)
125 real(kind=kreal), intent(in) :: cdsys(3,3)
126 real(kind=kreal), intent(in), optional :: time
127 real(kind=kreal), intent(in), optional :: dtime
128 real(kind=kreal), optional :: temp
129 real(kind=kreal), optional :: tempn
130
131 if( gauss%pMaterial%mtype==neohooke .or. gauss%pMaterial%mtype==mooneyrivlin ) then
132 call calupdateelasticmooneyrivlin( gauss%pMaterial, secttype, strain, stress )
133 elseif( gauss%pMaterial%mtype==arrudaboyce ) then ! Arruda-Boyce Hyperelastic material
134 call calupdateelasticarrudaboyce( gauss%pMaterial, secttype, strain, stress )
135 elseif( gauss%pMaterial%mtype==mooneyrivlin_aniso ) then
136 call calupdateelasticmooneyrivlinaniso( gauss%pMaterial, secttype, strain, stress, cdsys )
137 elseif( gauss%pMaterial%mtype==userhyperelastic .or. gauss%pMaterial%mtype==userelastic ) then ! user-defined
138 call uelasticupdate( gauss%pMaterial%variables(101:), strain, stress )
139 elseif( isviscoelastic( gauss%pMaterial%mtype) ) then
140 if( .not. present(dtime) ) stop "error in viscoelastic update!"
141 if( present(temp) .and. present(tempn) ) then
142 call updateviscoelastic( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, dtime, temp, tempn )
143 else
144 call updateviscoelastic( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, dtime )
145 endif
146 elseif ( gauss%pMaterial%mtype==norton ) then
147 if( .not. present(dtime) ) stop "error in viscoelastic update!"
148 if( present(temp) ) then
149 call update_iso_creep( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, gauss%plstrain, dtime, time, temp )
150 else
151 call update_iso_creep( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, gauss%plstrain, dtime, time )
152 endif
153 elseif ( gauss%pMaterial%mtype==usermaterial) then ! user-defined
154 call uupdate( gauss%pMaterial%name, gauss%pMaterial%variables(101:), &
155 strain, stress, gauss%fstatus, dtime, time )
156 end if
157
158 end subroutine stressupdate
159
161 subroutine mat_c2d( cijkl, dij, itype )
162 real(kind=kreal), intent(in) :: cijkl(3,3,3,3)
163 real(kind=kreal), intent(out) :: dij(6,6)
164 integer, intent(in) :: itype
165
166 dij(:,:) = 0.d0
167 select case( itype )
168 case( d3 )
169 dij(1,1) = cijkl(1,1,1,1) ! ---
170 dij(1,2) = cijkl(1,1,2,2)
171 dij(1,3) = cijkl(1,1,3,3)
172 dij(1,4) = cijkl(1,1,1,2)
173 dij(1,5) = cijkl(1,1,2,3)
174 dij(1,6) = cijkl(1,1,3,1)
175 dij(2,1) = cijkl(2,2,1,1) ! ---
176 dij(2,2) = cijkl(2,2,2,2)
177 dij(2,3) = cijkl(2,2,3,3)
178 dij(2,4) = cijkl(2,2,1,2)
179 dij(2,5) = cijkl(2,2,2,3)
180 dij(2,6) = cijkl(2,2,3,1)
181 dij(3,1) = cijkl(3,3,1,1) ! ---
182 dij(3,2) = cijkl(3,3,2,2)
183 dij(3,3) = cijkl(3,3,3,3)
184 dij(3,4) = cijkl(3,3,1,2)
185 dij(3,5) = cijkl(3,3,2,3)
186 dij(3,6) = cijkl(3,3,3,1)
187 dij(4,1) = cijkl(1,2,1,1) ! ---
188 dij(4,2) = cijkl(1,2,2,2)
189 dij(4,3) = cijkl(1,2,3,3)
190 dij(4,4) = cijkl(1,2,1,2)
191 dij(4,5) = cijkl(1,2,2,3)
192 dij(4,6) = cijkl(1,2,3,1)
193 dij(5,1) = cijkl(2,3,1,1) ! ---
194 dij(5,2) = cijkl(2,3,2,2)
195 dij(5,3) = cijkl(2,3,3,3)
196 dij(5,4) = cijkl(2,3,1,2)
197 dij(5,5) = cijkl(2,3,2,3)
198 dij(5,6) = cijkl(2,3,3,1)
199 dij(6,1) = cijkl(3,1,1,1) ! ---
200 dij(6,2) = cijkl(3,1,2,2)
201 dij(6,3) = cijkl(3,1,3,3)
202 dij(6,4) = cijkl(3,1,1,2)
203 dij(6,5) = cijkl(3,1,2,3)
204 dij(6,6) = cijkl(3,1,3,1)
205 !
206 case( planestress, planestrain )
207 dij(1,1) = cijkl(1,1,1,1) ! ---
208 dij(1,2) = cijkl(1,1,2,2)
209 dij(1,3) = cijkl(1,1,1,2)
210 dij(2,1) = cijkl(2,2,1,1) ! ---
211 dij(2,2) = cijkl(2,2,2,2)
212 dij(2,3) = cijkl(2,2,1,2)
213 dij(3,1) = cijkl(1,2,1,1) ! ---
214 dij(3,2) = cijkl(1,2,2,2)
215 dij(3,3) = cijkl(1,2,1,2)
216 case( axissymetric )
217 dij(1,1) = cijkl(1,1,1,1)
218 dij(1,2) = cijkl(1,1,2,2)
219 dij(1,3) = cijkl(1,1,1,2)
220 dij(1,4) = cijkl(1,1,3,3)
221 dij(2,1) = cijkl(2,2,1,1)
222 dij(2,2) = cijkl(2,2,2,2)
223 dij(2,3) = cijkl(2,2,1,2)
224 dij(2,4) = cijkl(2,2,3,3)
225 dij(3,1) = cijkl(1,2,1,1)
226 dij(3,2) = cijkl(1,2,2,2)
227 dij(3,3) = cijkl(1,2,1,2)
228 dij(3,4) = cijkl(1,2,3,3)
229 dij(4,1) = cijkl(3,3,1,1)
230 dij(4,2) = cijkl(3,3,2,2)
231 dij(4,3) = cijkl(3,3,1,2)
232 dij(4,4) = cijkl(3,3,3,3)
233 case( shell )
234 end select
235
236 end subroutine mat_c2d
237
238 ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
239 !####################################################################
240 subroutine matlmatrix_shell &
241 (gauss, secttype, d, &
242 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
243 alpha, n_layer)
244 !####################################################################
245
246 type(tgaussstatus), intent(in) :: gauss
247 integer, intent(in) :: sectType, n_layer
248 real(kind = kreal), intent(out) :: d(:, :)
249 real(kind = kreal), intent(in) :: e1_hat(3), e2_hat(3), e3_hat(3)
250 real(kind = kreal), intent(in) :: cg1(3), cg2(3), cg3(3)
251 real(kind = kreal), intent(out) :: alpha
252
253 !--------------------------------------------------------------------
254
255 real(kind = kreal) :: c(3, 3, 3, 3)
256 type(tmaterial), pointer :: matl
257
258 !--------------------------------------------------------------------
259
260 matl => gauss%pMaterial
261
262 !--------------------------------------------------------------------
263
264 if( iselastic(matl%mtype) ) then
266 (matl, secttype, c, &
267 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
268 alpha, n_layer)
269
270 call mat_c2d_shell(c, d, secttype)
271 else
272 stop "Material type not supported!"
273 end if
274
275 !--------------------------------------------------------------------
276
277 return
278
279 !####################################################################
280 end subroutine matlmatrix_shell
281 !####################################################################
282 ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
283
284
285 ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
286 !####################################################################
287 subroutine mat_c2d_shell(c, D, itype)
288 !####################################################################
289
290 real(kind = kreal), intent(in) :: c(:, :, :, :)
291 real(kind = kreal), intent(out) :: d(:, :)
292 integer, intent(in) :: itype
293
294 !--------------------------------------------------------------------
295
296 integer :: index_i(5), index_j(5), &
297 index_k(5), index_l(5)
298 integer :: i, j, k, l
299 integer :: is, js
300
301 !--------------------------------------------------------------------
302
303 index_i(1) = 1
304 index_i(2) = 2
305 index_i(3) = 1
306 index_i(4) = 2
307 index_i(5) = 3
308
309 index_j(1) = 1
310 index_j(2) = 2
311 index_j(3) = 2
312 index_j(4) = 3
313 index_j(5) = 1
314
315 index_k(1) = 1
316 index_k(2) = 2
317 index_k(3) = 1
318 index_k(4) = 2
319 index_k(5) = 3
320
321 index_l(1) = 1
322 index_l(2) = 2
323 index_l(3) = 2
324 index_l(4) = 3
325 index_l(5) = 1
326
327 !--------------------------------------------------------------------
328
329 d(:, :) = 0.0d0
330
331 !--------------------------------------------------------------------
332
333 select case( itype )
334 case( shell )
335
336 do js = 1, 5
337
338 do is = 1, 5
339
340 i = index_i(is)
341 j = index_j(is)
342 k = index_k(js)
343 l = index_l(js)
344
345 d(is, js) = c(i, j, k, l)
346
347 end do
348
349 end do
350
351 end select
352
353 !--------------------------------------------------------------------
354
355 return
356
357 !####################################################################
358 end subroutine mat_c2d_shell
359 !####################################################################
360 ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
361
362end module m_matmatrix
I/O and Utility.
Definition: hecmw_util_f.F90:7
This module provides functions for elastic material.
subroutine calelasticmatrix_ortho(matl, secttype, bij, dmat, temp)
Calculate orthotropic elastic matrix.
subroutine linearelastic_shell(matl, secttype, c, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
subroutine calelasticmatrix(matl, secttype, d, temp)
Calculate isotropic elastic matrix.
This module provide functions for elastoplastic calculation.
subroutine calelastoplasticmatrix(matl, secttype, stress, istat, extval, plstrain, d, temperature)
This subroutine calculates elastoplastic constitutive relation.
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine stressupdate(gauss, secttype, strain, stress, cdsys, time, dtime, temp, tempn)
Update strain and stress for elastic and hyperelastic materials.
subroutine matlmatrix_shell(gauss, secttype, d, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
subroutine matlmatrix(gauss, secttype, matrix, time, dtime, cdsys, temperature, isep)
Calculate constituive matrix.
integer function getnlgeomflag(gauss)
Fetch the nlgeom flag of the material.
subroutine mat_c2d_shell(c, d, itype)
subroutine mat_c2d(cijkl, dij, itype)
Transfer rank 4 constituive matrix to rank 2 form.
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
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
subroutine calupdateelasticarrudaboyce(matl, secttype, dstrain, dstress)
This subroutine provides to update stress and strain for Arrude-Royce material.
subroutine calelasticmooneyrivlin(matl, secttype, cijkl, strain)
This subroutine provides elastic tangent coefficient for Mooney-Rivlin hyperelastic material.
subroutine calelasticarrudaboyce(matl, secttype, cijkl, strain)
This subroutine provides elastic tangent coefficient for Arruda-Boyce hyperelastic material.
subroutine calelasticmooneyrivlinaniso(matl, secttype, cijkl, strain, cdsys)
This subroutine provides elastic tangent coefficient for anisotropic Mooney-Rivlin hyperelastic mater...
subroutine calupdateelasticmooneyrivlinaniso(matl, secttype, strain, stress, cdsys)
This subroutine provides to update stress and strain for anisotropic Mooney-Rivlin material.
subroutine calupdateelasticmooneyrivlin(matl, secttype, strain, stress)
This subroutine provides to update stress and strain for Mooney-Rivlin material.
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter mooneyrivlin
Definition: material.f90:65
integer(kind=kint), parameter mooneyrivlin_aniso
Definition: material.f90:68
integer(kind=kint), parameter planestress
Definition: material.f90:77
integer function getelastictype(mtype)
Get elastic type.
Definition: material.f90:282
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter planestrain
Definition: material.f90:78
integer(kind=kint), parameter elastic
Definition: material.f90:58
integer(kind=kint), parameter arrudaboyce
Definition: material.f90:66
integer(kind=kint), parameter shell
Definition: material.f90:80
integer(kind=kint), parameter norton
Definition: material.f90:71
integer(kind=kint), parameter axissymetric
Definition: material.f90:79
integer(kind=kint), parameter userelastic
Definition: material.f90:60
integer(kind=kint), parameter neohooke
Definition: material.f90:64
integer(kind=kint), parameter usermaterial
Definition: material.f90:56
logical function iselastic(mtype)
If it is an elastic material?
Definition: material.f90:327
integer(kind=kint), parameter userhyperelastic
Definition: material.f90:67
logical function isviscoelastic(mtype)
If it is an viscoelastic material?
Definition: material.f90:354
logical function iselastoplastic(mtype)
If it is an elastoplastic material?
Definition: material.f90:336
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
This module provides interface for elastic or hyperelastic calculation.
Definition: uelastic.f90:7
subroutine uelasticmatrix(matl, strain, d)
This subroutine calculates constitutive relation.
Definition: uelastic.f90:13
subroutine uelasticupdate(matl, strain, stress)
This subroutine calculate updated strain and stress.
Definition: uelastic.f90:41
This subroutine read in used-defined material properties tangent.
Definition: umat.f90:7
subroutine uupdate(mname, matl, strain, stress, fstat, dtime, ttime, temperature)
This subroutine calculate strain and stress increment.
Definition: umat.f90:31
subroutine umatlmatrix(mname, matl, strain, stress, fstat, d, dtime, ttime, temperature)
This subroutine calculates constitutive matrix.
Definition: umat.f90:17
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
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