FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
heat_LIB_CONDUCTIVITY.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!-------------------------------------------------------------------------------
6contains
7
8 subroutine heat_get_coefficient(Tpoi, imat, coef, ntab, temp, funcA, funcB)
9 use hecmw
10 implicit none
11 real(kind=kreal) :: coef(3)
12 integer(kind=kint) :: i, in, imat, itab, ntab
13 real(kind=kreal) :: tpoi, temp(ntab), funca(ntab+1), funcb(ntab+1)
14
15 itab = 0
16 if(tpoi < temp(1)) then
17 itab = 1
18 elseif(temp(ntab) <= tpoi)then
19 itab = ntab + 1
20 else
21 do in = 1, ntab - 1
22 if(temp(in) <= tpoi .and. tpoi < temp(in+1))then
23 itab = in + 1
24 exit
25 endif
26 enddo
27 endif
28
29 coef(1) = funca(itab)*tpoi + funcb(itab)
30 coef(2) = coef(1)
31 coef(3) = coef(1)
32 end subroutine heat_get_coefficient
33
34 subroutine heat_conductivity_c1(etype, nn, ecoord, temperature, IMAT, surf, stiff, &
35 ntab, temp, funcA, funcB)
36 use hecmw
37 implicit none
38 integer(kind=kint), intent(in) :: etype
39 integer(kind=kint), intent(in) :: nn
40 real(kind=kreal), intent(in) :: temperature(nn)
41 real(kind=kreal), intent(out) :: stiff(:,:)
42 integer(kind=kint), parameter :: ndof = 1
43 integer(kind=kint) :: i, j, IMAT
44 integer(kind=kint) :: ntab
45 real(kind=kreal) :: ecoord(3, nn)
46 real(kind=kreal) :: dx, dy, dz, surf, val, temp_i, length
47 real(kind=kreal) :: cc(3)
48 real(kind=kreal) :: temp(ntab), funca(ntab+1), funcb(ntab+1)
49
50 if(etype == 611)then
51 ecoord(1,2) = ecoord(1,3); ecoord(2,2) = ecoord(2,3); ecoord(3,2) = ecoord(3,3)
52 endif
53
54 dx = ecoord(1,2) - ecoord(1,1)
55 dy = ecoord(2,2) - ecoord(2,1)
56 dz = ecoord(3,2) - ecoord(3,1)
57 length = dsqrt(dx*dx + dy*dy + dz*dz)
58 val = surf*length
59
60 temp_i = (temperature(1) + temperature(2))*0.5d0
61 call heat_get_coefficient(temp_i, imat, cc, ntab, temp, funca, funcb)
62
63 stiff(1,1) = val*cc(1)
64 stiff(2,1) = -val*cc(1)
65 stiff(1,2) = -val*cc(1)
66 stiff(2,2) = val*cc(1)
67 end subroutine heat_conductivity_c1
68
69 subroutine heat_conductivity_c2(etype, nn, ecoord, temperature, IMAT, thick, stiff, &
70 ntab, temp, funcA, funcB)
71 use mmechgauss
72 use m_matmatrix
73 use elementinfo
74 implicit none
75 !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
76 integer(kind=kint), intent(in) :: etype
77 integer(kind=kint), intent(in) :: nn
78 real(kind=kreal), intent(in) :: ecoord(2,nn)
79 real(kind=kreal), intent(out) :: stiff(:,:)
80 real(kind=kreal), intent(in) :: temperature(nn)
81 type(tmaterial), pointer :: matl
82 integer(kind=kint) :: i, j, lx, imat, ntab
83 real(kind=kreal) :: naturalcoord(2)
84 real(kind=kreal) :: func(nn), thick, temp_i
85 real(kind=kreal) :: det, wg, rho, diag_stiff, total_stiff
86 real(kind=kreal) :: d(2,2), n(2,nn), dn(2,nn)
87 real(kind=kreal) :: gderiv(nn,2)
88 real(kind=kreal) :: cc(3)
89 real(kind=kreal) :: temp(ntab), funca(ntab+1), funcb(ntab+1)
90
91 stiff = 0.0d0
92 !matl => gausses(1)%pMaterial
93
94 do lx = 1, numofquadpoints(etype)
95 call getquadpoint(etype, lx, naturalcoord)
96 call getshapefunc(etype, naturalcoord, func)
97 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
98
99 temp_i = dot_product(func, temperature)
100 call heat_get_coefficient(temp_i, imat, cc, ntab, temp, funca, funcb)
101
102 d = 0.0d0
103 d(1,1) = cc(1)*thick
104 d(2,2) = cc(1)*thick
105 wg = getweight(etype, lx)*det
106
107 dn = matmul(d, transpose(gderiv))
108 forall(i = 1:nn, j = 1:nn)
109 stiff(i,j) = stiff(i,j) + dot_product(gderiv(i,:), dn(:,j))*wg
110 end forall
111 enddo
112 end subroutine heat_conductivity_c2
113
114 subroutine heat_conductivity_c3(etype, nn, ecoord, temperature, IMAT, stiff, &
115 ntab, temp, funcA, funcB)
116 use mmechgauss
117 use m_matmatrix
118 use elementinfo
119 implicit none
120 !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
121 integer(kind=kint), intent(in) :: etype
122 integer(kind=kint), intent(in) :: nn
123 real(kind=kreal), intent(in) :: ecoord(3,nn)
124 real(kind=kreal), intent(out) :: stiff(:,:)
125 real(kind=kreal), intent(in) :: temperature(nn)
126 type(tmaterial), pointer :: matl
127 integer(kind=kint) :: i, j, lx, imat, ntab
128 real(kind=kreal) :: naturalcoord(3)
129 real(kind=kreal) :: func(nn), temp_i
130 real(kind=kreal) :: det, wg, rho, diag_stiff, total_stiff
131 real(kind=kreal) :: d(3, 3), n(3, nn), dn(3, nn)
132 real(kind=kreal) :: gderiv(nn, 3)
133 real(kind=kreal) :: spe(3)
134 real(kind=kreal) :: temp(ntab), funca(ntab+1), funcb(ntab+1)
135
136 stiff = 0.0d0
137 !matl => gausses(1)%pMaterial
138
139 do lx = 1, numofquadpoints(etype)
140 call getquadpoint(etype, lx, naturalcoord)
141 call getshapefunc(etype, naturalcoord, func)
142 call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
143
144 temp_i = dot_product(func, temperature)
145 call heat_get_coefficient(temp_i, imat, spe, ntab, temp, funca, funcb)
146
147 d = 0.0d0
148 d(1,1) = spe(1)
149 d(2,2) = spe(1)
150 d(3,3) = spe(1)
151 wg = getweight(etype, lx)*det
152
153 dn = matmul(d, transpose(gderiv))
154 forall(i = 1:nn, j = 1:nn)
155 stiff(i,j) = stiff(i,j) + dot_product(gderiv(i,:), dn(:,j))*wg
156 end forall
157 enddo
158 end subroutine heat_conductivity_c3
159
161 subroutine heat_conductivity_shell_731(etype, nn, ecoord, TT, IMAT, thick, SS, stiff, &
162 ntab, temp, funcA, funcB)
163 use mmechgauss
164 use m_matmatrix
165 use elementinfo
167 implicit none
168 !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
169 integer(kind=kint), intent(in) :: etype
170 integer(kind=kint), intent(in) :: nn
171 real(kind=kreal), intent(in) :: ecoord(3,nn)
172 real(kind=kreal), intent(out) :: stiff(:,:)
173 real(kind=kreal), intent(inout) :: ss(:)
174 real(kind=kreal), intent(inout) :: tt(nn)
175 type(tmaterial), pointer :: matl
176 integer(kind=kint) :: i, j, lx, imat, ntab
177 integer(kind=kint) :: ig1, ig2, ig3, inod
178 real(kind=kreal) :: surf, thick, temp_i
179 real(kind=kreal) :: ri, rm, rp, si, sm, sp, ti, valx, valy, valz, var
180 real(kind=kreal) :: xj11, xj12, xj13, xj21, xj22, xj23, xj31, xj32, xj33, xsum
181 real(kind=kreal) :: det, wg, rho, diag_stiff, total_stiff
182 real(kind=kreal) :: cc(3), ctemp, dum
183 real(kind=kreal) :: temp(ntab), funca(ntab+1), funcb(ntab+1)
184 real(kind=kreal) :: xg(2), wgt(2), h(4), hr(4), hs(4)
185 real(kind=kreal) :: cod(3, 4)
186 real(kind=kreal) :: g1(3), g2(3), g3(3), e1(3), e2(3), e3(3), ref(3)
187 real(kind=kreal) :: en(3, 4), the(3, 3), amat(3, 3), bv(3), wk(3)
188 real(kind=kreal) :: dtdx(4), dtdy(4)
189 data xg/-0.5773502691896d0,0.5773502691896d0/
190 data wgt/1.0d0,1.0d0/
191
192 do i = 1, nn
193 cod(1,i) = ecoord(1,i)
194 cod(2,i) = ecoord(2,i)
195 cod(3,i) = ecoord(3,i)
196 enddo
197
198 tt(nn) = tt(nn-1)
199 cod(1,nn) = cod(1,nn-1)
200 cod(2,nn) = cod(2,nn-1)
201 cod(3,nn) = cod(3,nn-1)
202
203 !* SET REFFRENSE VECTOR TO DETERMINE LOCAL COORDINATE SYSTEM
204 ref(1) = 0.25*( cod(1,2) + cod(1,3) - cod(1,1) - cod(1,4) )
205 ref(2) = 0.25*( cod(2,2) + cod(2,3) - cod(2,1) - cod(2,4) )
206 ref(3) = 0.25*( cod(3,2) + cod(3,3) - cod(3,1) - cod(3,4) )
207
208 g1(1) = cod(1,1) - cod(1,2)
209 g1(2) = cod(2,1) - cod(2,2)
210 g1(3) = cod(3,1) - cod(3,2)
211 g2(1) = cod(1,2) - cod(1,3)
212 g2(2) = cod(2,2) - cod(2,3)
213 g2(3) = cod(3,2) - cod(3,3)
214
215 !* G3()=G1() X G2()
216 g3(1) = g1(2)*g2(3) - g1(3)*g2(2)
217 g3(2) = g1(3)*g2(1) - g1(1)*g2(3)
218 g3(3) = g1(1)*g2(2) - g1(2)*g2(1)
219
220 !* SET BASE VECTOR IN LOCAL CARTESIAN COORDINATE
221 xsum = dsqrt( g3(1)**2 + g3(2)**2 + g3(3)**2 )
222
223 do i = 1, nn
224 en(1,i) = g3(1) / xsum
225 en(2,i) = g3(2) / xsum
226 en(3,i) = g3(3) / xsum
227 enddo
228
229 !* LOOP FOR GAUSS INTEGRATION POINT
230 do ig3 = 1, 2
231 ti = xg(ig3)
232 do ig2 = 1, 2
233 si = xg(ig2)
234 do ig1 = 1, 2
235 ri = xg(ig1)
236
237 rp = 1.0 + ri
238 sp = 1.0 + si
239 rm = 1.0 - ri
240 sm = 1.0 - si
241
242 h(1) = 0.25*rm*sm
243 h(2) = 0.25*rp*sm
244 h(3) = 0.25*rp*sp
245 h(4) = 0.25*rm*sp
246
247 hr(1) = -0.25*sm
248 hr(2) = 0.25*sm
249 hr(3) = 0.25*sp
250 hr(4) = -0.25*sp
251
252 hs(1) = -0.25*rm
253 hs(2) = -0.25*rp
254 hs(3) = 0.25*rp
255 hs(4) = 0.25*rm
256
257 !* COVARIANT BASE VECTOR AT A GAUSS INTEGRARION POINT
258 do i = 1, 3
259 g1(i) = 0.0
260 g2(i) = 0.0
261 g3(i) = 0.0
262
263 do inod = 1, nn
264 var = cod(i,inod) + thick*0.5*ti*en(i,inod)
265 g1(i) = g1(i) + hr(inod)*var
266 g2(i) = g2(i) + hs(inod)*var
267 g3(i) = g3(i) + thick*0.5*h(inod)*en(i,inod)
268 enddo
269 enddo
270
271 !*JACOBI MATRIX
272 xj11 = g1(1)
273 xj12 = g1(2)
274 xj13 = g1(3)
275 xj21 = g2(1)
276 xj22 = g2(2)
277 xj23 = g2(3)
278 xj31 = g3(1)
279 xj32 = g3(2)
280 xj33 = g3(3)
281
282 !*DETERMINANT OF JACOBIAN
283 det = xj11*xj22*xj33 &
284 + xj12*xj23*xj31 &
285 + xj13*xj21*xj32 &
286 - xj13*xj22*xj31 &
287 - xj12*xj21*xj33 &
288 - xj11*xj23*xj32
289
290 !* INVERSION OF JACOBIAN
291 dum = 1.0 / det
292 amat(1,1) = dum*( xj22*xj33 - xj23*xj32 )
293 amat(2,1) = dum*( -xj21*xj33 + xj23*xj31 )
294 amat(3,1) = dum*( xj21*xj32 - xj22*xj31 )
295 amat(1,2) = dum*( -xj12*xj33 + xj13*xj32 )
296 amat(2,2) = dum*( xj11*xj33 - xj13*xj31 )
297 amat(3,2) = dum*( -xj11*xj32 + xj12*xj31 )
298 amat(1,3) = dum*( xj12*xj23 - xj13*xj22 )
299 amat(2,3) = dum*( -xj11*xj23 + xj13*xj21 )
300 amat(3,3) = dum*( xj11*xj22 - xj12*xj21 )
301
302 !* SET BASE VECTOR IN LOCAL CARTESIAN COORDINATE
303 xsum = dsqrt( g3(1)**2 + g3(2)**2 + g3(3)**2 )
304 e3(1) = g3(1) / xsum
305 e3(2) = g3(2) / xsum
306 e3(3) = g3(3) / xsum
307 e2(1) = -ref(2)*e3(3) + ref(3)*e3(2)
308 e2(2) = -ref(3)*e3(1) + ref(1)*e3(3)
309 e2(3) = -ref(1)*e3(2) + ref(2)*e3(1)
310 e1(1) = -e3(2)*e2(3) + e3(3)*e2(2)
311 e1(2) = -e3(3)*e2(1) + e3(1)*e2(3)
312 e1(3) = -e3(1)*e2(2) + e3(2)*e2(1)
313 xsum = dsqrt(e2(1)**2 + e2(2)**2 + e2(3)**2)
314
315 if ( xsum .GT. 1.e-15 ) then
316 e2(1) = e2(1) / xsum
317 e2(2) = e2(2) / xsum
318 e2(3) = e2(3) / xsum
319 e1(1) = -e3(2)*e2(3) + e3(3)*e2(2)
320 e1(2) = -e3(3)*e2(1) + e3(1)*e2(3)
321 e1(3) = -e3(1)*e2(2) + e3(2)*e2(1)
322 xsum = dsqrt( e1(1)**2 + e1(2)**2 + e1(3)**2 )
323 e1(1) = e1(1) / xsum
324 e1(2) = e1(2) / xsum
325 e1(3) = e1(3) / xsum
326 else
327 e1(1) = 0.d0
328 e1(2) = 0.d0
329 e1(3) = -1.d0
330 e2(1) = 0.d0
331 e2(2) = 1.d0
332 e2(3) = 0.d0
333 end if
334
335 the(1,1) = e1(1)
336 the(1,2) = e1(2)
337 the(1,3) = e1(3)
338 the(2,1) = e2(1)
339 the(2,2) = e2(2)
340 the(2,3) = e2(3)
341 the(3,1) = e3(1)
342 the(3,2) = e3(2)
343 the(3,3) = e3(3)
344
345 do i = 1, nn
346 bv(1) = hr(i)
347 bv(2) = hs(i)
348 bv(3) = 0.0
349 wk(1) = amat(1,1)*bv(1) &
350 + amat(1,2)*bv(2) &
351 + amat(1,3)*bv(3)
352 wk(2) = amat(2,1)*bv(1) &
353 + amat(2,2)*bv(2) &
354 + amat(2,3)*bv(3)
355 wk(3) = amat(3,1)*bv(1) &
356 + amat(3,2)*bv(2) &
357 + amat(3,3)*bv(3)
358 bv(1) = the(1,1)*wk(1) &
359 + the(1,2)*wk(2) &
360 + the(1,3)*wk(3)
361 bv(2) = the(2,1)*wk(1) &
362 + the(2,2)*wk(2) &
363 + the(2,3)*wk(3)
364 bv(3) = the(3,1)*wk(1) &
365 + the(3,2)*wk(2) &
366 + the(3,3)*wk(3)
367 dtdx(i) = bv(1)
368 dtdy(i) = bv(2)
369 enddo
370
371 !*CONDUCTIVITY AT CURRENT TEMPERATURE
372 ctemp=0.0
373 do i = 1, nn
374 ctemp = ctemp + h(i)*tt(i)
375 enddo
376 call heat_get_coefficient( ctemp,imat,cc,ntab,temp,funca,funcb )
377
378 !* SET INTEGRATION WEIGHT
379 valx = cc(1)*wgt(ig1)*wgt(ig2)*wgt(ig3)*det
380 valy = cc(2)*wgt(ig1)*wgt(ig2)*wgt(ig3)*det
381 ss( 1) = ss( 1) + dtdx(1)*dtdx(1)*valx &
382 + dtdy(1)*dtdy(1)*valy
383 ss( 2) = ss( 2) + dtdx(1)*dtdx(2)*valx &
384 + dtdy(1)*dtdy(2)*valy
385 ss( 3) = ss( 3) + dtdx(1)*dtdx(3)*valx &
386 + dtdy(1)*dtdy(3)*valy
387 ss( 4) = ss( 4) + dtdx(1)*dtdx(4)*valx &
388 + dtdy(1)*dtdy(4)*valy
389 ss( 6) = ss( 6) + dtdx(2)*dtdx(2)*valx &
390 + dtdy(2)*dtdy(2)*valy
391 ss( 7) = ss( 7) + dtdx(2)*dtdx(3)*valx &
392 + dtdy(2)*dtdy(3)*valy
393 ss( 8) = ss( 8) + dtdx(2)*dtdx(4)*valx &
394 + dtdy(2)*dtdy(4)*valy
395 ss(11) = ss(11) + dtdx(3)*dtdx(3)*valx &
396 + dtdy(3)*dtdy(3)*valy
397 ss(12) = ss(12) + dtdx(3)*dtdx(4)*valx &
398 + dtdy(3)*dtdy(4)*valy
399 ss(16) = ss(16) + dtdx(4)*dtdx(4)*valx &
400 + dtdy(4)*dtdy(4)*valy
401 enddo
402 enddo
403 enddo
404
405 ss( 3) = ss( 3) + ss( 4)
406 ss( 7) = ss( 7) + ss( 8)
407 ss(11) = ss(11) + ss(16) + 2.0*ss(12)
408 ss( 4) = 0.0d0
409 ss( 8) = 0.0d0
410 ss(12) = 0.0d0
411 ss(16) = 0.0d0
412
413 ss( 5) = ss( 2)
414 ss( 9) = ss( 3)
415 ss(10) = ss( 7)
416 ss(13) = ss( 4)
417 ss(14) = ss( 8)
418 ss(15) = ss(12)
419 end subroutine heat_conductivity_shell_731
420
421 subroutine heat_conductivity_shell_741(etype, nn, ecoord, TT, IMAT, thick, SS, stiff, &
422 ntab, temp, funcA, funcB)
423 use mmechgauss
424 use m_matmatrix
425 use elementinfo
426 implicit none
427 !type(tGaussStatus), intent(in) :: gausses(:) !< status of qudrature points
428 integer(kind=kint), intent(in) :: etype
429 integer(kind=kint), intent(in) :: nn
430 real(kind=kreal), intent(in) :: ecoord(3,nn)
431 real(kind=kreal), intent(out) :: stiff(:,:)
432 real(kind=kreal), intent(inout) :: ss(:)
433 real(kind=kreal), intent(in) :: tt(nn)
434 type(tmaterial), pointer :: matl
435 integer(kind=kint) :: i, j, lx, ly, imat, ntab
436 integer(kind=kint) :: ig1, ig2, ig3, inod
437 real(kind=kreal) :: ti, valx, valy, valz, var
438 real(kind=kreal) :: xj11, xj12, xj13, xj21, xj22, xj23, xj31, xj32, xj33
439 real(kind=kreal) :: naturalcoord(2), xsum
440 real(kind=kreal) :: func(nn), thick, temp_i
441 real(kind=kreal) :: det, wg, rho, diag_stiff, total_stiff
442 real(kind=kreal) :: d(1,1), n(1, nn), dn(1, nn)
443 real(kind=kreal) :: gderiv(nn,2)
444 real(kind=kreal) :: cc(3)
445 real(kind=kreal) :: temp(ntab), funca(ntab+1), funcb(ntab+1)
446 real(kind=kreal) :: xg(2), ri, si, rp, sp, rm, sm, hr(4), hs(4)
447 real(kind=kreal) :: xr, xs, yr, ys, zr, zs
448 real(kind=kreal) :: h(4), x(4), y(4), z(4)
449 real(kind=kreal) :: wgt(2), ctemp, dum
450 real(kind=kreal) :: cod(3, 4)
451 real(kind=kreal) :: g1(3), g2(3), g3(3), e1(3), e2(3), e3(3), ref(3)
452 real(kind=kreal) :: en(3, 4), the(3, 3), amat(3, 3), bv(3), wk(3)
453 real(kind=kreal) :: dtdx(4), dtdy(4)
454 data xg/-0.5773502691896d0,0.5773502691896d0/
455 data wgt/1.0d0,1.0d0/
456
457 do i = 1, nn
458 cod(1,i) = ecoord(1,i)
459 cod(2,i) = ecoord(2,i)
460 cod(3,i) = ecoord(3,i)
461 enddo
462
463 !* SET REFFRENSE VECTOR TO DETERMINE LOCAL COORDINATE SYSTEM
464 ref(1) = 0.25*( cod(1,2) + cod(1,3) - cod(1,1) - cod(1,4) )
465 ref(2) = 0.25*( cod(2,2) + cod(2,3) - cod(2,1) - cod(2,4) )
466 ref(3) = 0.25*( cod(3,2) + cod(3,3) - cod(3,1) - cod(3,4) )
467
468 do i = 1, nn
469 if ( i .EQ. 1 ) then
470 g1(1) = -0.5*cod(1,1) + 0.5*cod(1,2)
471 g1(2) = -0.5*cod(2,1) + 0.5*cod(2,2)
472 g1(3) = -0.5*cod(3,1) + 0.5*cod(3,2)
473 g2(1) = -0.5*cod(1,1) + 0.5*cod(1,4)
474 g2(2) = -0.5*cod(2,1) + 0.5*cod(2,4)
475 g2(3) = -0.5*cod(3,1) + 0.5*cod(3,4)
476 else if( i .EQ. 2 ) then
477 g1(1) = -0.5*cod(1,1) + 0.5*cod(1,2)
478 g1(2) = -0.5*cod(2,1) + 0.5*cod(2,2)
479 g1(3) = -0.5*cod(3,1) + 0.5*cod(3,2)
480 g2(1) = -0.5*cod(1,2) + 0.5*cod(1,3)
481 g2(2) = -0.5*cod(2,2) + 0.5*cod(2,3)
482 g2(3) = -0.5*cod(3,2) + 0.5*cod(3,3)
483 else if( i .EQ. 3 ) then
484 g1(1) = 0.5*cod(1,3) - 0.5*cod(1,4)
485 g1(2) = 0.5*cod(2,3) - 0.5*cod(2,4)
486 g1(3) = 0.5*cod(3,3) - 0.5*cod(3,4)
487 g2(1) = -0.5*cod(1,2) + 0.5*cod(1,3)
488 g2(2) = -0.5*cod(2,2) + 0.5*cod(2,3)
489 g2(3) = -0.5*cod(3,2) + 0.5*cod(3,3)
490 else if( i .EQ. 4 ) then
491 g1(1) = 0.5*cod(1,3) - 0.5*cod(1,4)
492 g1(2) = 0.5*cod(2,3) - 0.5*cod(2,4)
493 g1(3) = 0.5*cod(3,3) - 0.5*cod(3,4)
494 g2(1) = -0.5*cod(1,1) + 0.5*cod(1,4)
495 g2(2) = -0.5*cod(2,1) + 0.5*cod(2,4)
496 g2(3) = -0.5*cod(3,1) + 0.5*cod(3,4)
497 end if
498
499 !* G3()=G1() X G2()
500 g3(1) = g1(2)*g2(3) - g1(3)*g2(2)
501 g3(2) = g1(3)*g2(1) - g1(1)*g2(3)
502 g3(3) = g1(1)*g2(2) - g1(2)*g2(1)
503
504 !* SET BASE VECTOR IN LOCAL CARTESIAN COORDINATE
505 xsum = dsqrt( g3(1)**2 + g3(2)**2 + g3(3)**2 )
506 e3(1) = g3(1) / xsum
507 e3(2) = g3(2) / xsum
508 e3(3) = g3(3) / xsum
509 en(1,i) = e3(1)
510 en(2,i) = e3(2)
511 en(3,i) = e3(3)
512 enddo
513
514 !* LOOP FOR GAUSS INTEGRATION POINT
515 do ig3 = 1, 2
516 ti = xg(ig3)
517 do ig2 = 1, 2
518 si = xg(ig2)
519 do ig1 = 1, 2
520 ri = xg(ig1)
521 rp = 1.0 + ri
522 sp = 1.0 + si
523 rm = 1.0 - ri
524 sm = 1.0 - si
525
526 h(1) = 0.25*rm*sm
527 h(2) = 0.25*rp*sm
528 h(3) = 0.25*rp*sp
529 h(4) = 0.25*rm*sp
530
531 hr(1) = -0.25*sm
532 hr(2) = 0.25*sm
533 hr(3) = 0.25*sp
534 hr(4) = -0.25*sp
535
536 hs(1) = -0.25*rm
537 hs(2) = -0.25*rp
538 hs(3) = 0.25*rp
539 hs(4) = 0.25*rm
540
541 !* COVARIANT BASE VECTOR AT A GAUSS INTEGRARION POINT
542 do i = 1, 3
543 g1(i) = 0.0
544 g2(i) = 0.0
545 g3(i) = 0.0
546 do inod = 1, nn
547 var = cod(i,inod) + thick*0.5*ti*en(i,inod)
548 g1(i) = g1(i) + hr(inod)*var
549 g2(i) = g2(i) + hs(inod)*var
550 g3(i) = g3(i) + thick*0.5*h(inod)*en(i,inod)
551 enddo
552 enddo
553
554 !*JACOBI MATRIX
555 xj11 = g1(1)
556 xj12 = g1(2)
557 xj13 = g1(3)
558 xj21 = g2(1)
559 xj22 = g2(2)
560 xj23 = g2(3)
561 xj31 = g3(1)
562 xj32 = g3(2)
563 xj33 = g3(3)
564
565 !*DETERMINANT OF JACOBIAN
566 det = xj11*xj22*xj33 &
567 + xj12*xj23*xj31 &
568 + xj13*xj21*xj32 &
569 - xj13*xj22*xj31 &
570 - xj12*xj21*xj33 &
571 - xj11*xj23*xj32
572
573 !* INVERSION OF JACOBIAN
574 dum = 1.0 / det
575 amat(1,1) = dum*( xj22*xj33 - xj23*xj32 )
576 amat(2,1) = dum*( -xj21*xj33 + xj23*xj31 )
577 amat(3,1) = dum*( xj21*xj32 - xj22*xj31 )
578 amat(1,2) = dum*( -xj12*xj33 + xj13*xj32 )
579 amat(2,2) = dum*( xj11*xj33 - xj13*xj31 )
580 amat(3,2) = dum*( -xj11*xj32 + xj12*xj31 )
581 amat(1,3) = dum*( xj12*xj23 - xj13*xj22 )
582 amat(2,3) = dum*( -xj11*xj23 + xj13*xj21 )
583 amat(3,3) = dum*( xj11*xj22 - xj12*xj21 )
584
585 !* SET BASE VECTOR IN LOCAL CARTESIAN COORDINATE
586 xsum = dsqrt( g3(1)**2 + g3(2)**2 + g3(3)**2 )
587 e3(1) = g3(1) / xsum
588 e3(2) = g3(2) / xsum
589 e3(3) = g3(3) / xsum
590 e2(1) = -ref(2)*e3(3) + ref(3)*e3(2)
591 e2(2) = -ref(3)*e3(1) + ref(1)*e3(3)
592 e2(3) = -ref(1)*e3(2) + ref(2)*e3(1)
593 e1(1) = -e3(2)*e2(3) + e3(3)*e2(2)
594 e1(2) = -e3(3)*e2(1) + e3(1)*e2(3)
595 e1(3) = -e3(1)*e2(2) + e3(2)*e2(1)
596 xsum = dsqrt( e2(1)**2 + e2(2)**2 + e2(3)**2 )
597
598 if ( xsum .GT. 1.e-15 ) then
599 e2(1) = e2(1) / xsum
600 e2(2) = e2(2) / xsum
601 e2(3) = e2(3) / xsum
602 e1(1) = -e3(2)*e2(3) + e3(3)*e2(2)
603 e1(2) = -e3(3)*e2(1) + e3(1)*e2(3)
604 e1(3) = -e3(1)*e2(2) + e3(2)*e2(1)
605 xsum = dsqrt( e1(1)**2 + e1(2)**2 + e1(3)**2 )
606
607 e1(1) = e1(1) / xsum
608 e1(2) = e1(2) / xsum
609 e1(3) = e1(3) / xsum
610 else
611 e1(1) = 0.d0
612 e1(2) = 0.d0
613 e1(3) = -1.d0
614 e2(1) = 0.d0
615 e2(2) = 1.d0
616 e2(3) = 0.d0
617 endif
618 the(1,1) = e1(1)
619 the(1,2) = e1(2)
620 the(1,3) = e1(3)
621 the(2,1) = e2(1)
622 the(2,2) = e2(2)
623 the(2,3) = e2(3)
624 the(3,1) = e3(1)
625 the(3,2) = e3(2)
626 the(3,3) = e3(3)
627 do i = 1, nn
628 bv(1) = hr(i)
629 bv(2) = hs(i)
630 bv(3) = 0.0
631 wk(1) = amat(1,1)*bv(1) &
632 + amat(1,2)*bv(2) &
633 + amat(1,3)*bv(3)
634 wk(2) = amat(2,1)*bv(1) &
635 + amat(2,2)*bv(2) &
636 + amat(2,3)*bv(3)
637 wk(3) = amat(3,1)*bv(1) &
638 + amat(3,2)*bv(2) &
639 + amat(3,3)*bv(3)
640 bv(1) = the(1,1)*wk(1) &
641 + the(1,2)*wk(2) &
642 + the(1,3)*wk(3)
643 bv(2) = the(2,1)*wk(1) &
644 + the(2,2)*wk(2) &
645 + the(2,3)*wk(3)
646 bv(3) = the(3,1)*wk(1) &
647 + the(3,2)*wk(2) &
648 + the(3,3)*wk(3)
649 dtdx(i) = bv(1)
650 dtdy(i) = bv(2)
651 enddo
652
653 !*CONDUCTIVITY AT CURRENT TEMPERATURE
654 ctemp=0.0
655 do i = 1, nn
656 ctemp = ctemp + h(i)*tt(i)
657 enddo
658 call heat_get_coefficient( ctemp,imat,cc,ntab,temp,funca,funcb )
659
660 !* SET INTEGRATION WEIGHT
661 valx = cc(1)*wgt(ig1)*wgt(ig2)*wgt(ig3)*det
662 valy = cc(2)*wgt(ig1)*wgt(ig2)*wgt(ig3)*det
663 ss( 1) = ss( 1) + dtdx(1)*dtdx(1)*valx &
664 + dtdy(1)*dtdy(1)*valy
665 ss( 2) = ss( 2) + dtdx(1)*dtdx(2)*valx &
666 + dtdy(1)*dtdy(2)*valy
667 ss( 3) = ss( 3) + dtdx(1)*dtdx(3)*valx &
668 + dtdy(1)*dtdy(3)*valy
669 ss( 4) = ss( 4) + dtdx(1)*dtdx(4)*valx &
670 + dtdy(1)*dtdy(4)*valy
671 ss( 6) = ss( 6) + dtdx(2)*dtdx(2)*valx &
672 + dtdy(2)*dtdy(2)*valy
673 ss( 7) = ss( 7) + dtdx(2)*dtdx(3)*valx &
674 + dtdy(2)*dtdy(3)*valy
675 ss( 8) = ss( 8) + dtdx(2)*dtdx(4)*valx &
676 + dtdy(2)*dtdy(4)*valy
677 ss(11) = ss(11) + dtdx(3)*dtdx(3)*valx &
678 + dtdy(3)*dtdy(3)*valy
679 ss(12) = ss(12) + dtdx(3)*dtdx(4)*valx &
680 + dtdy(3)*dtdy(4)*valy
681 ss(16) = ss(16) + dtdx(4)*dtdx(4)*valx &
682 + dtdy(4)*dtdy(4)*valy
683 enddo
684 enddo
685 enddo
686 ss( 5) = ss( 2)
687 ss( 9) = ss( 3)
688 ss(10) = ss( 7)
689 ss(13) = ss( 4)
690 ss(14) = ss( 8)
691 ss(15) = ss(12)
692 end subroutine heat_conductivity_shell_741
693
694 subroutine heat_conductivity_541(NN,ecoord,TEMP,TZERO,THICK,HH,RR1,RR2,SS,stiff )
695 use hecmw
696 implicit real(kind=kreal) (a - h, o - z)
697 real(kind=kreal), intent(in) :: ecoord(3,nn)
698 real(kind=kreal), intent(out) :: stiff(:,:)
699 dimension xxx(nn), yyy(nn), zzz(nn), temp(nn), ss(nn*nn)
700 dimension xx(4), yy(4), zz(4)
701
702 xx(1)=ecoord(1,1)
703 xx(2)=ecoord(1,2)
704 xx(3)=ecoord(1,3)
705 xx(4)=ecoord(1,4)
706 yy(1)=ecoord(2,1)
707 yy(2)=ecoord(2,2)
708 yy(3)=ecoord(2,3)
709 yy(4)=ecoord(2,4)
710 zz(1)=ecoord(3,1)
711 zz(2)=ecoord(3,2)
712 zz(3)=ecoord(3,3)
713 zz(4)=ecoord(3,4)
714 call heat_get_area(xx, yy, zz, sa)
715
716 xx(1)=ecoord(1,5)
717 xx(2)=ecoord(1,6)
718 xx(3)=ecoord(1,7)
719 xx(4)=ecoord(1,8)
720 yy(1)=ecoord(2,5)
721 yy(2)=ecoord(2,6)
722 yy(3)=ecoord(2,7)
723 yy(4)=ecoord(2,8)
724 zz(1)=ecoord(3,5)
725 zz(2)=ecoord(3,6)
726 zz(3)=ecoord(3,7)
727 zz(4)=ecoord(3,8)
728 call heat_get_area(xx, yy, zz, sb)
729
730 t1z=temp(1)-tzero
731 t2z=temp(2)-tzero
732 t3z=temp(3)-tzero
733 t4z=temp(4)-tzero
734 t5z=temp(5)-tzero
735 t6z=temp(6)-tzero
736 t7z=temp(7)-tzero
737 t8z=temp(8)-tzero
738
739 rrr1 = rr1**0.25
740 rrr2 = rr2**0.25
741
742 ha1= (( rrr1 * t1z )**2 + ( rrr2 * t5z )**2 ) &
743 & * ( rrr1 * t1z + rrr2 * t5z ) * rrr1
744 ha2= (( rrr1 * t2z )**2 + ( rrr2 * t6z )**2 ) &
745 & * ( rrr1 * t2z + rrr2 * t6z ) * rrr1
746 ha3= (( rrr1 * t3z )**2 + ( rrr2 * t7z )**2 ) &
747 & * ( rrr1 * t3z + rrr2 * t7z ) * rrr1
748 ha4= (( rrr1 * t4z )**2 + ( rrr2 * t8z )**2 ) &
749 & * ( rrr1 * t4z + rrr2 * t8z ) * rrr1
750
751 hb1= (( rrr1 * t1z )**2 + ( rrr2 * t5z )**2 ) &
752 & * ( rrr1 * t1z + rrr2 * t5z ) * rrr2
753 hb2= (( rrr1 * t2z )**2 + ( rrr2 * t6z )**2 ) &
754 & * ( rrr1 * t2z + rrr2 * t6z ) * rrr2
755 hb3= (( rrr1 * t3z )**2 + ( rrr2 * t7z )**2 ) &
756 & * ( rrr1 * t3z + rrr2 * t7z ) * rrr2
757 hb4= (( rrr1 * t4z )**2 + ( rrr2 * t8z )**2 ) &
758 & * ( rrr1 * t4z + rrr2 * t8z ) * rrr2
759
760 hhh = hh / thick
761
762 ss( 1) = ( hhh + ha1 ) * sa * 0.25
763 ss(10) = ( hhh + ha2 ) * sa * 0.25
764 ss(19) = ( hhh + ha3 ) * sa * 0.25
765 ss(28) = ( hhh + ha4 ) * sa * 0.25
766
767 ss(37) = ( hhh + hb1 ) * sb * 0.25
768 ss(46) = ( hhh + hb2 ) * sb * 0.25
769 ss(55) = ( hhh + hb3 ) * sb * 0.25
770 ss(64) = ( hhh + hb4 ) * sb * 0.25
771
772 sm = ( sa + sb ) * 0.5
773 hh1 = ( ha1 + hb1 ) * 0.5
774 hh2 = ( ha2 + hb2 ) * 0.5
775 hh3 = ( ha3 + hb3 ) * 0.5
776 hh4 = ( ha4 + hb4 ) * 0.5
777
778 ss(33) = -( hhh + hh1 ) * sm * 0.25
779 ss(42) = -( hhh + hh2 ) * sm * 0.25
780 ss(51) = -( hhh + hh3 ) * sm * 0.25
781 ss(60) = -( hhh + hh4 ) * sm * 0.25
782
783 ss( 5) = ss(33)
784 ss(14) = ss(42)
785 ss(23) = ss(51)
786 ss(32) = ss(60)
787
788 !C 1 ( 1- 1) 2 3 4 5 6 7 8
789 !C 2 3 9 (10- 3) 11 12 13 14 15 16
790 !C 4 5 6 17 18 (19- 6) 20 21 22 23 24
791 !C 7 8 9 10 ---> 25 26 27 (28-10) 29 30 31 32
792 !C 11 12 13 14 15 ---> 33 34 35 36 (37-15) 38 39 40
793 !C 16 17 18 19 20 21 41 42 43 44 45 (46-21) 47 48
794 !C 22 23 24 25 26 27 28 49 50 51 52 53 54 (55-28) 56
795 !C 29 30 31 32 33 34 35 36 57 58 59 60 61 62 63 (64-36)
796 end subroutine heat_conductivity_541
797
798 subroutine heat_get_area ( XX,YY,ZZ,AA )
799 use hecmw
800 implicit real(kind=kreal)(a-h,o-z)
801 dimension xx(4),yy(4),zz(4)
802 dimension xg(2),h(4),hr(4),hs(4)
803
804 pi=4.0*atan(1.0)
805
806 xg(1) =-0.5773502691896258d0
807 xg(2) =-xg(1)
808 aa=0.0
809
810 do lx=1,2
811 ri=xg(lx)
812 do ly=1,2
813 si=xg(ly)
814 rp=1.0+ri
815 sp=1.0+si
816 rm=1.0-ri
817 sm=1.0-si
818 !*INTERPOLATION FUNCTION
819 h(1)=0.25*rp*sp
820 h(2)=0.25*rm*sp
821 h(3)=0.25*rm*sm
822 h(4)=0.25*rp*sm
823 !*DERIVATIVE OF INTERPOLATION FUNCTION
824 !* FOR R-COORDINATE
825 hr(1)= .25*sp
826 hr(2)=-.25*sp
827 hr(3)=-.25*sm
828 hr(4)= .25*sm
829 !* FOR S-COORDINATE
830 hs(1)= .25*rp
831 hs(2)= .25*rm
832 hs(3)=-.25*rm
833 hs(4)=-.25*rp
834 !*JACOBI MATRIX
835 xr=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4)
836 xs=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4)
837 yr=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4)
838 ys=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4)
839 zr=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4)
840 zs=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4)
841
842 det=(yr*zs-zr*ys)**2+(zr*xs-xr*zs)**2+(xr*ys-yr*xs)**2
843 det=sqrt(det)
844 aa=aa+det
845 enddo
846 enddo
847 end subroutine heat_get_area
848
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 module contains subroutines used in 3d eigen analysis for.
Definition: dynamic_mass.f90:6
subroutine heat_conductivity_shell_741(etype, nn, ecoord, tt, imat, thick, ss, stiff, ntab, temp, funca, funcb)
subroutine heat_conductivity_c3(etype, nn, ecoord, temperature, imat, stiff, ntab, temp, funca, funcb)
subroutine heat_get_coefficient(tpoi, imat, coef, ntab, temp, funca, funcb)
subroutine heat_conductivity_c1(etype, nn, ecoord, temperature, imat, surf, stiff, ntab, temp, funca, funcb)
subroutine heat_conductivity_c2(etype, nn, ecoord, temperature, imat, thick, stiff, ntab, temp, funca, funcb)
subroutine heat_conductivity_541(nn, ecoord, temp, tzero, thick, hh, rr1, rr2, ss, stiff)
subroutine heat_conductivity_shell_731(etype, nn, ecoord, tt, imat, thick, ss, stiff, ntab, temp, funca, funcb)
CALCULATION 4 NODE SHELL ELEMENT.
subroutine heat_get_area(xx, yy, zz, aa)
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6