FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
heat_LIB_FILM.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 module m_heat_lib_film
8 contains
9  !***********************************************************************
10  ! FILM_231(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
11  ! FILM_232(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
12  ! FILM_241(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
13  ! FILM_242(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
14  ! FILM_341(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
15  ! FILM_342(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
16  ! FILM_351(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
17  ! FILM_352(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
18  ! FILM_361(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
19  ! FILM_362(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
20  ! FILM_731(NN,XX,YY,ZZ,LTYPE,HH,SINK,TERM1,TERM2)
21  ! FILM_741(NN,XX,YY,ZZ,LTYPE,HH,SINK,TERM1,TERM2)
22  !----------------------------------------------------------------------*
23  subroutine heat_film_231(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
24  !----------------------------------------------------------------------*
25  !**
26  !** SET 231 FILM
27  !**
28  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
29  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
30  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
31  use hecmw
32  implicit none
33  ! I/F VARIABLES
34  integer(kind=kint) NN, LTYPE, MM
35  integer(kind=kint) NOD(MM)
36  real(kind=kreal) thick, hh, sink
37  real(kind=kreal) xx(nn),yy(nn),zz(nn),term1(mm*mm),term2(mm)
38  ! LOCAL VARIABLES
39  integer(kind=kint) I, IC, IP, JP, LX
40  real(kind=kreal) ri, gx, gy, xsum
41  real(kind=kreal) xg(2),wgt(2),h(2),hr(2)
42  data wgt/1.0,1.0/
43  data xg/-0.5773502691896, 0.5773502691896/
44  !
45  if(ltype.EQ.1) then
46  nod(1)=1
47  nod(2)=2
48  else if(ltype.EQ.2) then
49  nod(1)=2
50  nod(2)=3
51  else if(ltype.EQ.3) then
52  nod(1)=3
53  nod(2)=1
54  endif
55  !
56  ic = 0
57  do ip = 1, 2
58  term2(ip) = 0.0
59  do jp = 1, 2
60  ic = ic + 1
61  term1(ic) = 0.0
62  !** INTEGRATION OVER SURFACE
63  do lx=1,2
64  ri=xg(lx)
65  h(1)=0.5*(1.0-ri)
66  h(2)=0.5*(1.0+ri)
67  hr(1)=-0.5
68  hr(2)= 0.5
69  gx=0.0
70  gy=0.0
71  do i=1,2
72  gx=gx+hr(i)*xx(nod(i))
73  gy=gy+hr(i)*yy(nod(i))
74  enddo
75  xsum = dsqrt( gx*gx+gy*gy )*thick
76  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
77  if( ip.EQ.jp ) then
78  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
79  endif
80  enddo
81  enddo
82  enddo
83  !
84  return
85 
86  end subroutine heat_film_231
87  !----------------------------------------------------------------------*
88  subroutine heat_film_232(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
89  !----------------------------------------------------------------------*
90  !**
91  !** SET 232 FILM
92  !**
93  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
94  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
95  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
96  use hecmw
97  implicit none
98  ! I/F VARIABLES
99  integer(kind=kint) :: NN, LTYPE, MM
100  integer(kind=kint) :: NOD(MM)
101  real(kind=kreal) :: thick, hh, sink
102  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
103  ! LOCAL VARIABLES
104  integer(kind=kint) :: I, IC, IP, JP, LX
105  real(kind=kreal) :: ri, gx, gy, xsum
106  real(kind=kreal) :: xg(3), wgt(3), h(3), hr(3)
107  !*************************
108  ! GAUSS INTEGRATION POINT
109  !*************************
110  xg(1) = -0.7745966692
111  xg(2) = 0.0
112  xg(3) = 0.7745966692
113  wgt(1) = 0.5555555555
114  wgt(2) = 0.8888888888
115  wgt(3) = 0.5555555555
116 
117  if(ltype.EQ.1) then
118  nod(1)=1
119  nod(2)=2
120  nod(3)=4
121  else if(ltype.EQ.2) then
122  nod(1)=2
123  nod(2)=3
124  nod(3)=5
125  else if(ltype.EQ.3) then
126  nod(1)=3
127  nod(2)=1
128  nod(3)=6
129  endif
130 
131  ic = 0
132  do ip = 1, 3
133  term2(ip) = 0.0
134  do jp = 1, 3
135  ic = ic + 1
136  term1(ic) = 0.0
137  !** INTEGRATION OVER SURFACE
138  do lx=1,3
139  ri=xg(lx)
140  h(1)=-0.5*(1.0-ri)*ri
141  h(2)= 0.5*(1.0+ri)*ri
142  h(3)=1.0-ri*ri
143  hr(1)= ri-0.5
144  hr(2)= ri+0.5
145  hr(3)=-2.0*ri
146  gx=0.0
147  gy=0.0
148  do i=1,3
149  gx=gx+hr(i)*xx(nod(i))
150  gy=gy+hr(i)*yy(nod(i))
151  enddo
152  xsum = dsqrt( gx*gx+gy*gy )*thick
153  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
154  if( ip.EQ.jp ) then
155  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
156  endif
157  enddo
158  enddo
159  enddo
160 
161  return
162 
163  end subroutine heat_film_232
164  !----------------------------------------------------------------------*
165  subroutine heat_film_241(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
166  !----------------------------------------------------------------------*
167  !**
168  !** SET 241 FILM
169  !**
170  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
171  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
172  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
173  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
174  use hecmw
175  implicit none
176  ! I/F VARIABLES
177  integer(kind=kint) :: NN, LTYPE, MM
178  integer(kind=kint) :: NOD(MM)
179  real(kind=kreal) :: thick, hh, sink
180  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
181  ! LOCAL VARIABLES
182  integer(kind=kint) :: I, IC, IP, JP, LX
183  real(kind=kreal) :: ri, gx, gy, xsum
184  real(kind=kreal) :: xg(2), wgt(2), h(2), hr(2)
185  data wgt/1.0, 1.0/
186  data xg/-0.5773502691896, 0.5773502691896/
187 
188  if(ltype.EQ.1) then
189  nod(1)=1
190  nod(2)=2
191  else if(ltype.EQ.2) then
192  nod(1)=2
193  nod(2)=3
194  else if(ltype.EQ.3) then
195  nod(1)=3
196  nod(2)=4
197  else if(ltype.EQ.4) then
198  nod(1)=4
199  nod(2)=1
200  endif
201 
202  ic = 0
203  do ip = 1, 2
204  term2(ip) = 0.0
205  do jp = 1, 2
206  ic = ic + 1
207  term1(ic) = 0.0
208  !** INTEGRATION OVER SURFACE
209  do lx=1,2
210  ri=xg(lx)
211  h(1)=0.5*(1.0-ri)
212  h(2)=0.5*(1.0+ri)
213  hr(1)=-0.5
214  hr(2)= 0.5
215  gx=0.0
216  gy=0.0
217  do i=1,2
218  gx=gx+hr(i)*xx(nod(i))
219  gy=gy+hr(i)*yy(nod(i))
220  enddo
221  xsum = dsqrt( gx*gx+gy*gy )*thick
222  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
223  if( ip.EQ.jp ) then
224  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
225  endif
226  enddo
227  enddo
228  enddo
229 
230  return
231 
232  end subroutine heat_film_241
233  !----------------------------------------------------------------------*
234  subroutine heat_film_242(NN,XX,YY,ZZ,THICK,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
235  !----------------------------------------------------------------------*
236  !**
237  !** SET 242 FILM
238  !**
239  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
240  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
241  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
242  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
243  use hecmw
244  implicit none
245  ! I/F VARIABLES
246  integer(kind=kint) :: NN, LTYPE, MM
247  integer(kind=kint) :: NOD(MM)
248  real(kind=kreal) :: thick, hh, sink
249  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
250  ! LOCAL VARIABLES
251  integer(kind=kint) :: I, IC, IP, JP, LX
252  real(kind=kreal) :: ri, gx, gy, xsum
253  real(kind=kreal) :: xg(3), wgt(3), h(3), hr(3)
254  !*************************
255  ! GAUSS INTEGRATION POINT
256  !*************************
257  xg(1) = -0.7745966692
258  xg(2) = 0.0
259  xg(3) = 0.7745966692
260  wgt(1) = 0.5555555555
261  wgt(2) = 0.8888888888
262  wgt(3) = 0.5555555555
263 
264  if(ltype.EQ.1) then
265  nod(1)=1
266  nod(2)=2
267  nod(3)=5
268  else if(ltype.EQ.2) then
269  nod(1)=2
270  nod(2)=3
271  nod(3)=6
272  else if(ltype.EQ.3) then
273  nod(1)=3
274  nod(2)=4
275  nod(3)=7
276  else if(ltype.EQ.4) then
277  nod(1)=4
278  nod(2)=1
279  nod(3)=8
280  endif
281 
282  ic = 0
283  do ip = 1, 3
284  term2(ip) = 0.0
285  do jp = 1, 3
286  ic = ic + 1
287  term1(ic) = 0.0
288  !** INTEGRATION OVER SURFACE
289  do lx=1,3
290  ri=xg(lx)
291  h(1)=-0.5*(1.0-ri)*ri
292  h(2)= 0.5*(1.0+ri)*ri
293  h(3)=1.0-ri*ri
294  hr(1)= ri-0.5
295  hr(2)= ri+0.5
296  hr(3)=-2.0*ri
297  gx=0.0
298  gy=0.0
299  do i=1,3
300  gx=gx+hr(i)*xx(nod(i))
301  gy=gy+hr(i)*yy(nod(i))
302  enddo
303  xsum = dsqrt( gx*gx+gy*gy )*thick
304  term1(ic) = term1(ic) - xsum*wgt(lx)*h(ip)*h(jp)*hh
305  if( ip.EQ.jp ) then
306  term2(ip) = term2(ip) - xsum*wgt(lx)*h(jp)*hh*sink
307  endif
308  enddo
309  enddo
310  enddo
311 
312  return
313 
314  end subroutine heat_film_242
315  !----------------------------------------------------------------------*
316  subroutine heat_film_341(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
317  !----------------------------------------------------------------------*
318  !**
319  !** SET 341 FILM
320  !**
321  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
322  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
323  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
324  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
325  use hecmw
326  implicit none
327  ! I/F VARIABLES
328  integer(kind=kint) :: NN, LTYPE, MM
329  integer(kind=kint) :: NOD(MM)
330  real(kind=kreal) :: hh, sink
331  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
332  ! LOCAL VARIABLES
333  integer(kind=kint) :: IC, IP, JP
334  real(kind=kreal) :: ax, ay, az, bx, by, bz, aa
335 
336  !
337  if ( ltype.EQ.1 ) then
338  nod(1) = 1
339  nod(2) = 2
340  nod(3) = 3
341  else if( ltype.EQ.2 ) then
342  nod(1) = 4
343  nod(2) = 2
344  nod(3) = 1
345  else if( ltype.EQ.3 ) then
346  nod(1) = 4
347  nod(2) = 3
348  nod(3) = 2
349  else if( ltype.EQ.4 ) then
350  nod(1) = 4
351  nod(2) = 1
352  nod(3) = 3
353  endif
354  !
355  ax = xx( nod(2) ) - xx( nod(1) )
356  ay = yy( nod(2) ) - yy( nod(1) )
357  az = zz( nod(2) ) - zz( nod(1) )
358  bx = xx( nod(3) ) - xx( nod(1) )
359  by = yy( nod(3) ) - yy( nod(1) )
360  bz = zz( nod(3) ) - zz( nod(1) )
361  !
362  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
363  !
364  ic = 0
365  do ip = 1, 3
366  term2(ip) = - aa*hh*sink
367  do jp = 1, 3
368  ic = ic + 1
369  term1(ic) = - aa*hh/3.0
370  enddo
371  enddo
372  !
373  return
374 
375  end subroutine heat_film_341
376  !----------------------------------------------------------------------*
377  subroutine heat_film_342(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
378  !----------------------------------------------------------------------*
379  !**
380  !** SET 342 FILM
381  !**
382  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
383  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
384  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
385  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
386  use hecmw
387  implicit none
388  ! I/F VARIABLES
389  integer(kind=kint) :: NN, LTYPE, MM
390  integer(kind=kint) :: NOD(MM)
391  real(kind=kreal) :: hh, sink
392  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
393  ! LOCAL VARIABLES
394  integer(kind=kint) :: I, IC, IP, JP, L1, L2
395  real(kind=kreal) :: x1, x2, x3, xl1, xl2
396  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
397  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
398  real(kind=kreal) :: xsum, det, wg
399  real(kind=kreal) :: xg(3), wgt(3), h(6), hl1(6), hl2(6), hl3(6)
400  !*************************
401  ! GAUSS INTEGRATION POINT
402  !*************************
403  xg(1) = -0.7745966692
404  xg(2) = 0.0
405  xg(3) = 0.7745966692
406  wgt(1) = 0.5555555555
407  wgt(2) = 0.8888888888
408  wgt(3) = 0.5555555555
409  !
410  if(ltype.EQ.1) then
411  nod(1)=1
412  nod(2)=2
413  nod(3)=3
414  nod(4)=5
415  nod(5)=6
416  nod(6)=7
417  else if(ltype.EQ.2) then
418  nod(1)=4
419  nod(2)=2
420  nod(3)=1
421  nod(4)=9
422  nod(5)=5
423  nod(6)=8
424  else if(ltype.EQ.3) then
425  nod(1)=4
426  nod(2)=3
427  nod(3)=2
428  nod(4)=10
429  nod(5)=6
430  nod(6)=9
431  else if(ltype.EQ.4) then
432  nod(1)=4
433  nod(2)=1
434  nod(3)=3
435  nod(4)=8
436  nod(5)=7
437  nod(6)=10
438  endif
439  !
440  ic = 0
441  do ip = 1, 6
442  term2(ip) = 0.0
443  do jp = 1, 6
444  ic = ic + 1
445  term1(ic) = 0.0
446  ! INTEGRATION OVER SURFACE
447  do l2=1,3
448  xl2=xg(l2)
449  x2 =(xl2+1.0)*0.5
450  do l1=1,3
451  xl1=xg(l1)
452  x1=0.5*(1.0-x2)*(xl1+1.0)
453  ! INTERPOLATION FUNCTION
454  x3=1.0-x1-x2
455  h(1)= x1*(2.0*x1-1.)
456  h(2)= x2*(2.0*x2-1.)
457  h(3)= x3*(2.0*x3-1.)
458  h(4)= 4.0*x1*x2
459  h(5)= 4.0*x2*x3
460  h(6)= 4.0*x1*x3
461  ! DERIVATIVE OF INTERPOLATION FUNCTION
462  ! FOR L1-COORDINATE
463  hl1(1)=4.0*x1-1.0
464  hl1(2)= 0.0
465  hl1(3)= 0.0
466  hl1(4)= 4.0*x2
467  hl1(5)= 0.0
468  hl1(6)= 4.0*x3
469  ! FOR L2-COORDINATE
470  hl2(1)= 0.0
471  hl2(2)= 4.0*x2-1.0
472  hl2(3)= 0.0
473  hl2(4)= 4.0*x1
474  hl2(5)= 4.0*x3
475  hl2(6)= 0.0
476  ! FOR L3-COORDINATE
477  hl3(1)= 0.0
478  hl3(2)= 0.0
479  hl3(3)= 4.0*x3-1.0
480  hl3(4)= 0.0
481  hl3(5)= 4.0*x2
482  hl3(6)= 4.0*x1
483  !
484  g1x=0.0
485  g1y=0.0
486  g1z=0.0
487  g2x=0.0
488  g2y=0.0
489  g2z=0.0
490  do i=1,6
491  g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
492  g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
493  g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
494  g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
495  g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
496  g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
497  enddo
498  g3x=g1y*g2z-g1z*g2y
499  g3y=g1z*g2x-g1x*g2z
500  g3z=g1x*g2y-g1y*g2x
501  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
502  g3x=g3x/xsum
503  g3y=g3y/xsum
504  g3z=g3z/xsum
505  !JACOBI MATRIX
506  xj11=g1x
507  xj12=g1y
508  xj13=g1z
509  xj21=g2x
510  xj22=g2y
511  xj23=g2z
512  xj31=g3x
513  xj32=g3y
514  xj33=g3z
515  !DETERMINANT OF JACOBIAN
516  det=xj11*xj22*xj33 &
517  +xj12*xj23*xj31 &
518  +xj13*xj21*xj32 &
519  -xj13*xj22*xj31 &
520  -xj12*xj21*xj33 &
521  -xj11*xj23*xj32
522  !
523  wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25*hh
524  term1(ic) = term1(ic)-wg*h(ip)*h(jp)
525  if( ip.EQ.jp ) term2(ip) = term2(ip)-wg*h(jp)*sink
526  !
527  enddo
528  enddo
529  enddo
530  enddo
531  !
532  return
533 
534  end subroutine heat_film_342
535  !----------------------------------------------------------------------*
536  subroutine heat_film_351(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
537  !----------------------------------------------------------------------*
538  !**
539  !** SET 351 FILM
540  !**
541  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
542  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
543  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
544  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
545  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
546  use hecmw
547  implicit none
548  ! I/F VARIABLES
549  integer(kind=kint) :: NN, LTYPE, MM
550  integer(kind=kint) :: NOD(MM)
551  real(kind=kreal) :: hh, sink
552  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
553  ! LOCAL VARIABLES
554  integer(kind=kint) :: ISUF, I, IC, IP, JP, IG1, IG2
555  real(kind=kreal) :: ri, si, xsum
556  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
557  real(kind=kreal) :: ax, ay, az, bx, by, bz, aa
558  real(kind=kreal) :: h(6), hr(6), hs(6), ht(6)
559  real(kind=kreal) :: xg(2), wgt(2)
560  !*************************
561  ! GAUSS INTEGRATION POINT
562  !*************************
563  data xg/-0.5773502691896,0.5773502691896/
564  data wgt/1.0,1.0/
565 
566  isuf = 0
567  !
568  if( ltype.EQ.1 ) then
569  nod(1) = 1
570  nod(2) = 2
571  nod(3) = 3
572  isuf = 1
573  else if( ltype.EQ.2 ) then
574  nod(1) = 6
575  nod(2) = 5
576  nod(3) = 4
577  isuf = 1
578  else if( ltype.EQ.3 ) then
579  nod(1) = 4
580  nod(2) = 5
581  nod(3) = 2
582  nod(4) = 1
583  isuf = 2
584  else if( ltype.EQ.4 ) then
585  nod(1) = 5
586  nod(2) = 6
587  nod(3) = 3
588  nod(4) = 2
589  isuf = 2
590  else if( ltype.EQ.5 ) then
591  nod(1) = 6
592  nod(2) = 4
593  nod(3) = 1
594  nod(4) = 3
595  isuf = 2
596  endif
597  !
598  if( isuf.EQ.2 ) then
599 
600  ic = 0
601  do ip = 1, 4
602  term2(ip) = 0.0
603  do jp = 1, 4
604  ic = ic + 1
605  term1(ic) = 0.0
606  !
607  do ig2=1,2
608  si=xg(ig2)
609  do ig1=1,2
610  ri=xg(ig1)
611  !
612  h(1)=0.25*(1.0-ri)*(1.0-si)
613  h(2)=0.25*(1.0+ri)*(1.0-si)
614  h(3)=0.25*(1.0+ri)*(1.0+si)
615  h(4)=0.25*(1.0-ri)*(1.0+si)
616  hr(1)=-.25*(1.0-si)
617  hr(2)= .25*(1.0-si)
618  hr(3)= .25*(1.0+si)
619  hr(4)=-.25*(1.0+si)
620  hs(1)=-.25*(1.0-ri)
621  hs(2)=-.25*(1.0+ri)
622  hs(3)= .25*(1.0+ri)
623  hs(4)= .25*(1.0-ri)
624  !
625  g1x=0.0
626  g1y=0.0
627  g1z=0.0
628  g2x=0.0
629  g2y=0.0
630  g2z=0.0
631  do i = 1, 4
632  g1x=g1x+hr(i)*xx(nod(i))
633  g1y=g1y+hr(i)*yy(nod(i))
634  g1z=g1z+hr(i)*zz(nod(i))
635  g2x=g2x+hs(i)*xx(nod(i))
636  g2y=g2y+hs(i)*yy(nod(i))
637  g2z=g2z+hs(i)*zz(nod(i))
638  enddo
639  g3x=g1y*g2z-g1z*g2y
640  g3y=g1z*g2x-g1x*g2z
641  g3z=g1x*g2y-g1y*g2x
642  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
643  !
644  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
645  if( ip.EQ.jp ) then
646  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
647  endif
648  enddo
649  enddo
650  enddo
651  enddo
652  elseif( isuf.EQ.1 ) then
653  ax = xx( nod(2) ) - xx( nod(1) )
654  ay = yy( nod(2) ) - yy( nod(1) )
655  az = zz( nod(2) ) - zz( nod(1) )
656  bx = xx( nod(3) ) - xx( nod(1) )
657  by = yy( nod(3) ) - yy( nod(1) )
658  bz = zz( nod(3) ) - zz( nod(1) )
659  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
660  !
661  ic = 0
662  do ip = 1, 3
663  term2(ip) = - aa*hh*sink
664  do jp = 1, 3
665  ic = ic + 1
666  term1(ic) = - aa*hh/3.0
667  enddo
668  enddo
669 
670  endif
671 
672  return
673 
674  end subroutine heat_film_351
675  !----------------------------------------------------------------------*
676  subroutine heat_film_352(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
677  !----------------------------------------------------------------------*
678  !**
679  !** SET 352 FILM
680  !**
681  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
682  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
683  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
684  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
685  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
686  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
687  use hecmw
688  implicit none
689  ! I/F VARIABLES
690  integer(kind=kint) :: NN, LTYPE, MM
691  integer(kind=kint) :: NOD(MM)
692  real(kind=kreal) :: hh, sink
693  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
694  ! LOCAL VARIABLES
695  integer(kind=kint) :: ISUF, I, IC, IP, JP, IG1, IG2, L1, L2
696  real(kind=kreal) :: ri, si, rp, sp, rm, sm
697  real(kind=kreal) :: xsum, det, wg
698  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
699  real(kind=kreal) :: x1, x2, x3, xl1, xl2
700  real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
701  real(kind=kreal) :: h(8), hr(8), hs(8), ht(8), hl1(6), hl2(6), hl3(6)
702  real(kind=kreal) :: xg(3), wgt(3)
703  !*************************
704  ! GAUSS INTEGRATION POINT
705  !*************************
706  xg(1) = -0.7745966692
707  xg(2) = 0.0
708  xg(3) = 0.7745966692
709  wgt(1) = 0.5555555555
710  wgt(2) = 0.8888888888
711  wgt(3) = 0.5555555555
712  isuf = 0
713  !
714  if( ltype.EQ.1 ) then
715  nod(1) = 1
716  nod(2) = 2
717  nod(3) = 3
718  nod(4) = 7
719  nod(5) = 8
720  nod(6) = 9
721  isuf = 1
722  else if( ltype.EQ.2 ) then
723  nod(1) = 6
724  nod(2) = 5
725  nod(3) = 4
726  nod(4) = 11
727  nod(5) = 10
728  nod(6) = 12
729  isuf = 1
730  else if( ltype.EQ.3 ) then
731  nod(1) = 4
732  nod(2) = 5
733  nod(3) = 2
734  nod(4) = 1
735  nod(5) = 10
736  nod(6) = 14
737  nod(7) = 7
738  nod(8) = 13
739  isuf = 2
740  else if( ltype.EQ.4 ) then
741  nod(1) = 5
742  nod(2) = 6
743  nod(3) = 3
744  nod(4) = 2
745  nod(5) = 11
746  nod(6) = 15
747  nod(7) = 8
748  nod(8) = 14
749  isuf = 2
750  else if( ltype.EQ.5 ) then
751  nod(1) = 6
752  nod(2) = 4
753  nod(3) = 1
754  nod(4) = 3
755  nod(5) = 12
756  nod(6) = 13
757  nod(7) = 9
758  nod(8) = 15
759  isuf = 2
760  endif
761  !
762  if( isuf.EQ.2 ) then
763  ic = 0
764  do ip = 1,8
765  term2(ip) = 0.0
766  do jp = 1,8
767  ic = ic + 1
768  term1(ic) = 0.0
769  do ig2=1,3
770  si=xg(ig2)
771  do ig1=1,3
772  ri=xg(ig1)
773  rp=1.0+ri
774  sp=1.0+si
775  rm=1.0-ri
776  sm=1.0-si
777  h(1)=0.25*rm*sm*(-1.0-ri-si)
778  h(2)=0.25*rp*sm*(-1.0+ri-si)
779  h(3)=0.25*rp*sp*(-1.0+ri+si)
780  h(4)=0.25*rm*sp*(-1.0-ri+si)
781  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
782  h(6)=0.5*(1.0-si*si)*(1.0+ri)
783  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
784  h(8)=0.5*(1.0-si*si)*(1.0-ri)
785  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
786  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
787  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
788  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
789  hr(5)=-ri*(1.0-si)
790  hr(6)= .5*(1.0-si*si)
791  hr(7)=-ri*(1.0+si)
792  hr(8)=-.5*(1.0-si*si)
793  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
794  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
795  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
796  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
797  hs(5)=-.5*(1.0-ri*ri)
798  hs(6)=-si*(1.0+ri)
799  hs(7)= .5*(1.0-ri*ri)
800  hs(8)=-si*(1.0-ri)
801  g1x=0.0
802  g1y=0.0
803  g1z=0.0
804  g2x=0.0
805  g2y=0.0
806  g2z=0.0
807  do i = 1,8
808  g1x=g1x+hr(i)*xx(nod(i))
809  g1y=g1y+hr(i)*yy(nod(i))
810  g1z=g1z+hr(i)*zz(nod(i))
811  g2x=g2x+hs(i)*xx(nod(i))
812  g2y=g2y+hs(i)*yy(nod(i))
813  g2z=g2z+hs(i)*zz(nod(i))
814  enddo
815  g3x=g1y*g2z-g1z*g2y
816  g3y=g1z*g2x-g1x*g2z
817  g3z=g1x*g2y-g1y*g2x
818  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
819  wg=xsum*wgt(ig1)*wgt(ig2)*hh
820  term1(ic)=term1(ic)-wg*h(ip)*h(jp)
821  if( ip.EQ.jp ) term2(ip)=term2(ip)-wg*h(jp)*sink
822  enddo
823  enddo
824  enddo
825  enddo
826  !
827  elseif( isuf.EQ.1 ) then
828  ic = 0
829  do ip = 1, 6
830  term2(ip) = 0.0
831  do jp = 1, 6
832  ic = ic + 1
833  term1(ic) = 0.0
834  ! INTEGRATION OVER SURFACE
835  do l2=1,3
836  xl2=xg(l2)
837  x2 =(xl2+1.0)*0.5
838  do l1=1,3
839  xl1=xg(l1)
840  x1=0.5*(1.0-x2)*(xl1+1.0)
841  ! INTERPOLATION FUNCTION
842  x3=1.0-x1-x2
843  h(1)= x1*(2.0*x1-1.)
844  h(2)= x2*(2.0*x2-1.)
845  h(3)= x3*(2.0*x3-1.)
846  h(4)= 4.0*x1*x2
847  h(5)= 4.0*x2*x3
848  h(6)= 4.0*x1*x3
849  ! DERIVATIVE OF INTERPOLATION FUNCTION
850  ! FOR L1-COORDINATE
851  hl1(1)=4.0*x1-1.0
852  hl1(2)= 0.0
853  hl1(3)= 0.0
854  hl1(4)= 4.0*x2
855  hl1(5)= 0.0
856  hl1(6)= 4.0*x3
857  ! FOR L2-COORDINATE
858  hl2(1)= 0.0
859  hl2(2)= 4.0*x2-1.0
860  hl2(3)= 0.0
861  hl2(4)= 4.0*x1
862  hl2(5)= 4.0*x3
863  hl2(6)= 0.0
864  ! FOR L3-COORDINATE
865  hl3(1)= 0.0
866  hl3(2)= 0.0
867  hl3(3)= 4.0*x3-1.0
868  hl3(4)= 0.0
869  hl3(5)= 4.0*x2
870  hl3(6)= 4.0*x1
871  !
872  g1x=0.0
873  g1y=0.0
874  g1z=0.0
875  g2x=0.0
876  g2y=0.0
877  g2z=0.0
878  do i=1,6
879  g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
880  g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
881  g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
882  g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
883  g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
884  g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
885  enddo
886  g3x=g1y*g2z-g1z*g2y
887  g3y=g1z*g2x-g1x*g2z
888  g3z=g1x*g2y-g1y*g2x
889  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
890  g3x=g3x/xsum
891  g3y=g3y/xsum
892  g3z=g3z/xsum
893  !JACOBI MATRIX
894  xj11=g1x
895  xj12=g1y
896  xj13=g1z
897  xj21=g2x
898  xj22=g2y
899  xj23=g2z
900  xj31=g3x
901  xj32=g3y
902  xj33=g3z
903  !DETERMINANT OF JACOBIAN
904  det=xj11*xj22*xj33 &
905  +xj12*xj23*xj31 &
906  +xj13*xj21*xj32 &
907  -xj13*xj22*xj31 &
908  -xj12*xj21*xj33 &
909  -xj11*xj23*xj32
910  wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25*hh
911  term1(ic) = term1(ic)-wg*h(ip)*h(jp)
912  if( ip.EQ.jp ) term2(ip) = term2(ip)-wg*h(jp)*sink
913  !
914  enddo
915  enddo
916  enddo
917  enddo
918  endif
919  !
920  return
921 
922  end subroutine heat_film_352
923  !----------------------------------------------------------------------*
924  subroutine heat_film_361(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
925  !----------------------------------------------------------------------*
926  !**
927  !** SET 361 FILM
928  !**
929  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
930  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
931  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
932  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
933  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
934  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
935  use hecmw
936  implicit none
937  ! I/F VARIABLES
938  integer(kind=kint) :: NN, LTYPE, MM
939  integer(kind=kint) :: NOD(MM)
940  real(kind=kreal) :: hh, sink
941  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
942  ! LOCAL VARIABLES
943  integer(kind=kint) :: I, IC, IP, JP, IG1, IG2
944  real(kind=kreal) :: ri, si, xsum
945  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
946  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
947  real(kind=kreal) :: xg(2), wgt(2)
948  !*************************
949  ! GAUSS INTEGRATION POINT
950  !*************************
951  data xg/-0.5773502691896,0.5773502691896/
952  data wgt/1.0,1.0/
953  !
954  if( ltype.EQ.1 ) then
955  nod(1) = 1
956  nod(2) = 2
957  nod(3) = 3
958  nod(4) = 4
959  else if( ltype.EQ.2 ) then
960  nod(1) = 8
961  nod(2) = 7
962  nod(3) = 6
963  nod(4) = 5
964  else if( ltype.EQ.3 ) then
965  nod(1) = 5
966  nod(2) = 6
967  nod(3) = 2
968  nod(4) = 1
969  else if( ltype.EQ.4 ) then
970  nod(1) = 6
971  nod(2) = 7
972  nod(3) = 3
973  nod(4) = 2
974  else if( ltype.EQ.5 ) then
975  nod(1) = 7
976  nod(2) = 8
977  nod(3) = 4
978  nod(4) = 3
979  else if( ltype.EQ.6 ) then
980  nod(1) = 8
981  nod(2) = 5
982  nod(3) = 1
983  nod(4) = 4
984  endif
985  !
986  ic = 0
987  do ip = 1, 4
988  term2(ip) = 0.0
989  do jp = 1, 4
990  ic = ic + 1
991  term1(ic) = 0.0
992  do ig2=1,2
993  si=xg(ig2)
994  do ig1=1,2
995  ri=xg(ig1)
996  h(1)=0.25*(1.0-ri)*(1.0-si)
997  h(2)=0.25*(1.0+ri)*(1.0-si)
998  h(3)=0.25*(1.0+ri)*(1.0+si)
999  h(4)=0.25*(1.0-ri)*(1.0+si)
1000  hr(1)=-.25*(1.0-si)
1001  hr(2)= .25*(1.0-si)
1002  hr(3)= .25*(1.0+si)
1003  hr(4)=-.25*(1.0+si)
1004  hs(1)=-.25*(1.0-ri)
1005  hs(2)=-.25*(1.0+ri)
1006  hs(3)= .25*(1.0+ri)
1007  hs(4)= .25*(1.0-ri)
1008  g1x=0.0
1009  g1y=0.0
1010  g1z=0.0
1011  g2x=0.0
1012  g2y=0.0
1013  g2z=0.0
1014  do i = 1, 4
1015  g1x=g1x+hr(i)*xx(nod(i))
1016  g1y=g1y+hr(i)*yy(nod(i))
1017  g1z=g1z+hr(i)*zz(nod(i))
1018  g2x=g2x+hs(i)*xx(nod(i))
1019  g2y=g2y+hs(i)*yy(nod(i))
1020  g2z=g2z+hs(i)*zz(nod(i))
1021  enddo
1022  g3x=g1y*g2z-g1z*g2y
1023  g3y=g1z*g2x-g1x*g2z
1024  g3z=g1x*g2y-g1y*g2x
1025  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1026  !
1027  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
1028  if( ip.EQ.jp ) then
1029  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
1030  endif
1031  !
1032  enddo
1033  enddo
1034  enddo
1035  enddo
1036 
1037  return
1038 
1039  end subroutine heat_film_361
1040  !----------------------------------------------------------------------*
1041  subroutine heat_film_362(NN,XX,YY,ZZ,LTYPE,HH,SINK,MM,TERM1,TERM2,NOD)
1042  !----------------------------------------------------------------------*
1043  !**
1044  !** SET 362 FILM
1045  !**
1046  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
1047  ! F2 LTYPE=2 : FILM IN NORMAL-DIRECTION FOR FACE-2
1048  ! F3 LTYPE=3 : FILM IN NORMAL-DIRECTION FOR FACE-3
1049  ! F4 LTYPE=4 : FILM IN NORMAL-DIRECTION FOR FACE-4
1050  ! F5 LTYPE=5 : FILM IN NORMAL-DIRECTION FOR FACE-5
1051  ! F6 LTYPE=6 : FILM IN NORMAL-DIRECTION FOR FACE-6
1052  use hecmw
1053  implicit none
1054  ! I/F VARIABLES
1055  integer(kind=kint) :: NN, LTYPE, MM
1056  integer(kind=kint) :: NOD(MM)
1057  real(kind=kreal) :: hh, sink
1058  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), term1(mm*mm), term2(mm)
1059  ! LOCAL VARIABLES
1060  integer(kind=kint) :: I, IC, IP, JP, IG1, IG2
1061  real(kind=kreal) :: ri, si, rp, sp, rm, sm
1062  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
1063  real(kind=kreal) :: xsum
1064  real(kind=kreal) :: h(8), hr(8), hs(8), ht(8)
1065  real(kind=kreal) :: xg(3), wgt(3)
1066  !*************************
1067  ! GAUSS INTEGRATION POINT
1068  !*************************
1069  xg(1) = -0.7745966692
1070  xg(2) = 0.0
1071  xg(3) = 0.7745966692
1072  wgt(1) = 0.5555555555
1073  wgt(2) = 0.8888888888
1074  wgt(3) = 0.5555555555
1075  !
1076  if( ltype.EQ.1 ) then
1077  nod(1) = 1
1078  nod(2) = 2
1079  nod(3) = 3
1080  nod(4) = 4
1081  nod(5) = 9
1082  nod(6) = 10
1083  nod(7) = 11
1084  nod(8) = 12
1085  else if( ltype.EQ.2 ) then
1086  nod(1) = 8
1087  nod(2) = 7
1088  nod(3) = 6
1089  nod(4) = 5
1090  nod(5) = 15
1091  nod(6) = 14
1092  nod(7) = 13
1093  nod(8) = 16
1094  else if( ltype.EQ.3 ) then
1095  nod(1) = 5
1096  nod(2) = 6
1097  nod(3) = 2
1098  nod(4) = 1
1099  nod(5) = 13
1100  nod(6) = 18
1101  nod(7) = 9
1102  nod(8) = 17
1103  else if( ltype.EQ.4 ) then
1104  nod(1) = 6
1105  nod(2) = 7
1106  nod(3) = 3
1107  nod(4) = 2
1108  nod(5) = 14
1109  nod(6) = 19
1110  nod(7) = 10
1111  nod(8) = 18
1112  else if( ltype.EQ.5 ) then
1113  nod(1) = 7
1114  nod(2) = 8
1115  nod(3) = 4
1116  nod(4) = 3
1117  nod(5) = 15
1118  nod(6) = 20
1119  nod(7) = 11
1120  nod(8) = 19
1121  else if( ltype.EQ.6 ) then
1122  nod(1) = 8
1123  nod(2) = 5
1124  nod(3) = 1
1125  nod(4) = 4
1126  nod(5) = 16
1127  nod(6) = 17
1128  nod(7) = 12
1129  nod(8) = 20
1130  endif
1131 
1132  ic = 0
1133  do ip = 1,8
1134  term2(ip) = 0.0
1135  do jp = 1,8
1136  ic = ic + 1
1137  term1(ic) = 0.0
1138  !
1139  do ig2=1,3
1140  si=xg(ig2)
1141  do ig1=1,3
1142  ri=xg(ig1)
1143 
1144  rp=1.0+ri
1145  sp=1.0+si
1146  rm=1.0-ri
1147  sm=1.0-si
1148  h(1)=0.25*rm*sm*(-1.0-ri-si)
1149  h(2)=0.25*rp*sm*(-1.0+ri-si)
1150  h(3)=0.25*rp*sp*(-1.0+ri+si)
1151  h(4)=0.25*rm*sp*(-1.0-ri+si)
1152  h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1153  h(6)=0.5*(1.0-si*si)*(1.0+ri)
1154  h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1155  h(8)=0.5*(1.0-si*si)*(1.0-ri)
1156  hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1157  hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1158  hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1159  hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1160  hr(5)=-ri*(1.0-si)
1161  hr(6)= .5*(1.0-si*si)
1162  hr(7)=-ri*(1.0+si)
1163  hr(8)=-.5*(1.0-si*si)
1164  hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1165  hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1166  hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1167  hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1168  hs(5)=-.5*(1.0-ri*ri)
1169  hs(6)=-si*(1.0+ri)
1170  hs(7)= .5*(1.0-ri*ri)
1171  hs(8)=-si*(1.0-ri)
1172 
1173  g1x=0.0
1174  g1y=0.0
1175  g1z=0.0
1176  g2x=0.0
1177  g2y=0.0
1178  g2z=0.0
1179 
1180  do i = 1,8
1181  g1x=g1x+hr(i)*xx(nod(i))
1182  g1y=g1y+hr(i)*yy(nod(i))
1183  g1z=g1z+hr(i)*zz(nod(i))
1184  g2x=g2x+hs(i)*xx(nod(i))
1185  g2y=g2y+hs(i)*yy(nod(i))
1186  g2z=g2z+hs(i)*zz(nod(i))
1187  enddo
1188 
1189  g3x=g1y*g2z-g1z*g2y
1190  g3y=g1z*g2x-g1x*g2z
1191  g3z=g1x*g2y-g1y*g2x
1192  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1193  !
1194  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
1195  if( ip.EQ.jp ) then
1196  term2(ip) = term2(ip)- xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
1197  endif
1198  !
1199  enddo
1200  enddo
1201  enddo
1202  enddo
1203 
1204  return
1205 
1206  end subroutine heat_film_362
1207  !----------------------------------------------------------------------*
1208  subroutine heat_film_731( NN,XX,YY,ZZ,HH,SINK,TERM1,TERM2 )
1209  !----------------------------------------------------------------------*
1210  !**
1211  !** SET 731 FILM
1212  !**
1213  ! F1 LTYPE=1 : SURFACE FILM
1214  use hecmw
1215  implicit none
1216  ! I/F VARIABLES
1217  integer(kind=kint) :: NN
1218  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), hh, sink, term1(nn*nn), term2(nn)
1219  ! LOCAL VARIABLES
1220  real(kind=kreal) :: ax, ay, az, bx, by, bz, aa
1221  integer(kind=kint) :: IC, IP, JP
1222  !
1223  ax = xx(2) - xx(1)
1224  ay = yy(2) - yy(1)
1225  az = zz(2) - zz(1)
1226  bx = xx(3) - xx(1)
1227  by = yy(3) - yy(1)
1228  bz = zz(3) - zz(1)
1229  aa = dsqrt( (ay*bz-az*by)**2+(az*bx-ax*bz)**2+(ax*by-ay*bx)**2 )/6.0
1230  !
1231  ic = 0
1232  do ip = 1, nn
1233  term2(ip) = - aa*hh*sink
1234  do jp = 1, nn
1235  ic = ic + 1
1236  term1(ic) = - aa*hh/3.0
1237  enddo
1238  enddo
1239 
1240  return
1241 
1242  end subroutine heat_film_731
1243  !----------------------------------------------------------------------*
1244  subroutine heat_film_741( NN,XX,YY,ZZ,HH,SINK,TERM1,TERM2 )
1245  !----------------------------------------------------------------------*
1246  !**
1247  !** SET 741 FILM
1248  !**
1249  ! F1 LTYPE=1 : FILM IN NORMAL-DIRECTION FOR FACE-1
1250  use hecmw
1251  implicit none
1252  ! I/F VARIABLES
1253  integer(kind=kint) :: NN
1254  real(kind=kreal) :: xx(nn), yy(nn), zz(nn), hh, sink, term1(nn*nn), term2(nn)
1255  ! LOCAL VARIABLES
1256  integer(kind=kint) :: I, IC, IP, JP, IG1, IG2
1257  real(kind=kreal) :: ri, si, xsum
1258  real(kind=kreal) :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, g3z
1259  real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
1260  real(kind=kreal) :: xg(2), wgt(2)
1261  !*************************
1262  ! GAUSS INTEGRATION POINT
1263  !*************************
1264  data xg/-0.5773502691896,0.5773502691896/
1265  data wgt/1.0,1.0/
1266  !
1267  ic = 0
1268  do ip = 1, nn
1269  term2(ip) = 0.0
1270  do jp = 1, nn
1271  ic = ic + 1
1272  term1(ic) = 0.0
1273  !
1274  do ig2=1,2
1275  si=xg(ig2)
1276  do ig1=1,2
1277  ri=xg(ig1)
1278  !
1279  h(1)=0.25*(1.0-ri)*(1.0-si)
1280  h(2)=0.25*(1.0+ri)*(1.0-si)
1281  h(3)=0.25*(1.0+ri)*(1.0+si)
1282  h(4)=0.25*(1.0-ri)*(1.0+si)
1283  hr(1)=-.25*(1.0-si)
1284  hr(2)= .25*(1.0-si)
1285  hr(3)= .25*(1.0+si)
1286  hr(4)=-.25*(1.0+si)
1287  hs(1)=-.25*(1.0-ri)
1288  hs(2)=-.25*(1.0+ri)
1289  hs(3)= .25*(1.0+ri)
1290  hs(4)= .25*(1.0-ri)
1291  !
1292  g1x=0.0
1293  g1y=0.0
1294  g1z=0.0
1295  g2x=0.0
1296  g2y=0.0
1297  g2z=0.0
1298  do i = 1, nn
1299  g1x=g1x+hr(i)*xx(i)
1300  g1y=g1y+hr(i)*yy(i)
1301  g1z=g1z+hr(i)*zz(i)
1302  g2x=g2x+hs(i)*xx(i)
1303  g2y=g2y+hs(i)*yy(i)
1304  g2z=g2z+hs(i)*zz(i)
1305  enddo
1306  g3x=g1y*g2z-g1z*g2y
1307  g3y=g1z*g2x-g1x*g2z
1308  g3z=g1x*g2y-g1y*g2x
1309  xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1310  !
1311  term1(ic)=term1(ic)-xsum*wgt(ig1)*wgt(ig2)*h(ip)*h(jp)*hh
1312  if( ip.EQ.jp ) then
1313  term2(ip)=term2(ip)-xsum*wgt(ig1)*wgt(ig2)*h(jp)*hh*sink
1314  endif
1315  !
1316  enddo
1317  enddo
1318  enddo
1319  enddo
1320 
1321  return
1322 
1323  end subroutine heat_film_741
1324 end module m_heat_lib_film
m_heat_lib_film::heat_film_741
subroutine heat_film_741(NN, XX, YY, ZZ, HH, SINK, TERM1, TERM2)
Definition: heat_LIB_FILM.f90:1245
m_heat_lib_film::heat_film_361
subroutine heat_film_361(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:925
m_heat_lib_film::heat_film_731
subroutine heat_film_731(NN, XX, YY, ZZ, HH, SINK, TERM1, TERM2)
Definition: heat_LIB_FILM.f90:1209
m_heat_lib_film::heat_film_351
subroutine heat_film_351(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:537
m_heat_lib_film::heat_film_341
subroutine heat_film_341(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:317
m_heat_lib_film::heat_film_231
subroutine heat_film_231(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:24
m_heat_lib_film
This module provides subroutines to generate heat transfer boundary.
Definition: heat_LIB_FILM.f90:7
m_heat_lib_film::heat_film_352
subroutine heat_film_352(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:677
hecmw
Definition: hecmw.f90:6
m_heat_lib_film::heat_film_241
subroutine heat_film_241(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:166
m_heat_lib_film::heat_film_362
subroutine heat_film_362(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:1042
m_heat_lib_film::heat_film_242
subroutine heat_film_242(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:235
m_heat_lib_film::heat_film_232
subroutine heat_film_232(NN, XX, YY, ZZ, THICK, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:89
m_heat_lib_film::heat_film_342
subroutine heat_film_342(NN, XX, YY, ZZ, LTYPE, HH, SINK, MM, TERM1, TERM2, NOD)
Definition: heat_LIB_FILM.f90:378