FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
6 contains
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
166  use m_dynamic_mass
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 
849 end module m_heat_lib_conductivity
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_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_c3(etype, nn, ecoord, temperature, IMAT, stiff, ntab, temp, funcA, funcB)
subroutine heat_conductivity_c2(etype, nn, ecoord, temperature, IMAT, thick, stiff, ntab, temp, funcA, funcB)
subroutine heat_conductivity_shell_731(etype, nn, ecoord, TT, IMAT, thick, SS, stiff, ntab, temp, funcA, funcB)
CALCULATION 4 NODE SHELL ELEMENT.
subroutine heat_conductivity_541(NN, ecoord, TEMP, TZERO, THICK, HH, RR1, RR2, SS, stiff)
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