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