CoCalc -- USF_GetFrictionHeating.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/UserFunctions/USF_GetFrictionHeating.F90
³¹⁹⁶ views
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 program is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU General Public License
9
! *  as published by the Free Software Foundation; either version 2
10
! *  of the License, or (at your option) any later version.
11
! * 
12
! *  This program 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
15
! *  GNU General Public License for more details.
16
! *
17
! *  You should have received a copy of the GNU General Public License
18
! *  along with this program (in file fem/GPL-2); if not, write to the 
19
! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
20
! *  Boston, MA 02110-1301, USA.
21
! *
22
! *****************************************************************************/
23
!
24
!/******************************************************************************
25
! *
26
! *  Module containing a functions for friction heat
27
! *
28
! ******************************************************************************
29
! *
30
! *  Authors: Thomas Zwinger, Juha Ruokolainen,  Joe Todd
31
! *  Email:   [email protected]
32
! *  Web:     http://www.csc.fi/elmer
33
! *  Address: CSC - IT Center for Science Ltd.
34
! *           Keilaranta 14
35
! *           02101 Espoo, Finland 
36
! *
37
! *  Original Date: 08 Jun 1997
38
! *  Current date: 20 July 2012 (Martina)
39
! *
40
! *****************************************************************************/
41

42

43
! contains function for adding friction heat to geothermal heat flux at the basis and exporting it
44
! ATTENTION: consistent units for capacity and conductivy are needed.
45
! ATTENTION: The keyword "Friction Heat" does make reference to Strainheating, not Friction heat!
46
! usage:
47

48

49

50

51
FUNCTION getFrictionHeat(  Model, Node, DummyInput)RESULT(frictionheat)
52
  
53
  USE DefUtils
54

55
  IMPLICIT NONE
56
  
57
  !------------------------------------------------------------------------------
58
  TYPE(Model_t) :: Model
59
  INTEGER :: Node
60
  REAL(KIND=dp) :: DummyInput, frictionHeat
61
  !----------------------------------------------------------------------------
62
  
63
  INTEGER :: DIM, i, j, Ind(3,3)
64
  REAL(KIND=dp), POINTER :: FlowValues(:),NormalValues(:),StressValues(:)
65
  REAL(KIND=dp) :: normal(3), velo(3), un, ut, Sig(3,3), Sn(3), snn, snt
66
  INTEGER, POINTER :: FlowPerm(:),StressPerm(:), NormalPerm(:)
67
  LOGICAL :: FirstTime=.TRUE.,UnFoundFatal=.TRUE.
68
  TYPE(Variable_t), POINTER :: FlowVar,StressVariable, NormalVar
69
  CHARACTER(LEN=MAX_NAME_LEN) :: FunctionName
70
  
71
  SAVE FirstTime, DIM, FunctionName,Ind
72
  
73
  IF (FirstTime) THEN
74
     WRITE(FunctionName, '(A)') 'getFrictionHeat'
75
     DIM = CoordinateSystemDimension()
76
     FirstTime = .FALSE.
77
     DO i=1, 3
78
        Ind(i,i) = i
79
     END DO
80
     Ind(1,2) = 4
81
     Ind(2,1) = 4
82
     Ind(2,3) = 5
83
     Ind(3,2) = 5
84
     Ind(3,1) = 6
85
     Ind(1,3) = 6
86
  END IF
87
  
88
  ! Get the variable velocity
89
  !---------------------------
90
  FlowVar => VariableGet( Model % Variables, 'Flow Solution',UnFoundFatal=UnFoundFatal)
91
  FlowPerm    => FlowVar % Perm
92
  FlowValues  => FlowVar % Values
93
  
94
  ! Get the stress variable
95
  !------------------------
96
  StressVariable => VariableGet( Model % Variables, 'Stress',UnFoundFatal=UnFoundFatal)
97
  StressPerm    => StressVariable % Perm
98
  StressValues  => StressVariable % Values
99
  
100
  ! Get the variable for normal vector
101
  NormalVar =>  VariableGet(Model % Variables,'Normal Vector',UnFoundFatal=UnFoundFatal)
102
  NormalPerm => NormalVar % Perm
103
  NormalValues => NormalVar % Values
104
  
105
  DO i=1, DIM
106
     normal(i) = NormalValues(DIM*(NormalPerm(Node)-1) + i)      
107
     velo(i) = FlowValues( (DIM+1)*(FlowPerm(Node)-1) + i )
108
  END DO
109
  
110
  !Tangential velocity
