1!/*****************************************************************************/
2! *
3! *  Elmer, A Finite Element Software for Multiphysical Problems
4! *
5! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6! *
7! * This library is free software; you can redistribute it and/or
8! * modify it under the terms of the GNU Lesser General Public
9! * License as published by the Free Software Foundation; either
10! * version 2.1 of the License, or (at your option) any later version.
11! *
12! * This library is distributed in the hope that it will be useful,
13! * but WITHOUT ANY WARRANTY; without even the implied warranty of
14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15! * Lesser General Public License for more details.
16! *
17! * You should have received a copy of the GNU Lesser General Public
18! * License along with this library (in file ../LGPL-2.1); if not, write
19! * to the Free Software Foundation, Inc., 51 Franklin Street,
20! * Fifth Floor, Boston, MA  02110-1301  USA
21! *
22! *****************************************************************************/
23!
24!/******************************************************************************
25! *
26! *  Authors: Juha Ruokolainen
27! *  Email:   Juha.Ruokolainen@csc.fi
28! *  Web:     http://www.csc.fi/elmer
29! *  Address: CSC - IT Center for Science Ltd.
30! *           Keilaranta 14
31! *           02101 Espoo, Finland
32! *
33! *  Original Date: 02 Jun 1997
34! *
35! *****************************************************************************/
36
37!> \ingroup ElmerLib
38!> \{
39
40!-----------------------------------------------------------------------------
41!>  Utility routines for radiation computation
42!-----------------------------------------------------------------------------
43
44MODULE Radiation
45
46   USE ElementUtils
47   USE CoordinateSystems
48   USE DefUtils
49
50   IMPLICIT NONE
51
52CONTAINS
53
54!------------------------------------------------------------------------------
55   FUNCTION ComputeRadiationLoad( Model, Mesh, Element, Temperature, &
56                 Reorder, Emissivity, AngleFraction, Areas, Emiss ) RESULT(T)
57!------------------------------------------------------------------------------
58     TYPE(Mesh_t), POINTER :: Mesh
59     TYPE(Model_t) :: Model
60     TYPE(Element_t)  :: Element
61     INTEGER :: Reorder(:)
62     REAL(KIND=dp), OPTIONAL :: AngleFraction, Areas(:), Emiss(:)
63     REAL(KIND=dp) :: T
64     REAL(KIND=dp) :: Temperature(:), Emissivity
65
66     REAL(KIND=dp) :: Asum
67     TYPE(Element_t),POINTER  :: RadElement
68     INTEGER :: i,j,n, bindex,nf
69     REAL(KIND=dp), POINTER :: Vals(:)
70     INTEGER, POINTER :: Cols(:)
71     REAL(KIND=dp) :: A1,A2,Emissivity1
72     LOGICAL :: Found
73!------------------------------------------------------------------------------
74
75     IF( .NOT. ASSOCIATED( Element % BoundaryInfo % GebhardtFactors ) ) THEN
76       CALL Fatal('ComputeRadiationLoad','Gebhardt factors not calculated for boundary!')
77     END IF
78
79     nf = Element % BoundaryInfo % GebhardtFactors % NumberOfFactors
80
81     IF(PRESENT(Areas) .AND. PRESENT(Emiss) ) THEN
82
83       bindex = Element % ElementIndex - Mesh % NumberOfBulkElements
84       A1  = Emiss(bIndex)**2
85
86       Cols => Element % BoundaryInfo % GebhardtFactors % Elements
87       Vals => Element % BoundaryInfo % GebhardtFactors % Factors
88
89       T = 0._dp
90       Asum = 0._dp
91       DO i=1,nf
92         RadElement => Mesh % Elements(Cols(i))
93         n = RadElement % TYPE % NumberOfNodes
94
95         bindex = Cols(i) - Mesh % NumberOfBulkElements
96         A2 = Emiss(bindex)
97
98         T=T+A2*Vals(i)*SUM(Temperature(Reorder(RadElement % NodeIndexes))/n)**4
99         Asum = Asum + A2*Vals(i)
100       END DO
101     ELSE
102       A1 = Emissivity**2
103
104       Cols => Element % BoundaryInfo % GebhardtFactors % Elements
105       Vals => Element % BoundaryInfo % GebhardtFactors % Factors
106
107       T = 0.0_dp
108       Asum = 0.0_dp
109       DO i=1,nf
110         RadElement => Mesh % Elements(Cols(i))
111         n = RadElement % TYPE % NumberOfNodes
112
113         Emissivity1 = SUM(ListGetReal( Model % BCs(RadElement % &
114              BoundaryInfo % Constraint) % Values, 'Emissivity', &
115              n, RAdElement % NodeIndexes, Found) ) / n
116         IF(.NOT. Found) THEN
117            Emissivity1 = SUM(GetParentMatProp('Emissivity',RadElement)) / n
118         END IF
119
120         A2 = Emissivity1
121         T = T + A2 * Vals(i) * &
122           SUM(Temperature(Reorder(RadElement % NodeIndexes))/n)**4
123
124         Asum = Asum + A2 * Vals(i)
125       END DO
126     END IF
127
128     T = (T/A1)**(1._dp/4._dp)
129
130     IF(PRESENT(AngleFraction)) AngleFraction = Asum
131!------------------------------------------------------------------------------
132   END FUNCTION ComputeRadiationLoad
133!------------------------------------------------------------------------------
134
135
136
137!------------------------------------------------------------------------------
138   FUNCTION ComputeRadiationCoeff( Model,Mesh,Element,k ) RESULT(T)
139!------------------------------------------------------------------------------
140
141     TYPE(Mesh_t), POINTER :: Mesh
142     TYPE(Model_t)  :: Model
143     TYPE(Element_t) :: Element
144     INTEGER :: k
145!------------------------------------------------------------------------------
146
147     REAL(KIND=dp) :: T
148
149     TYPE(Element_t),POINTER  :: CurrentElement
150     INTEGER :: i,j,n
151     LOGICAL :: Found
152
153     REAL(KIND=dp) :: Area,Emissivity
154!------------------------------------------------------------------------------
155
156     CurrentElement => Model % Elements( &
157             Element % BoundaryInfo % GebhardtFactors % Elements(k) )
158     n = CurrentElement % TYPE % NumberOfNodes
159
160     Emissivity = SUM(ListGetReal(Model % BCs(CurrentElement % &
161        BoundaryInfo % Constraint) % Values, 'Emissivity', &
162        n, CurrentElement % NodeIndexes, Found)) / n
163     IF(.NOT. Found) THEN
164        Emissivity = SUM(GetParentMatProp('Emissivity',CurrentElement)) / n
165     END IF
166
167     Area = Emissivity * ElementArea( Mesh,CurrentElement, n)
168
169     T =  ABS(Element % BoundaryInfo % GebhardtFactors % Factors(k)) * Area
170
171   END FUNCTION ComputeRadiationCoeff
172!------------------------------------------------------------------------------
173
174!------------------------------------------------------------------------------
175END MODULE Radiation
176!------------------------------------------------------------------------------
177
178!> \}
179