FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_C3D8.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!-------------------------------------------------------------------------------
8
9 use hecmw, only : kint, kreal
10 use elementinfo
11
12 implicit none
13
14contains
15
16
22 !----------------------------------------------------------------------*
23 subroutine stf_c3d8bbar &
24 (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
25 time, tincr, u, temperature)
26 !----------------------------------------------------------------------*
27
28 use mmechgauss
29 use m_matmatrix
31 use m_static_lib_3d, only: geomat_c3
32
33 !---------------------------------------------------------------------
34
35 integer(kind=kint), intent(in) :: etype
36 integer(kind=kint), intent(in) :: nn
37 real(kind=kreal), intent(in) :: ecoord(3, nn)
38 type(tgaussstatus), intent(in) :: gausses(:)
39 real(kind=kreal), intent(out) :: stiff(:,:)
40 integer(kind=kint), intent(in) :: cdsys_id
41 real(kind=kreal), intent(inout) :: coords(3, 3)
42 real(kind=kreal), intent(in) :: time
43 real(kind=kreal), intent(in) :: tincr
44 real(kind=kreal), intent(in), optional :: u(:, :)
45 real(kind=kreal), intent(in), optional :: temperature(nn)
46
47 !---------------------------------------------------------------------
48
49 integer(kind=kint) :: flag
50 integer(kind=kint), parameter :: ndof = 3
51 real(kind=kreal) :: d(6, 6),b(6, ndof*nn),db(6, ndof*nn)
52 real(kind=kreal) :: gderiv(nn, 3),stress(6),mat(6, 6)
53 real(kind=kreal) :: det, wg, temp, spfunc(nn)
54 integer(kind=kint) :: i, j, lx, serr
55 real(kind=kreal) :: naturalcoord(3)
56 real(kind=kreal) :: gdispderiv(3, 3)
57 real(kind=kreal) :: b1(6, ndof*nn), bbar(nn, 3)
58 real(kind=kreal) :: smat(9, 9), elem(3, nn)
59 real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
60 real(kind=kreal) :: b4, b6, b8, coordsys(3, 3)
61
62 !---------------------------------------------------------------------
63
64 stiff(:, :) = 0.0d0
65 ! we suppose the same material type in the element
66 flag = gausses(1)%pMaterial%nlgeom_flag
67 if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
68 elem(:, :) = ecoord(:, :)
69 if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
70
71 ! dilatation component at centroid
72 naturalcoord = 0.0d0
73 call getglobalderiv(etype, nn, naturalcoord, elem, det, bbar)
74
75 do lx = 1, numofquadpoints(etype)
76
77 call getquadpoint( etype, lx, naturalcoord(:) )
78 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
79
80 if( cdsys_id > 0 ) then
81 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
82 if( serr == -1 ) stop "Fail to setup local coordinate"
83 if( serr == -2 ) then
84 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
85 end if
86 end if
87
88 if( present(temperature) ) then
89 call getshapefunc( etype, naturalcoord, spfunc )
90 temp = dot_product( temperature, spfunc )
91 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
92 else
93 call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
94 end if
95
96 if( flag == updatelag ) then
97 call geomat_c3( gausses(lx)%stress, mat )
98 d(:, :) = d(:, :)-mat
99 endif
100
101 wg = getweight(etype, lx)*det
102 b(1:6, 1:nn*ndof) = 0.0d0
103 do j = 1, nn
104 b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
105 b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
106 b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
107 b(1, 3*j-2) = gderiv(j, 1)+b4
108 b(1, 3*j-1) = b6
109 b(1, 3*j ) = b8
110
111 b(2, 3*j-2) = b4
112 b(2, 3*j-1) = gderiv(j, 2)+b6
113 b(2, 3*j ) = b8
114
115 b(3, 3*j-2) = b4
116 b(3, 3*j-1) = b6
117 b(3, 3*j ) = gderiv(j, 3)+b8
118
119 b(4, 3*j-2) = gderiv(j, 2)
120 b(4, 3*j-1) = gderiv(j, 1)
121 b(5, 3*j-1) = gderiv(j, 3)
122 b(5, 3*j ) = gderiv(j, 2)
123 b(6, 3*j-2) = gderiv(j, 3)
124 b(6, 3*j ) = gderiv(j, 1)
125 end do
126
127 ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
128 if( flag == totallag ) then
129 ! ---dudx(i,j) ==> gdispderiv(i,j)
130 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
131 b1(1:6, 1:nn*ndof) = 0.0d0
132 do j = 1, nn
133 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
134 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
135 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
136 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
137 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
138 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
139 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
140 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
141 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
142 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
143 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
144 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
145 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
146 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
147 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
148 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
149 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
150 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
151 end do
152 ! ---BL = BL0 + BL1
153 do j = 1, nn*ndof
154 b(:, j) = b(:, j)+b1(:, j)
155 end do
156 end if
157
158 db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
159 forall( i=1:nn*ndof, j=1:nn*ndof )
160 stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
161 end forall
162
163 ! calculate the initial stress matrix
164 if( flag == totallag .or. flag == updatelag ) then
165 stress(1:6) = gausses(lx)%stress
166 bn(1:9, 1:nn*ndof) = 0.0d0
167 do j = 1, nn
168 bn(1, 3*j-2) = gderiv(j, 1)
169 bn(2, 3*j-1) = gderiv(j, 1)
170 bn(3, 3*j ) = gderiv(j, 1)
171 bn(4, 3*j-2) = gderiv(j, 2)
172 bn(5, 3*j-1) = gderiv(j, 2)
173 bn(6, 3*j ) = gderiv(j, 2)
174 bn(7, 3*j-2) = gderiv(j, 3)
175 bn(8, 3*j-1) = gderiv(j, 3)
176 bn(9, 3*j ) = gderiv(j, 3)
177 end do
178 smat(:, :) = 0.0d0
179 do j= 1, 3
180 smat(j , j ) = stress(1)
181 smat(j , j+3) = stress(4)
182 smat(j , j+6) = stress(6)
183 smat(j+3, j ) = stress(4)
184 smat(j+3, j+3) = stress(2)
185 smat(j+3, j+6) = stress(5)
186 smat(j+6, j ) = stress(6)
187 smat(j+6, j+3) = stress(5)
188 smat(j+6, j+6) = stress(3)
189 end do
190 sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
191 forall( i=1:nn*ndof, j=1:nn*ndof )
192 stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
193 end forall
194 end if
195
196 end do ! gauss roop
197
198 end subroutine stf_c3d8bbar
199
200
202 !----------------------------------------------------------------------*
203 subroutine update_c3d8bbar &
204 (etype, nn, ecoord, u, du, cdsys_id, coords, &
205 qf ,gausses, iter, time, tincr, tt,t0, tn )
206 !----------------------------------------------------------------------*
207
208 use m_fstr
209 use mmaterial
210 use mmechgauss
211 use m_matmatrix
213 use mhyperelastic
214 use m_utilities
216
217 !---------------------------------------------------------------------
218
219 integer(kind=kint), intent(in) :: etype
220 integer(kind=kint), intent(in) :: nn
221 real(kind=kreal), intent(in) :: ecoord(3, nn)
222 real(kind=kreal), intent(in) :: u(3, nn)
223 real(kind=kreal), intent(in) :: du(3, nn)
224 integer(kind=kint), intent(in) :: cdsys_id
225 real(kind=kreal), intent(inout) :: coords(3, 3)
226 real(kind=kreal), intent(out) :: qf(nn*3)
227 type(tgaussstatus), intent(inout) :: gausses(:)
228 integer, intent(in) :: iter
229 real(kind=kreal), intent(in) :: time
230 real(kind=kreal), intent(in) :: tincr
231 real(kind=kreal), intent(in), optional :: tt(nn)
232 real(kind=kreal), intent(in), optional :: t0(nn)
233 real(kind=kreal), intent(in), optional :: tn(nn)
234
235 !---------------------------------------------------------------------
236
237 integer(kind=kint) :: flag
238 integer(kind=kint), parameter :: ndof = 3
239 real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
240 real(kind=kreal) :: gderiv(nn, 3), gderiv1(nn,3), gdispderiv(3, 3), f(3,3), det, det1, wg
241 integer(kind=kint) :: j, lx, mtype, serr
242 real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
243 real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
244 real(kind=kreal) :: dstrain(6)
245 real(kind=kreal) :: dvol, vol0, bbar(nn, 3), derivdum(1:ndof, 1:ndof), bbar2(nn, 3)
246 real(kind=kreal) :: b4, b6, b8, ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
247 logical :: ierr, matlaniso
248
249 !---------------------------------------------------------------------
250
251 qf(:) = 0.0d0
252 ! we suppose the same material type in the element
253 flag = gausses(1)%pMaterial%nlgeom_flag
254 elem(:, :) = ecoord(:, :)
255 totaldisp(:, :) = u(:, :)+du(:, :)
256 if( flag == updatelag ) then
257 elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
258 elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
259 ! elem = elem1
260 totaldisp(:, :) = du(:, :)
261 end if
262
263 matlaniso = .false.
264 if( cdsys_id > 0 .AND. present(tt) ) then
265 ina = tt(1)
266 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
267 if( .not. ierr ) matlaniso = .true.
268 end if
269
270 ! dilatation at centroid
271 naturalcoord = 0.0d0
272 call getglobalderiv(etype, nn, naturalcoord, elem, det, bbar)
273 derivdum = matmul( totaldisp(1:ndof, 1:nn), bbar(1:nn, 1:ndof) )
274 vol0 = ( derivdum(1, 1)+derivdum(2, 2)+derivdum(3, 3) )/3.0d0
275 if( flag == updatelag ) call getglobalderiv(etype, nn, naturalcoord, elem1, det, bbar2)
276
277 do lx = 1, numofquadpoints(etype)
278
279 mtype = gausses(lx)%pMaterial%mtype
280
281 call getquadpoint( etype, lx, naturalcoord(:) )
282 call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
283
284 if( cdsys_id > 0 ) then
285 call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
286 if( serr == -1 ) stop "Fail to setup local coordinate"
287 if( serr == -2 ) then
288 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
289 end if
290 end if
291
292 gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
293 dvol = vol0-( gdispderiv(1, 1)+gdispderiv(2, 2)+gdispderiv(3, 3) )/3.0d0
294 !gdispderiv(1, 1) = gdispderiv(1, 1)+dvol
295 !gdispderiv(2, 2) = gdispderiv(2, 2)+dvol
296 !gdispderiv(3, 3) = gdispderiv(3, 3)+dvol
297
298 ! ========================================================
299 ! UPDATE STRAIN and STRESS
300 ! ========================================================
301
302 ! Thermal Strain
303 epsth = 0.0d0
304 if( present(tt) .AND. present(t0) ) then
305 call getshapefunc(etype, naturalcoord, spfunc)
306 ttc = dot_product(tt, spfunc)
307 tt0 = dot_product(t0, spfunc)
308 ttn = dot_product(tn, spfunc)
309 call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
310 end if
311
312 ! Update strain
313 ! Small strain
314 dstrain(1) = gdispderiv(1, 1)+dvol
315 dstrain(2) = gdispderiv(2, 2)+dvol
316 dstrain(3) = gdispderiv(3, 3)+dvol
317 dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
318 dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
319 dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
320 dstrain(:) = dstrain(:)-epsth(:)
321
322 f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
323 if( flag == infinitesimal ) then
324 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
325
326 else if( flag == totallag ) then
327 ! Green-Lagrange strain
328 dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
329 dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
330 dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
331 dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
332 +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
333 dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
334 +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
335 dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
336 +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
337
338 gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
339 f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
340
341 else if( flag == updatelag ) then
342 rot = 0.0d0
343 rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
344 rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
345 rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
346
347 gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
348 call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
349 f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn)+du(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
350
351 end if
352
353 ! Update stress
354 if( present(tt) .AND. present(t0) ) then
355 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
356 else
357 call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
358 end if
359
360 ! ========================================================
361 ! calculate the internal force ( equivalent nodal force )
362 ! ========================================================
363 ! Small strain
364 b(1:6, 1:nn*ndof) = 0.0d0
365 do j = 1, nn
366 b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
367 b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
368 b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
369 b(1, 3*j-2) = gderiv(j, 1)+b4
370 b(1, 3*j-1) = b6
371 b(1, 3*j ) = b8
372
373 b(2, 3*j-2) = b4
374 b(2, 3*j-1) = gderiv(j, 2)+b6
375 b(2, 3*j ) = b8
376
377 b(3, 3*j-2) = b4
378 b(3, 3*j-1) = b6
379 b(3, 3*j ) = gderiv(j, 3)+b8
380
381 b(4, 3*j-2) = gderiv(j, 2)
382 b(4, 3*j-1) = gderiv(j, 1)
383 b(5, 3*j-1) = gderiv(j, 3)
384 b(5, 3*j ) = gderiv(j, 2)
385 b(6, 3*j-2) = gderiv(j, 3)
386 b(6, 3*j ) = gderiv(j, 1)
387 end do
388
389 if( flag == infinitesimal ) then
390
391 else if( flag == totallag ) then
392
393 b1(1:6, 1:nn*ndof) = 0.0d0
394 do j = 1, nn
395 b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
396 b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
397 b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
398 b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
399 b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
400 b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
401 b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
402 b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
403 b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
404 b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
405 b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
406 b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
407 b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
408 b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
409 b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
410 b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
411 b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
412 b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
413 end do
414 ! BL = BL0 + BL1
415 do j = 1, nn*ndof
416 b(:, j) = b(:, j)+b1(:, j)
417 end do
418
419 else if( flag == updatelag ) then
420
421 call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
422 b(1:6, 1:nn*ndof) = 0.0d0
423 do j = 1, nn
424 b4 = ( bbar2(j, 1)-gderiv(j, 1) )/3.0d0
425 b6 = ( bbar2(j, 2)-gderiv(j, 2) )/3.0d0
426 b8 = ( bbar2(j, 3)-gderiv(j, 3) )/3.0d0
427 b(1, 3*j-2) = gderiv(j, 1)+b4
428 b(1, 3*j-1) = b6
429 b(1, 3*j ) = b8
430
431 b(2, 3*j-2) = b4
432 b(2, 3*j-1) = gderiv(j, 2)+b6
433 b(2, 3*j ) = b8
434
435 b(3, 3*j-2) = b4
436 b(3, 3*j-1) = b6
437 b(3, 3*j ) = gderiv(j, 3)+b8
438
439 b(4, 3*j-2) = gderiv(j, 2)
440 b(4, 3*j-1) = gderiv(j, 1)
441 b(5, 3*j-1) = gderiv(j, 3)
442 b(5, 3*j ) = gderiv(j, 2)
443 b(6, 3*j-2) = gderiv(j, 3)
444 b(6, 3*j ) = gderiv(j, 1)
445 end do
446 end if
447
448 !! calculate the Internal Force
449 wg = getweight(etype, lx)*det
450 qf(1:nn*ndof) = qf(1:nn*ndof) &
451 +matmul( gausses(lx)%stress(1:6), b(1:6, 1:nn*ndof) )*wg
452
453 end do
454
455 end subroutine update_c3d8bbar
456
458 !----------------------------------------------------------------------*
459 subroutine tload_c3d8bbar &
460 (etype, nn, xx, yy, zz, tt, t0, &
461 gausses, vect, cdsys_id, coords)
462 !----------------------------------------------------------------------*
463
464 use m_fstr
465 use mmechgauss
466 use m_matmatrix
467 use m_utilities
468
469 !---------------------------------------------------------------------
470
471 integer(kind=kint), parameter :: ndof = 3
472 integer(kind=kint), intent(in) :: etype, nn
473 type(tgaussstatus), intent(in) :: gausses(:)
474 real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
475 real(kind=kreal), intent(in) :: tt(nn), t0(nn)
476 real(kind=kreal), intent(out) :: vect(nn*ndof)
477 integer(kind=kint), intent(in) :: cdsys_ID
478 real(kind=kreal), intent(inout) :: coords(3, 3)
479
480 !---------------------------------------------------------------------
481
482 real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
483 real(kind=kreal) :: b4, b6, b8, det, ecoord(3, nn)
484 integer(kind=kint) :: j, LX, serr
485 real(kind=kreal) :: estrain(6), sgm(6), h(nn)
486 real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
487 real(kind=kreal) :: wg, outa(1), ina(1), bbar(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
488 real(kind=kreal) :: tempc, temp0, tempc0, temp00, thermal_eps, tm(6,6)
489 logical :: ierr, matlaniso
490
491 !---------------------------------------------------------------------
492
493 matlaniso = .false.
494
495 if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
496 ina = tt(1)
497 call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
498 if( .not. ierr ) matlaniso = .true.
499 end if
500
501 vect(:) = 0.0d0
502
503 ecoord(1, :) = xx(:)
504 ecoord(2, :) = yy(:)
505 ecoord(3, :) = zz(:)
506
507 naturalcoord = 0.0d0
508 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, bbar)
509 call getshapefunc( etype, naturalcoord, h(1:nn) )
510 tempc0 = dot_product( h(1:nn), tt(1:nn) )
511 temp00 = dot_product( h(1:nn), t0(1:nn) )
512
513 ! LOOP FOR INTEGRATION POINTS
514 do lx = 1, numofquadpoints(etype)
515
516 call getquadpoint( etype, lx, naturalcoord(:) )
517 call getshapefunc( etype, naturalcoord, h(1:nn) )
518 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
519
520 if( matlaniso ) then
521 call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
522 if( serr == -1 ) stop "Fail to setup local coordinate"
523 if( serr == -2 ) then
524 write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
525 end if
526 end if
527
528 ! WEIGHT VALUE AT GAUSSIAN POINT
529 wg = getweight(etype, lx)*det
530 b(1:6, 1:nn*ndof) = 0.0d0
531 do j = 1, nn
532 b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
533 b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
534 b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
535 b(1, 3*j-2) = gderiv(j, 1)+b4
536 b(1, 3*j-1) = b6
537 b(1, 3*j ) = b8
538
539 b(2, 3*j-2) = b4
540 b(2, 3*j-1) = gderiv(j, 2)+b6
541 b(2, 3*j ) = b8
542
543 b(3, 3*j-2) = b4
544 b(3, 3*j-1) = b6
545 b(3, 3*j ) = gderiv(j, 3)+b8
546
547 b(4, 3*j-2) = gderiv(j, 2)
548 b(4, 3*j-1) = gderiv(j, 1)
549 b(5, 3*j-1) = gderiv(j, 3)
550 b(5, 3*j ) = gderiv(j, 2)
551 b(6, 3*j-2) = gderiv(j, 3)
552 b(6, 3*j ) = gderiv(j, 1)
553 end do
554
555 tempc = dot_product( h(1:nn),tt(1:nn) )
556 temp0 = dot_product( h(1:nn),t0(1:nn) )
557
558 call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
559
560 ina(1) = tempc
561 if( matlaniso ) then
562 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
563 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
564 else
565 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
566 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
567 alp = outa(1)
568 end if
569 ina(1) = temp0
570 if( matlaniso ) then
571 call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
572 if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
573 else
574 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
575 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
576 alp0 = outa(1)
577 end if
578
579 !**
580 !** THERMAL strain
581 !**
582 if( matlaniso ) then
583 do j = 1, 3
584 estrain(j) = alpo(j)*(tempc0-ref_temp)-alpo0(j)*(temp00-ref_temp)
585 end do
586 estrain(4:6) = 0.0d0
587 call transformation(coordsys, tm)
588 estrain(:) = matmul( estrain(:), tm ) ! to global coord
589 estrain(4:6) = estrain(4:6)*2.0d0
590 else
591 thermal_eps = alp*(tempc0-ref_temp)-alp0*(temp00-ref_temp)
592 estrain(1:3) = thermal_eps
593 estrain(4:6) = 0.0d0
594 end if
595
596 !**
597 !** SET SGM {s}=[D]{e}
598 !**
599 sgm(:) = matmul( d(:, :), estrain(:) )
600
601 !**
602 !** CALCULATE LOAD {F}=[B]T{e}
603 !**
604 vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
605
606 end do
607
608 end subroutine tload_c3d8bbar
609
610
611end module m_static_lib_c3d8
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, secttype, matrix, time, dtime, cdsys, temperature, isep)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine update_stress3d(flag, gauss, rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn)
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, epsth)
This module contains several strategy to free locking problem in Eight-node hexagonal element.
subroutine stf_c3d8bbar(etype, nn, ecoord, gausses, stiff, cdsys_id, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using b-bar method.
subroutine update_c3d8bbar(etype, nn, ecoord, u, du, cdsys_id, coords, qf, gausses, iter, time, tincr, tt, t0, tn)
Update Strain stress of this element.
subroutine tload_c3d8bbar(etype, nn, xx, yy, zz, tt, t0, gausses, vect, cdsys_id, coords)
This subroutien calculate thermal loading.
This module provides aux functions.
Definition: utilities.f90:6
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13