111
  un = SUM(velo(1:DIM)*(normal(1:DIM))) 
112
  ut = SQRT(SUM( (velo(1:DIM))**2.0_dp ) - (un**2.0_dp) )
113

114
  !Tangential Stress
115
  DO i=1, DIM
116
     DO j= 1, DIM
117
        Sig(i,j) =  &
118
             StressValues( 2*DIM *(StressPerm(Node)-1) + Ind(i,j) )
119
     END DO
120
  END DO
121
  DO i=1, DIM
122
     Sn(i) = SUM(Sig(i,1:DIM)*normal(1:DIM)) 
123
  END DO
124
  Snn = SUM( Sn(1:DIM) * normal(1:DIM) )
125
  Snt = SQRT( SUM(Sn(1:DIM)**2.0_dp) - (Snn**2.0_dp))
126

127
  frictionHeat =ut*Snt
128
END FUNCTION getFrictionHeat
129

130
!/******************************************************************************
131
! *
132
! *  Module containing a functions for friction heat based
133
! *       on residuals from (Navier-)Stokes solver
134
! *  This function should be preferably used to compute the heat production
135
! *  at the base, as it utilizes the natural way to couple the flow
136
! *  and temperature solution
137
! *
138
! ******************************************************************************
139
! *
140
! *  Authors: Thomas Zwinger, Juha Ruokolainen, Martina Schäfer, Olivier Gagliardini
141
! *  Email:   [email protected]
142
! *  Web:     http://www.csc.fi/elmer
143
! *  Address: CSC - IT Center for Science Ltd.
144
! *           Keilaranta 14
145
! *           02101 Espoo, Finland 
146
! *
147
! *  Original Date: 08 Jun 1997
148
! *  Current date:  28 August 2014 (Martina/Thomas)
149
! *
150
! *****************************************************************************/
151

152
FUNCTION getFrictionLoads(  Model, Node, DummyInput )RESULT(frictionLoad)
153
  
154
  USE DefUtils
155

156
  IMPLICIT NONE
157
  
158
  !------------------------------------------------------------------------------
159
  TYPE(Model_t) :: Model
160
  INTEGER :: Node
161
  REAL(KIND=dp) :: DummyInput, frictionLoad
162
  !----------------------------------------------------------------------------
163

164
  INTEGER :: DIM, i, other_body_id
165
  REAL(KIND=dp), POINTER :: FlowValues(:),FlowLoadValues(:),NormalValues(:),MaskValues(:)
166
  REAL(KIND=dp) :: normal(3), velo(3), normalvelocity, flowload(3), tangvelocity(3)
167
  INTEGER, POINTER :: FlowPerm(:),FlowLoadPerm(:), NormalPerm(:), MaskPerm(:)
168
  LOGICAL :: FirstTime=.TRUE., GotIt,UnFoundFatal,UseMask = .FALSE., Warned=.FALSE.
169
  TYPE(Variable_t), POINTER :: FlowVar,FlowLoadVar, NormalVar, MaskVar
170
  TYPE(ValueList_t), POINTER :: Equation
171
  CHARACTER(LEN=MAX_NAME_LEN) :: FlowSolutionName, FlowLoadsName, MaskName
172
  CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: FunctionName='USF_GetFrictionHeating(getFrictionLoads)'
173
  TYPE(Element_t), POINTER ::  BoundaryElement, ParentElement
174
  
175
  SAVE FirstTime,  DIM, UseMask
176

177
  
178
  IF (FirstTime) THEN
179
    !WRITE(FunctionName,'(A)') 'USF_GetFrictionHeating(getFrictionLoads)'    
180

181
    DIM = CoordinateSystemDimension()
182
  END IF
183
   
184
  ! Get variable names from Equation section
185
  !-----------------------------------------
186
  BoundaryElement => Model % CurrentElement
187
  other_body_id = BoundaryElement % BoundaryInfo % outbody
188
  IF (other_body_id < 1) THEN ! only one body in calculation
189
    ParentElement => BoundaryElement % BoundaryInfo % Right
190
    IF ( .NOT. ASSOCIATED(ParentElement) ) ParentElement => BoundaryElement % BoundaryInfo % Left
191
  ELSE ! we are dealing with a body-body boundary and assume that the normal is pointing outwards
192
    ParentElement => BoundaryElement % BoundaryInfo % Right
193
    IF (ParentElement % BodyId == other_body_id) ParentElement => BoundaryElement % BoundaryInfo % Left
194
  END IF
195
  Equation => GetEquation(Element=ParentElement,Found=GotIt)
