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