FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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!-------------------------------------------------------------------------------
8contains
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
1303end module m_heat_lib_radiate
Definition: hecmw.f90:6
This module provides subroutines to generate heat radiate boundary.
subroutine heat_radiate_342(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_242(nn, xx, yy, zz, thick, temp, ltype, rr, sink, tzero, mm, term1, term2, nod)
subroutine heat_radiate_352(nn, xx, yy, zz, 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_731(nn, xx, yy, zz, temp, ltype, rr, sink, tzero, term1, term2)
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_241(nn, xx, yy, zz, thick, temp, ltype, rr, sink, tzero, mm, term1, term2, nod)
subroutine heat_radiate_362(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)