196
  IF (.NOT.ASSOCIATED(Equation) .OR. .NOT.GotIt) THEN
197
    IF (FirstTime) THEN
198
      WRITE (Message,'(A,I3)') 'No "Equation" found. Using default values for variables'
199
      CALL WARN(FunctionName,Message)    
200
      Warned = .TRUE.
201
    END IF
202
    WRITE(FlowSolutionName,'(A)') 'Flow Solution'
203
    WRITE(FlowLoadsName,'(A)') TRIM(FlowSolutionName)//' Loads'
204
    UseMask = .FALSE.
205
  ELSE
206
    FlowSolutionName = GetString( Equation , 'Flow Solution Name', GotIt )
207
    IF (.NOT.GotIt) THEN
208
      WRITE(FlowSolutionName,'(A)') 'Flow Solution'
209
      IF (FirstTime) THEN
210
        WRITE(Message,'(A,A)') 'Using default name for flow solution: ', &
211
             FlowSolutionName
212
        CALL WARN(FunctionName,Message)
213
        Warned = .TRUE.
214
      END IF
215
    END IF
216
    FlowLoadsName = GetString( Equation , 'Flow Loads Name', GotIt )  
217
    IF (.NOT. GotIt) THEN
218
      WRITE(FlowLoadsName,'(A)') TRIM(FlowSolutionName)//' Loads'
219
      IF (FirstTime) THEN
220
        WRITE(Message,'(A,A)') 'Using default name for flow solution loads: ', &
221
             FlowLoadsName
222
        CALL WARN(FunctionName,Message)
223
        Warned = .TRUE.
224
      END IF
225
    END IF
226
    MaskName = GetString( Equation , 'Friction Load Mask', UseMask )    
227
    IF (UseMask) THEN
228
      WRITE (Message, '(A,A)') '>Friction Load Mask< found and set to ', &
229
           TRIM(MaskName)
230
      CALL INFO(FunctionName, Message, Level=20)
231
    END IF
232
  END IF
233
  IF (Warned .AND. FirstTime) &
234
       CALL WARN(FunctionName,"All Warnings will be further omitted")
235
  FirstTime = .FALSE.
236
  
237
  ! Get the variable velocity
238
  !---------------------------
239
  FlowVar => VariableGet( Model % Variables, TRIM(FlowSolutionName),UnFoundFatal=UnFoundFatal)
240
  FlowPerm    => FlowVar % Perm
241
  FlowValues  => FlowVar % Values
242

243
  
244
  ! Get the Stokes loads
245
  !---------------------------
246
  FlowLoadVar => VariableGet( Model % Variables, TRIM(FlowLoadsName),UnFoundFatal=UnFoundFatal)
247
  FlowLoadPerm    => FlowLoadVar % Perm
248
  FlowLoadValues  => FlowLoadVar % Values
249
  
250

251
  ! Get the variable for normal vector
252
  !-----------------------------------
253
  NormalVar =>  VariableGet(Model % Variables,'Normal Vector',UnFoundFatal=UnFoundFatal)
254
  NormalPerm => NormalVar % Perm
255
  NormalValues => NormalVar % Values
256

257
  IF (UseMask) THEN
258
    MaskVar => VariableGet( Model % Variables, TRIM(MaskName),UnFoundFatal=UnFoundFatal)
259
    MaskPerm    => MaskVar % Perm
260
    MaskValues  => MaskVar % Values
261
  END IF
262

263
  IF (UseMask) THEN
264
    IF ( MaskValues(MaskPerm(Node)).LE.0.0_dp ) THEN
265
      frictionLoad = 0.0_dp
266
    END IF
267
  ELSE
268
    DO i=1, DIM
269
      normal(i) = NormalValues(DIM*(NormalPerm(Node)-1) + i)      
270
      velo(i) = FlowValues( (DIM+1)*(FlowPerm(Node)-1) + i )
271
      flowload(i) = FlowLoadValues( (DIM+1)*(FlowLoadPerm(Node)-1) + i )
272
    END DO
273
 
274
    normalvelocity    = SUM( velo(1:DIM) * normal(1:DIM) )
275

276
    DO i=1, DIM
277
      tangvelocity(i) = velo(i) - normalvelocity * normal(i)
278
    END DO
279
    
280
    frictionLoad = &
281
         MAX( (-1.0_dp * (SUM(tangvelocity(1:DIM) * flowLoad(1:DIM)))), 0.0_dp)
282
  END IF
283

284
END FUNCTION getFrictionLoads
285

286
Product

Resources

Company
