1!/*****************************************************************************/ 2! * Elmer/Ice, a glaciological add-on to Elmer 3! * http://elmerice.elmerfem.org 4! * 5! * This program is free software; you can redistribute it and/or 6! * modify it under the terms of the GNU General Public License 7! * as published by the Free Software Foundation; either version 2 8! * of the License, or (at your option) any later version. 9! * 10! * This program is distributed in the hope that it will be useful, 11! * but WITHOUT ANY WARRANTY; without even the implied warranty of 12! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! * GNU General Public License for more details. 14! * 15! * You should have received a copy of the GNU General Public License 16! * along with this program (in file fem/GPL-2); if not, write to the 17! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 18! * Boston, MA 02110-1301, USA. 19! * 20! *****************************************************************************/ 21! 22!/****************************************************************************** 23! * 24! * Module containing a solver for enthalpy equation in glaciology 25! * 26! ****************************************************************************** 27! * 28! * Authors: Adrien Gilbert, 29! * 30! * Original Date: Jun 2014 31! * 32! *****************************************************************************/ 33 34!------------------------------------------------------------------------------ 35 SUBROUTINE EnthalpySolver( Model,Solver,Timestep,TransientSimulation ) 36!------------------------------------------------------------------------------ 37 USE DiffuseConvective 38 USE DiffuseConvectiveGeneral 39 USE Differentials 40 USE Radiation 41 USE MaterialModels 42 USE Adaptive 43 USE DefUtils 44 45!------------------------------------------------------------------------------ 46 IMPLICIT NONE 47!------------------------------------------------------------------------------ 48 INTEGER, PARAMETER :: PHASE_SPATIAL_1 = 1 49 INTEGER, PARAMETER :: PHASE_SPATIAL_2 = 2 50 INTEGER, PARAMETER :: PHASE_TEMPORAL = 3 51 52 TYPE(Model_t) :: Model 53 TYPE(Solver_t), TARGET :: Solver 54 55 LOGICAL :: TransientSimulation 56 REAL(KIND=dp) :: Timestep 57!------------------------------------------------------------------------------ 58! Local variables 59!------------------------------------------------------------------------------ 60 TYPE(Matrix_t), POINTER :: StiffMatrix 61 62 INTEGER :: i,j,k,l,m,n,nd,t,tt,iter,k1,k2,body_id,eq_id,istat,LocalNodes,bf_id 63 64 TYPE(Nodes_t) :: ElementNodes 65 TYPE(Element_t),POINTER :: Element,Parent,RadiationElement 66 67 REAL(KIND=dp) :: RelativeChange, & 68 Norm,PrevNorm,Text,S,C,C1,Emissivity,StefanBoltzmann, & 69 ReferencePressure=0.0d0, SpecificHeatRatio 70 71 CHARACTER(LEN=MAX_NAME_LEN) :: RadiationFlag,ConvectionFlag 72 73 INTEGER :: PhaseChangeModel 74 CHARACTER(LEN=MAX_NAME_LEN) :: PhaseModel, StabilizeFlag, VarName 75 76 INTEGER, POINTER :: NodeIndexes(:) 77 LOGICAL :: Stabilize = .TRUE., Bubbles = .TRUE., UseBubbles,NewtonLinearization = .FALSE., & 78 Found, GotIt, HeatFluxBC, HeatGapBC, GotMeltPoint, IsRadiation, InfBC, UnFoundFatal=.TRUE. 79! Which compressibility model is used 80 CHARACTER(LEN=MAX_NAME_LEN) :: CompressibilityFlag, ConvectionField 81 INTEGER :: CompressibilityModel 82 83 LOGICAL :: AllocationsDone = .FALSE.,PhaseSpatial=.FALSE., & 84 PhaseChange=.FALSE., CheckLatentHeatRelease=.FALSE., FirstTime, & 85 SmartHeaterControl, IntegralHeaterControl, HeaterControlLocal, SmartTolReached=.FALSE., & 86 TransientHeaterControl, SmartHeaterAverage, ConstantBulk, SaveBulk, & 87 TransientAssembly 88 LOGICAL, POINTER :: SmartHeaters(:), IntegralHeaters(:) 89 90 TYPE(Variable_t), POINTER :: TempSol,FlowSol,HeatSol,CurrentSol, MeshSol, DensitySol 91 TYPE(ValueList_t), POINTER :: Equation,Material,SolverParams,BodyForce,BC,Constants 92 93 INTEGER, POINTER :: EnthalpyPerm(:),FlowPerm(:),CurrentPerm(:),MeshPerm(:) 94 95 INTEGER :: NSDOFs,NewtonIter,NonlinearIter,MDOFs, & 96 SmartHeaterBC, SmartHeaterNode, DoneTime=0 97 REAL(KIND=dp) :: NonlinearTol,NewtonTol,SmartTol,Relax, & 98 SaveRelax,dt,dt0,CumulativeTime, VisibleFraction, PowerScaling=1.0, PrevPowerScaling=1.0, & 99 PowerRelax, PowerTimeScale, PowerSensitivity, xave, yave, Normal(3), & 100 dist, mindist, ControlPoint(3) 101 102 REAL(KIND=dp), POINTER :: Enthalpy_h(:),PrevEnthalpy_h(:),FlowSolution(:), & 103 ElectricCurrent(:), PhaseChangeIntervals(:,:),ForceVector(:), & 104 PrevSolution(:), HC(:), Hwrk(:,:,:),MeshVelocity(:), XX(:), YY(:),ForceHeater(:),& 105 RealWork(:,:) 106 107 REAL(KIND=dp), ALLOCATABLE :: vals(:) 108 REAL(KIND=dp) :: Jx,Jy,Jz,JAbs, Power, MeltPoint, IntHeatSource 109 110 REAL(KIND=dp) :: Tref,Tm,P,hm,hi,L_heat,A_cap,B_cap,Ptriple,Psurf,beta 111 112 CHARACTER(LEN=MAX_NAME_LEN) :: PressureName,WaterName,PhaseEnthName,TempName 113 114 TYPE(Variable_t), POINTER :: PressureVariable,PhaseChangeEnthVar,WaterVar,TemphomoVar 115 REAL(KIND=dp), POINTER :: PressureValues(:),PhaseChangeEnthValues(:) 116 INTEGER, POINTER :: PressurePerm(:),PhaseChangeEnthPerm(:) 117 118 INTEGER, ALLOCATABLE, SAVE :: Indexes(:), SaveIndexes(:) 119 120 REAL(KIND=dp), ALLOCATABLE :: MASS(:,:), & 121 STIFF(:,:), LOAD(:), HeatConductivity(:,:,:), & 122 FORCE(:), U(:), V(:), W(:), MU(:,:),TimeForce(:), & 123 Density(:), LatentHeat(:), HeatTransferCoeff(:), & 124 HeatCapacity(:),WaterDiffusivity(:), Enthalpy(:), Viscosity(:), LocalEnthalpy_h(:), & 125 NodalEmissivity(:), ElectricConductivity(:), Permeability(:), Work(:), C0(:), & 126 Pressure(:), dPressuredt(:), GasConstant(:),AText(:), HeaterArea(:), & 127 HeaterTarget(:), HeaterScaling(:), HeaterDensity(:), HeaterSource(:), & 128 HeatExpansionCoeff(:), ReferenceEnthalpy_h(:), PressureCoeff(:), & 129 PhaseVelocity(:,:), HeatConductivityIso(:), & 130 PerfusionRate(:), PerfusionDensity(:), PerfusionHeatCapacity(:), PerfusionRefEnthalpy_h(:) 131 132 SAVE U, V, W, MU, MASS, STIFF, LOAD, PressureCoeff, & 133 FORCE, ElementNodes, HeatConductivity, HeatCapacity,WaterDiffusivity, HeatTransferCoeff, & 134 Enthalpy, Density, LatentHeat, PhaseVelocity, AllocationsDone, Viscosity, TimeForce, & 135 LocalNodes, LocalEnthalpy_h, Work, ElectricConductivity, & 136 NodalEmissivity, Permeability, C0, dPressuredt, Pressure, & 137 GasConstant,AText,Hwrk, XX, YY, ForceHeater, Power, HeaterArea, HeaterTarget, & 138 HeaterScaling, HeaterDensity, HeaterSource, SmartHeaters, IntegralHeaters, SmartTolReached, & 139 ReferenceEnthalpy_h, HeatExpansionCoeff, PrevPowerScaling, PowerScaling, & 140 MeltPoint, DoneTime, SmartHeaterNode, SmartHeaterBC, SmartHeaterAverage, & 141 HeatConductivityIso, & 142 PerfusionRate, PerfusionDensity, PerfusionHeatCapacity, PerfusionRefEnthalpy_h 143 144 145 INTERFACE 146 FUNCTION HeatBoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm ) RESULT(Indicator) 147 USE Types 148 TYPE(Element_t), POINTER :: Edge 149 TYPE(Model_t) :: Model 150 TYPE(Mesh_t), POINTER :: Mesh 151 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 152 INTEGER :: Perm(:) 153 END FUNCTION HeatBoundaryResidual 154 155 FUNCTION HeatEdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT(Indicator) 156 USE Types 157 TYPE(Element_t), POINTER :: Edge 158 TYPE(Model_t) :: Model 159 TYPE(Mesh_t), POINTER :: Mesh 160 REAL(KIND=dp) :: Quant(:), Indicator(2) 161 INTEGER :: Perm(:) 162 END FUNCTION HeatEdgeResidual 163 164 FUNCTION HeatInsideResidual( Model,Element,Mesh,Quant,Perm, Fnorm ) RESULT(Indicator) 165 USE Types 166 TYPE(Element_t), POINTER :: Element 167 TYPE(Model_t) :: Model 168 TYPE(Mesh_t), POINTER :: Mesh 169 REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm 170 INTEGER :: Perm(:) 171 END FUNCTION HeatInsideResidual 172 END INTERFACE 173 174#ifdef USE_ISO_C_BINDINGS 175 REAL(KIND=dp) :: at,at0,totat,st,totst,t1 176#else 177 REAL(KIND=dp) :: at,at0,totat,st,totst,t1,CPUTime,RealTime 178#endif 179 180!------------------------------------------------------------------------------ 181! The View and Gebhardt factors may change. If this is necessary, this is 182! done within this subroutine. The routine is called in the 183! start as it may affect the matrix topology. 184! Newton lineariarization option is needed only when there is radiation. 185!------------------------------------------------------------------------------ 186 IsRadiation = ListCheckPresentAnyBC( Model,'Radiation') 187 IF( IsRadiation ) THEN 188 CALL RadiationFactors( Solver, .FALSE.) 189 END IF 190 191!------------------------------------------------------------------------------ 192! Get variables needed for solution 193!------------------------------------------------------------------------------ 194 195 IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN 196 197 StiffMatrix => GetMatrix() 198 ForceVector => Solver % Matrix % RHS 199 200 TempSol => Solver % Variable 201 EnthalpyPerm => TempSol % Perm 202 Enthalpy_h => TempSol % Values 203 VarName = GetVarName( TempSol ) 204 205 LocalNodes = COUNT( EnthalpyPerm > 0 ) 206 IF ( LocalNodes <= 0 ) RETURN 207 208 SolverParams => GetSolverParams() 209 ConvectionField = GetString( SolverParams, 'Enthalpy_h Convection Field', Found ) 210 211 IF ( Found ) THEN 212 FlowSol => VariableGet( Solver % Mesh % Variables, ConvectionField ) 213 ELSE 214 FlowSol => VariableGet( Solver % Mesh % Variables, 'Flow Solution' ) 215 END IF 216 217 IF ( ASSOCIATED( FlowSol ) ) THEN 218 FlowPerm => FlowSol % Perm 219 NSDOFs = FlowSol % DOFs 220 FlowSolution => FlowSol % Values 221 END IF 222 223 DensitySol => VariableGet( Solver % Mesh % Variables, 'Enthalpy Density' ) 224 225 ! Check whether we have some heater controls. This will affect initialization stuff. 226 SmartHeaterControl = ListCheckPresentAnyBodyForce( Model,'Smart Heater Control') 227 IntegralHeaterControl = ListCheckPresentAnyBodyForce( Model,'Integral Heat Source') 228 229!------------------------------------------------------------------------------ 230! Allocate some permanent storage, this is done first time only 231!------------------------------------------------------------------------------ 232 IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed ) THEN 233 N = Solver % Mesh % MaxElementDOFs 234 235 IF ( AllocationsDone ) THEN 236 DEALLOCATE( & 237 U, V, W, MU, & 238 Pressure, & 239 dPressureDt, & 240 PressureCoeff, & 241 ElementNodes % x, & 242 ElementNodes % y, & 243 ElementNodes % z, & 244 Density,Work, & 245 LatentHeat, & 246 PhaseVelocity, & 247 ElectricConductivity, & 248 Permeability, & 249 Viscosity,C0, & 250 HeatTransferCoeff, & 251 HeatExpansionCoeff, & 252 ReferenceEnthalpy_h, & 253 MASS, & 254 LocalEnthalpy_h, & 255 HeatCapacity,Enthalpy,WaterDiffusivity, & 256 NodalEmissivity, & 257 GasConstant, AText, & 258 HeatConductivity, & 259 STIFF,LOAD, & 260 Indexes, SaveIndexes, & 261 FORCE, TimeForce, & 262 HeatConductivityIso, & 263 PerfusionRate, & 264 PerfusionDensity, & 265 PerfusionHeatCapacity, & 266 PerfusionRefEnthalpy_h ) 267 END IF 268 269 ALLOCATE( & 270 Indexes(N), SaveIndexes(N), & 271 U( N ), V( N ), W( N ), & 272 MU( 3,N ), & 273 Pressure( N ), & 274 dPressureDt( N ), & 275 PressureCoeff( N ), & 276 ElementNodes % x( N ), & 277 ElementNodes % y( N ), & 278 ElementNodes % z( N ), & 279 Density( N ),Work( N ), & 280 LatentHeat( N ), & 281 PhaseVelocity(3, N), & 282 ElectricConductivity(N), & 283 Permeability(N), & 284 Viscosity(N),C0(N), & 285 HeatTransferCoeff( N ), & 286 HeatExpansionCoeff( N ), & 287 ReferenceEnthalpy_h( N ), & 288 MASS( 2*N,2*N ), & 289 LocalEnthalpy_h( N ), & 290 HeatCapacity( N ),Enthalpy( N ), WaterDiffusivity( N ), & 291 NodalEmissivity( N ), & 292 GasConstant( N ),AText( N ), & 293 HeatConductivity( 3,3,N ), & 294 HeatConductivityIso( N ), & 295 STIFF( 2*N,2*N ),LOAD( N ), & 296 FORCE( 2*N ), TimeForce(2*N), & 297 PerfusionRate( N ), & 298 PerfusionDensity( N ), & 299 PerfusionHeatCapacity( N ), & 300 PerfusionRefEnthalpy_h( N ), & 301 STAT=istat ) 302 303 IF ( istat /= 0 ) THEN 304 CALL Fatal( 'HeatSolve', 'Memory allocation error' ) 305 END IF 306 307 308 IF ( SmartHeaterControl .OR. IntegralHeaterControl) THEN 309 n = Model % NumberOfBodyForces 310 IF ( AllocationsDone ) DEALLOCATE( HeaterArea, HeaterDensity, HeaterSource, & 311 HeaterScaling, HeaterTarget, SmartHeaters, IntegralHeaters ) 312 ALLOCATE( HeaterArea(n), HeaterDensity(n), HeaterSource(n), & 313 HeaterScaling(n), HeaterTarget(n), SmartHeaters(n), & 314 IntegralHeaters(n) ) 315 IF ( istat /= 0 ) THEN 316 CALL Fatal( 'HeatSolve', 'Memory allocation error' ) 317 END IF 318 SmartHeaters = .FALSE. 319 IntegralHeaters = .FALSE. 320 END IF 321 322 IF( SmartHeaterControl ) THEN 323 IF ( AllocationsDone ) DEALLOCATE( XX, YY, ForceHeater ) 324 n = SIZE( Enthalpy_h ) 325 ALLOCATE( XX( n ), YY(n), ForceHeater( n ), STAT=istat ) 326 IF ( istat /= 0 ) THEN 327 CALL Fatal( 'HeatSolve', 'Memory allocation error' ) 328 END IF 329 XX = 0.0d0 330 YY = 0.0d0 331 ForceHeater = 0.0d0 332 END IF 333 334 NULLIFY( Hwrk ) 335 AllocationsDone = .TRUE. 336 END IF 337 338!------------------------------------------------------------------------------ 339! Do some additional initialization, and go for it 340!------------------------------------------------------------------------------ 341 dt = Timestep 342 Constants => GetConstants() 343 IF( IsRadiation ) THEN 344 StefanBoltzmann = ListGetConstReal( Model % Constants, & 345 'Stefan Boltzmann' ) 346 END IF 347 348!------------------------------------------------------------------------------ 349 Stabilize = GetLogical( SolverParams,'Stabilize',Found ) 350 351 UseBubbles = GetLogical( SolverParams,'Bubbles',Found ) 352 IF ( .NOT.Found ) UseBubbles = .TRUE. 353 354 StabilizeFlag = GetString( SolverParams, & 355 'Stabilization Method',Found ) 356 357 SELECT CASE(StabilizeFlag) 358 CASE('vms') 359 Stabilize = .FALSE. 360 UseBubbles= .FALSE. 361 CASE('stabilized') 362 Stabilize = .TRUE. 363 UseBubbles = .FALSE. 364 CASE('bubbles') 365 Stabilize = .FALSE. 366 UseBubbles = .TRUE. 367 END SELECT 368 369 NonlinearIter = GetInteger( SolverParams, & 370 'Nonlinear System Max Iterations', Found ) 371 IF ( .NOT.Found ) NonlinearIter = 1 372 373 NonlinearTol = GetConstReal( SolverParams, & 374 'Nonlinear System Convergence Tolerance', Found ) 375 376 IF( IsRadiation ) THEN 377 NewtonTol = GetConstReal( SolverParams, & 378 'Nonlinear System Newton After Tolerance', Found ) 379 NewtonIter = GetInteger( SolverParams, & 380 'Nonlinear System Newton After Iterations', Found ) 381 ELSE 382 NewtonTol = 1.0_dp 383 NewtonIter = 0 384 END IF 385 IF ( NewtonIter == 0) NewtonLinearization = .TRUE. 386 387 Relax = GetCReal( SolverParams, & 388 'Nonlinear System Relaxation Factor',Found ) 389 IF ( .NOT.Found ) Relax = 1 390 391 TransientAssembly = TransientSimulation 392 dt0 = ListGetConstReal(SolverParams,'Steady State Transition Timestep',Found) 393 IF(.NOT. Found) dt0 = ListGetConstReal(SolverParams,'Smart Heater Time Scale',Found) 394 395 IF(Found .AND. dt > dt0) TransientAssembly = .FALSE. 396 397 398!------------------------------------------------------------------------------ 399 400 TransientHeaterControl = .FALSE. 401 IF(SmartHeaterControl) THEN 402 403 ! Mark the smart heaters 404 SmartHeaters = .FALSE. 405 bf_id = 0 406 DO i = 1,Model % NumberOfBodyForces 407 IF( ListGetLogical( Model % BodyForces(i) % Values, & 408 'Smart Heater Control', Found ) ) THEN 409 SmartHeaters(i) = .TRUE. 410 bf_id = i 411 END IF 412 END DO 413 414 ! Find the BC that controls the heater 415 ! If not found assume that smart heater is related to phase change 416 MeltPoint = GetCReal( Model % BodyForces(bf_id) % Values,& 417 'Smart Heater Enthalpy_h',GotMeltPoint) 418 419 SmartHeaterAverage = .FALSE. 420 SmartHeaterNode = ListGetInteger( Model % BodyForces(bf_id) % Values,& 421 'Smart Heater Control Node',GotIt) 422 IF(.NOT. GotIt) THEN 423 RealWork => ListGetConstRealArray( Model % BodyForces(bf_id) % Values,& 424 'Smart Heater Control Point',GotIt) 425 IF( GotIt ) THEN 426 ControlPoint(1:3) = RealWork(1:3,1) 427 428 mindist = HUGE( mindist ) 429 DO l=1,Model % NumberOfNodes 430 IF( EnthalpyPerm(l) == 0 ) CYCLE 431 432 jx = Model % Mesh % Nodes % x(l) 433 jy = Model % Mesh % Nodes % y(l) 434 jz = Model % Mesh % Nodes % z(l) 435 436 dist = (ControlPoint(1)-jx)**2 + (ControlPoint(2)-jy)**2 + (ControlPoint(3)-jz)**2 437 IF( dist < mindist ) THEN 438 mindist = dist 439 SmartHeaterNode = l 440 END IF 441 END DO 442 END IF 443 444 WRITE(Message,*) 'Found Control Point at distance:',SQRT(mindist) 445 CALL Info('HeatSolve',Message) 446 WRITE(Message,*) 'Control Point Index:',SmartHeaterNode 447 CALL Info('HeatSolve',Message) 448 END IF 449 450 IF( .NOT. GotMeltPoint .OR. SmartHeaterNode == 0) THEN 451 GotIt = .FALSE. 452 Found = .FALSE. 453 SmartHeaterBC = 0 454 455 DO i=1,Model % NumberOfBCs 456 GotIt = ListGetLogical( Model % BCs(i) % Values,'Smart Heater Boundary', Found ) 457 IF(GotIt) THEN 458 SmartHeaterBC = i 459 EXIT 460 END IF 461 END DO 462 IF(.NOT. GotIt) THEN 463 DO i=1,Model % NumberOfBCs 464 GotIt = ListGetLogical( Model % BCs(i) % Values,'Phase Change', Found ) 465 IF(GotIt) THEN 466 SmartHeaterBC = i 467 EXIT 468 END IF 469 END DO 470 END IF 471 IF(SmartHeaterBC == 0) THEN 472 CALL Fatal('HeatSolve','Smart Heater Boundary / Phase Change is undefined') 473 END IF 474 475 MeltPoint = GetCReal( Model % BCs(SmartHeaterBC) % Values,& 476 'Smart Heater Enthalpy_h',Found) 477 IF(.NOT. Found) THEN 478 DO k=1, Model % NumberOfMaterials 479 MeltPoint = GetCReal( Model % Materials(k) % Values, & 480 'Melting Point', Found ) 481 IF(Found) EXIT 482 END DO 483 IF(.NOT. Found) THEN 484 CALL Fatal('HeatSolver','Smart Heater Enthalpy_h / Melting Point is undefined') 485 END IF 486 END IF 487 488 ! Find the node related to Enthalpy_h control 489 SmartHeaterAverage = ListGetLogical(Solver % Values,'Smart Heater Average', Found) 490 IF(.NOT. SmartHeaterAverage) THEN 491 jx = -HUGE(jx) 492 DO k = Model % Mesh % NumberOfBulkElements + 1, & 493 Model % Mesh % NumberOfBulkElements + Model % Mesh % NumberOfBoundaryElements 494 495 Element => Model % Mesh % Elements(k) 496 497 IF ( Element % BoundaryInfo % Constraint == SmartHeaterBC ) THEN 498 DO l=1,Element % TYPE % NumberOfNodes 499 IF ( Model % Mesh % Nodes % x(Element % NodeIndexes(l)) >= jx ) THEN 500 j = Element % NodeIndexes(l) 501 jx = Model % Mesh % Nodes % x(Element % NodeIndexes(l)) 502 END IF 503 END DO 504 END IF 505 END DO 506 SmartHeaterNode = j 507 END IF 508 END IF 509 510 SmartTol = GetConstReal( SolverParams, & 511 'Smart Heater Control After Tolerance', Found ) 512 IF(.NOT. Found) THEN 513 SmartTolReached = .TRUE. 514 SmartTol = 1.0 515 END IF 516 517 PowerTimeScale = ListGetConstReal(Solver % Values, & 518 'Smart Heater Time Scale',Found) 519 520 IF(TransientSimulation .AND. dt < PowerTimeScale) THEN 521 TransientHeaterControl = .TRUE. 522 CALL Info( 'HeatSolve', 'Using Transient Heater Control') 523 ELSE 524 TransientHeaterControl = .FALSE. 525 CALL Info( 'HeatSolve', 'Using Steady-state Heater Control') 526 END IF 527 528 IF(Solver % DoneTime /= DoneTime) THEN 529 PrevPowerScaling = PowerScaling 530 DoneTime = Solver % DoneTime 531 END IF 532 END IF 533 534 IF( IntegralHeaterControl) THEN 535 CALL Info( 'HeatSolve', 'Using Integral Heater Control') 536 IntegralHeaters = .FALSE. 537 DO i = 1,Model % NumberOfBodyForces 538 IntegralHeaters(i) = ListCheckPresent( Model % BodyForces(i) % Values, & 539 'Integral Heat Source') 540 END DO 541 END IF 542 543!------------------------------------------------------------------------------ 544 545 ConstantBulk = GetLogical( SolverParams, 'Constant Bulk System', Found ) 546 SaveBulk = ConstantBulk .OR. GetLogical( SolverParams, 'Save Bulk System', Found ) 547 SaveBulk = ConstantBulk .OR. GetLogical( SolverParams, 'Calculate Loads', Found ) 548 549!------------------------------------------------------------------------------ 550 551 SaveRelax = Relax 552 CumulativeTime = 0.0d0 553 554!------------------------------------------------------------------------------ 555 FirstTime = .TRUE. 556 ALLOCATE(PrevSolution(LocalNodes)) 557 558 DO WHILE( CumulativeTime < Timestep-1.0d-12 .OR. .NOT. TransientSimulation ) 559!------------------------------------------------------------------------------ 560! The first time around this has been done by the caller... 561!------------------------------------------------------------------------------ 562 IF ( TransientSimulation .AND. .NOT.FirstTime ) & 563 CALL InitializeTimestep(Solver) 564 FirstTime = .FALSE. 565!------------------------------------------------------------------------------ 566! Save current solution 567!------------------------------------------------------------------------------ 568 PrevSolution = Enthalpy_h(1:LocalNodes) 569 IF ( TransientSimulation ) THEN 570 PrevEnthalpy_h => Solver % Variable % PrevValues(:,1) 571 END IF 572!------------------------------------------------------------------------------ 573 574 totat = 0.0d0 575 totst = 0.0d0 576 577 Norm = Solver % Variable % Norm 578 579 580!Calculate Phase change enthalpy ================================================================== 581!================================================================================================== 582 583 A_cap = GetConstReal(Model % Constants, "Enthalpy Heat Capacity A",GotIt) 584 IF (.NOT.GotIt) THEN 585 CALL WARN('EnthalpySolver', 'No Keyword >Enthalpy Heat Capacity A< & 586 & defined in model constants. Using >7.253< as default (Paterson, 1994).') 587 A_cap = 7.253 588 END IF 589 590 B_cap = GetConstReal(Model % Constants, "Enthalpy Heat Capacity B",GotIt) 591 IF (.NOT.GotIt) THEN 592 CALL WARN('EnthalpySolver', 'No Keyword >Enthalpy Heat Capacity B< & 593 & defined in model constants. Using >146.3< as default (Paterson, 1994).') 594 B_cap = 146.3 595 END IF 596 597 PressureName = GetString(Constants,'Pressure Variable', GotIt) 598 IF (.NOT.GotIt) THEN 599 CALL WARN('EnthalpySolver', 'No Keyword >Pressure Variable< defined. Using >Pressure< as default.') 600 WRITE(PressureName,'(A)') 'Pressure' 601 ELSE 602 WRITE(Message,'(a,a)') 'Variable Name for pressure: ', PressureName 603 CALL INFO('EnthalpySolve',Message,Level=12) 604 END IF 605 606 Tref = GetConstReal(Model % Constants, "T_ref_enthalpy",GotIt) 607 IF (.NOT.GotIt) THEN 608 CALL WARN('EnthalpySolver', 'No Keyword >Reference Enthalpy< defined. Using >200K< as default.') 609 Tref = 200.0 610 END IF 611 612 beta = GetConstReal(Model % Constants, "beta_clapeyron",GotIt) 613 IF (.NOT.GotIt) THEN 614 CALL WARN('EnthalpySolver', 'No Keyword >beta_clapeyron< defined. Using >9.74 10-2 K Mpa-1< as default.') 615 beta = 0.0974 616 END IF 617 618 Psurf = GetConstReal(Model % Constants, "P_surf",GotIt) 619 IF (.NOT.GotIt) THEN 620 CALL WARN('EnthalpySolver', 'No Keyword >Psurf< defined. Using >1.013 10-1 MPa < as default.') 621 Psurf = 0.1013 622 END IF 623 624 Ptriple = GetConstReal(Model % Constants,"P_triple",GotIt) 625 IF (.NOT.GotIt) THEN 626 CALL WARN('EnthalpySolver', 'No Keyword >P_triple< defined. Using >0.061173 MPa < as default.') 627 Ptriple = 0.061173 628 END IF 629 630 PressureVariable => VariableGet(Solver % Mesh %Variables,PressureName) 631 IF ( ASSOCIATED( PressureVariable ) ) THEN 632 PressurePerm => PressureVariable % Perm 633 PressureValues => PressureVariable % Values 634 END IF 635 636 PhaseEnthName = GetString(SolverParams , 'Exported Variable 1', Found ) 637 638 IF (.NOT.Found) & 639 CALL FATAL('EnthalpySolver','No value > Exported Variable 1 (enthalpy phase change) < found in Solver') 640 641 WaterVar => VariableGet(Solver % Mesh % Variables, WaterName) 642 TemphomoVar => VariableGet(Solver % Mesh % Variables,TempName) 643 644 PhaseChangeEnthVar => & 645 VariableGet(Solver % Mesh % Variables,PhaseEnthName) 646 IF ( ASSOCIATED( PhaseChangeEnthVar ) ) THEN 647 PhaseChangeEnthPerm => PhaseChangeEnthVar % Perm 648 PhaseChangeEnthValues => PhaseChangeEnthVar % Values 649 END IF 650 651do i=1,model % NumberOfNodes 652 P=PressureValues (PressurePerm(i))+Psurf 653 Tm=273.16-beta*(P-Ptriple)! Clapeyron 654 PhaseChangeEnthValues (PhaseChangeEnthPerm(i)) = A_cap*0.5*Tm*Tm+B_cap*Tm-Tref*(A_cap*0.5*Tref+B_cap) ! use CP(T)=B_cap+A_cap*T 655enddo 656 657!===================================================================================================================== 658!===================================================================================================================== 659 660 661 DO iter=1,NonlinearIter 662 at = CPUTime() 663 at0 = RealTime() 664 665 CALL Info( 'HeatSolve', ' ', Level=4 ) 666 CALL Info( 'HeatSolve', ' ', Level=4 ) 667 CALL Info( 'HeatSolve', '-------------------------------------',Level=4 ) 668 WRITE( Message,* ) 'Enthalpy_h ITERATION', iter 669 CALL Info( 'HeatSolve', Message, Level=4 ) 670 CALL Info( 'HeatSolve', '-------------------------------------',Level=4 ) 671 CALL Info( 'HeatSolve', ' ', Level=4 ) 672 CALL Info( 'HeatSolve', 'Starting Assembly...', Level=4 ) 673 674 IF ( ConstantBulk .AND. ASSOCIATED(Solver % Matrix % BulkValues) ) THEN 675 Solver % Matrix % Values = Solver % Matrix % BulkValues 676 Solver % Matrix % RHS = Solver % Matrix % BulkRHS 677 GOTO 1000 678 END IF 679 680!------------------------------------------------------------------------------ 681 CALL DefaultInitialize() 682!------------------------------------------------------------------------------ 683 684 685 686 687 IF ( SmartHeaterControl .OR. IntegralHeaterControl ) THEN 688 IF( SmartHeaterControl) ForceHeater = 0.0d0 689 HeaterArea = 0.0d0 690 HeaterSource = 0.0d0 691 HeaterScaling = 1.0d0 692 HeaterDensity = 0.0d0 693 HeaterTarget = 0.0d0 694 HeaterControlLocal = .FALSE. 695 696 DO t=1,Solver % NumberOfActiveElements 697 Element => GetActiveElement(t) 698 BodyForce => GetBodyForce() 699 700 IF ( .NOT. ASSOCIATED( BodyForce ) ) CYCLE 701 bf_id = GetBodyForceId() 702 703 IF( .NOT. (SmartHeaters(bf_id) .OR. IntegralHeaters(bf_id) ) ) CYCLE 704 705 n = GetElementNOFNodes() 706 707 Material => GetMaterial() 708 709 Density(1:n) = GetReal( Material, 'Enthalpy Density' ) 710 Load(1:n) = GetReal( BodyForce, 'Heat Source' ) 711 712 s = ElementArea( Solver % Mesh, Element, n ) 713 714 IF( CurrentCoordinateSystem() == AxisSymmetric .OR. & 715 CurrentCoordinateSystem() == CylindricSymmetric ) s = 2 * PI * s 716 717 HeaterSource(bf_id) = HeaterSource(bf_id) + s * SUM(Density(1:n) * Load(1:n)) / n 718 HeaterArea(bf_id) = HeaterArea(bf_id) + s 719 HeaterDensity(bf_id) = HeaterDensity(bf_id) + s * SUM( Density(1:n) ) / n 720 END DO 721 722 DO i = 1,Model % NumberOfBodyForces 723 IF( IntegralHeaters(i) .OR. SmartHeaters(i) ) THEN 724 HeaterDensity(i) = HeaterDensity(i) / HeaterArea(i) 725 END IF 726 IF(IntegralHeaters(i)) THEN 727 HeaterTarget(i) = GetCReal( Model % BodyForces(i) % Values, & 728 'Integral Heat Source', Found ) 729 HeaterScaling(i) = HeaterTarget(i) / HeaterSource(i) 730 END IF 731 END DO 732 END IF 733 734!------------------------------------------------------------------------------ 735 body_id = -1 736 NULLIFY(Material) 737!------------------------------------------------------------------------------ 738! Bulk elements 739!------------------------------------------------------------------------------ 740 CALL StartAdvanceOutput( 'HeatSolve', 'Assembly:' ) 741 DO t=1,Solver % NumberOfActiveElements 742 743 CALL AdvanceOutput(t,GetNOFActive()) 744!------------------------------------------------------------------------------ 745! Check if this element belongs to a body where Enthalpy_h 746! should be calculated 747!------------------------------------------------------------------------------ 748 Element => GetActiveElement(t) 749 750!------------------------------------------------------------------------------ 751 IF ( Element % BodyId /= body_id ) THEN 752!------------------------------------------------------------------------------ 753 Equation => GetEquation() 754 ConvectionFlag = GetString( Equation, 'Convection', Found ) 755 756 Material => GetMaterial() 757!------------------------------------------------------------------------------ 758 CompressibilityFlag = GetString( Material, & 759 'Compressibility Model', Found) 760 IF ( .NOT.Found ) CompressibilityModel = Incompressible 761 762 SELECT CASE( CompressibilityFlag ) 763 764 CASE( 'incompressible' ) 765 CompressibilityModel = Incompressible 766 767 CASE( 'user defined' ) 768 CompressibilityModel = UserDefined1 769 770 CASE( 'perfect gas', 'perfect gas equation 1' ) 771 CompressibilityModel = PerfectGas1 772 773 CASE( 'thermal' ) 774 CompressibilityModel = Thermal 775 776 CASE DEFAULT 777 CompressibilityModel = Incompressible 778 END SELECT 779!------------------------------------------------------------------------------ 780 781 PhaseModel = GetString( Equation, 'Phase Change Model',Found ) 782 IF(.NOT. Found) PhaseModel = GetString( Material, 'Phase Change Model',Found ) 783 784 PhaseChange = Found .AND. (PhaseModel(1:4) /= 'none') 785 IF ( PhaseChange ) THEN 786 CheckLatentHeatRelease = GetLogical( Equation, & 787 'Check Latent Heat Release',Found ) 788 END IF 789 END IF 790!------------------------------------------------------------------------------ 791 792 n = GetElementNOFNodes() 793 CALL GetElementNodes( ElementNodes ) 794 795 CALL GetScalarLocalSolution( LocalEnthalpy_h ) 796!------------------------------------------------------------------------------ 797! Get element material parameters 798!------------------------------------------------------------------------------ 799 HeatCapacity(1:n) = 1.0 800 WaterDiffusivity(1:n) = GetReal( Material, 'Enthalpy Water Diffusivity', Found ) 801 802IF (.NOT.Found) THEN 803 CALL FATAL('EnthalpySolver', 'No value for >Enthalpy Water Diffusivity< defined in material.') 804ENDIF 805 806 807 CALL ListGetRealArray( Material,'Enthalpy Heat Diffusivity',Hwrk,n, & 808 Element % NodeIndexes ) 809 810 HeatConductivity = 0.0d0 811 IF ( SIZE(Hwrk,1) == 1 ) THEN 812 DO i=1,3 813 DO j=1,n 814 IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(j)))<& 815 &Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(j)))) then 816 HeatConductivity( i,i,j ) = WaterDiffusivity(j) 817 ELSE 818 HeatConductivity( i,i,j ) = Hwrk( 1,1,j ) 819 ENDIF 820 ENDDO 821 END DO 822 ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN 823 DO i=1,MIN(3,SIZE(Hwrk,1)) 824 DO j=1,n 825 IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(j)))<& 826 &Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(j)))) then 827 HeatConductivity(i,i,j) = WaterDiffusivity(j) 828 ELSE 829 HeatConductivity(i,i,j) = Hwrk(i,1,j) 830 ENDIF 831 ENDDO 832 END DO 833 ELSE 834 DO i=1,MIN(3,SIZE(Hwrk,1)) 835 DO j=1,MIN(3,SIZE(Hwrk,2)) 836 DO k=1,n 837 IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(k)))<& 838 &Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(k)))) then 839 HeatConductivity( i,j,k ) = WaterDiffusivity(k) 840 ELSE 841 HeatConductivity( i,j,k ) = Hwrk(i,j,k) 842 ENDIF 843 844 ENDDO 845 END DO 846 END DO 847 END IF 848!------------------------------------------------------------------------------ 849 850 IF ( CompressibilityModel == PerfectGas1 ) THEN 851 852 ! Read Specific Heat Ratio: 853 !-------------------------- 854 SpecificHeatRatio = GetConstReal( Material, & 855 'Specific Heat Ratio', Found ) 856 IF ( .NOT.Found ) SpecificHeatRatio = 5.d0/3.d0 857 858 ! For an ideal gas, \gamma, c_p and R are really a constant 859 ! GasConstant is an array only since HeatCapacity formally is: 860 !------------------------------------------------------------- 861 GasConstant(1:n) = ( SpecificHeatRatio - 1.d0 ) * & 862 HeatCapacity(1:n) / SpecificHeatRatio 863 864 PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found ) 865 IF ( .NOT. Found ) PressureCoeff(1:n) = 1.0_dp 866 ELSE IF ( CompressibilityModel == Thermal ) THEN 867 ReferenceEnthalpy_h(1:n) = GetReal( Material, 'Reference Enthalpy_h' ) 868 HeatExpansionCoeff(1:n) = GetReal( Material, 'Heat Expansion Coefficient' ) 869 870 Density(1:n) = GetReal( Material,'Enthalpy Density' ) 871 Density(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n) * & 872 ( LocalEnthalpy_h(1:n) - ReferenceEnthalpy_h(1:n) ) ) 873 874 PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found ) 875 IF ( .NOT. Found ) & 876 PressureCoeff(1:n) = LocalEnthalpy_h(1:n) * HeatExpansionCoeff(1:n) / & 877 ( 1-HeatExpansionCoeff(1:n)*( & 878 LocalEnthalpy_h(1:n)-ReferenceEnthalpy_h(1:n)) ) 879 ELSE IF ( CompressibilityModel == UserDefined1 ) THEN 880 IF ( ASSOCIATED( DensitySol ) ) THEN 881 CALL GetScalarLocalSolution( Density, 'Enthalpy Density' ) 882 ELSE 883 Density(1:n) = GetReal( Material,'Enthalpy Density' ) 884 END IF 885 PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found ) 886 IF ( .NOT. Found ) PressureCoeff(1:n) = 0.0_dp 887 ELSE 888 PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found ) 889 IF ( .NOT. Found ) PressureCoeff(1:n) = 0.0_dp 890 Density(1:n) = GetReal( Material, 'Enthalpy Density' ) 891 END IF 892 893!------------------------------------------------------------------------------ 894! Take pressure deviation p_d as the dependent variable, p = p_0 + p_d 895! for PerfectGas, read p_0 896!------------------------------------------------------------------------------ 897 IF ( CompressibilityModel /= Incompressible ) THEN 898 ReferencePressure = ListGetConstReal( Material, & 899 'Reference Pressure', Found,UnFoundFatal=UnFoundFatal) 900 !Previous default value: ReferencePressure = 0.0d0 901 END IF 902!------------------------------------------------------------------------------ 903 904 HeaterControlLocal = .FALSE. 905 Load = 0.0D0 906 Pressure = 0.0d0 907 dPressuredt = 0.0d0 908!------------------------------------------------------------------------------ 909! Check for convection model 910!------------------------------------------------------------------------------ 911 C1 = 1.0D0 912 U = 0._dp 913 V = 0._dp 914 W = 0._dp 915 916 MU = 0.0d0 917 CALL GetVectorLocalSolution( MU, 'Mesh Velocity' ) 918 919 IF ( ConvectionFlag == 'constant' ) THEN 920 921 U(1:n) = GetReal( Material, 'Convection Velocity 1', Found ) 922 IF ( .NOT. Found ) & 923 U(1:n) = GetReal( Equation, 'Convection Velocity 1', Found ) 924 V(1:n) = GetReal( Material, 'Convection Velocity 2', Found ) 925 IF ( .NOT. Found ) & 926 V(1:n) = GetReal( Equation, 'Convection Velocity 2', Found ) 927 W(1:n) = GetReal( Material, 'Convection Velocity 3', Found ) 928 IF ( .NOT. Found ) & 929 W(1:n) = GetReal( Equation, 'Convection Velocity 3', Found ) 930 931 ELSE IF ( ConvectionFlag == 'computed' .AND. & 932 ASSOCIATED(FlowSolution) ) THEN 933 DO i=1,n 934 k = FlowPerm(Element % NodeIndexes(i)) 935 IF ( k > 0 ) THEN 936!------------------------------------------------------------------------------ 937 Pressure(i) = FlowSolution(NSDOFs*k) + ReferencePressure 938 SELECT CASE( CompressibilityModel ) 939 CASE( PerfectGas1 ) 940 Density(i) = Pressure(i) / & 941 ( GasConstant(i) * LocalEnthalpy_h(i) ) 942 END SELECT 943 IF ( TransientSimulation ) THEN 944 dPressureDt(i) = ( FlowSolution(NSDOFs*k) - & 945 FlowSol % PrevValues(NSDOFs*k,1) ) / dt 946 END IF 947!------------------------------------------------------------------------------ 948 949 SELECT CASE( NSDOFs ) 950 CASE(3) 951 U(i) = FlowSolution( NSDOFs*k-2 ) 952 V(i) = FlowSolution( NSDOFs*k-1 ) 953 W(i) = 0.0D0 954 955 CASE(4) 956 U(i) = FlowSolution( NSDOFs*k-3 ) 957 V(i) = FlowSolution( NSDOFs*k-2 ) 958 W(i) = FlowSolution( NSDOFs*k-1 ) 959 END SELECT 960 ELSE 961 U(i) = 0.0d0 962 V(i) = 0.0d0 963 W(i) = 0.0d0 964 END IF 965 END DO 966 ELSE IF (ConvectionFlag=='computed' ) THEN 967 CALL Warn( 'HeatSolver', 'Convection model specified ' // & 968 'but no associated flow field present?' ) 969 ELSE 970 IF ( ALL(MU==0) ) C1 = 0.0D0 971 END IF 972!------------------------------------------------------------------------------ 973! Check if modelling Phase Change with Eulerian approach 974!------------------------------------------------------------------------------ 975 PhaseSpatial = .FALSE. 976 IF ( PhaseChange ) THEN 977 CALL EffectiveHeatCapacity() 978 ELSE 979 HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n) 980 END IF 981 982 Viscosity = 0.0d0 983!------------------------------------------------------------------------------ 984! Add body forces, if any 985!------------------------------------------------------------------------------ 986 BodyForce => GetBodyForce() 987 IF ( ASSOCIATED( BodyForce ) ) THEN 988 bf_id = GetBodyForceId() 989!------------------------------------------------------------------------------ 990! Frictional viscous heating 991!------------------------------------------------------------------------------ 992 IF ( GetLogical( BodyForce, 'Friction Heat',Found) ) THEN 993 Viscosity(1:n) = GetReal( Material,'Viscosity' ) 994 END IF 995!------------------------------------------------------------------------------ 996! Given heat source 997!------------------------------------------------------------------------------ 998 Load(1:n) = Density(1:n) * GetReal( BodyForce, 'Heat Source', Found ) 999 1000 IF ( SmartHeaterControl .AND. NewtonLinearization .AND. SmartTolReached) THEN 1001 IF( SmartHeaters(bf_id) ) THEN 1002 HeaterControlLocal = .TRUE. 1003 IF( TransientHeaterControl ) THEN 1004 Load(1:n) = PrevPowerScaling * Load(1:n) 1005 HeaterScaling(bf_id) = PrevPowerScaling 1006 END IF 1007 END IF 1008 END IF 1009 1010 IF ( IntegralHeaterControl ) THEN 1011 IF( IntegralHeaters(bf_id) ) THEN 1012 Load(1:n) = Load(1:n) * HeaterScaling(bf_id) 1013 END IF 1014 END IF 1015 END IF 1016 1017 C0 = 0.0_dp 1018!------------------------------------------------------------------------------ 1019! Note at this point HeatCapacity = \rho * c_p OR \rho * (c_p - R) 1020! and C1 = 0 (diffusion) or 1 (convection) 1021!------------------------------------------------------------------------------ 1022 1023!------------------------------------------------------------------------------ 1024! Perfusion (added as suggested by Matthias Zenker) 1025!------------------------------------------------------------------------------ 1026 PerfusionRate(1:n) = GetReal( BodyForce, 'Perfusion Rate', Found ) 1027 1028 IF ( Found ) THEN 1029 PerfusionRefEnthalpy_h(1:n) = GetReal( BodyForce, 'Perfusion Reference Enthalpy_h' ) 1030 PerfusionDensity(1:n) = GetReal( BodyForce, 'Perfusion Density' ) 1031 PerfusionHeatCapacity(1:n) = GetReal( BodyForce, 'Perfusion Heat Capacity' ) 1032 C0(1:n) = PerfusionHeatCapacity(1:n) * PerfusionRate(1:n) * PerfusionDensity(1:n) 1033 Load(1:n) = Load(1:n) + C0(1:n) * PerfusionRefEnthalpy_h(1:n) 1034 END IF 1035 1036!------------------------------------------------------------------------------ 1037! Get element local matrices, and RHS vectors 1038!------------------------------------------------------------------------------ 1039 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1040!------------------------------------------------------------------------------ 1041 CALL DiffuseConvectiveCompose( & 1042 MASS, STIFF, FORCE, LOAD, & 1043 HeatCapacity, C0, C1*HeatCapacity(1:n), HeatConductivity, & 1044 PhaseSpatial, LocalEnthalpy_h, Enthalpy, U, V, W, & 1045 MU(1,1:n),MU(2,1:n),MU(3,1:n), Viscosity, Density, Pressure, & 1046 dPressureDt, PressureCoeff, CompressibilityModel /= Incompressible, & 1047 Stabilize, UseBubbles, Element, n, ElementNodes ) 1048 1049!------------------------------------------------------------------------------ 1050 ELSE 1051!------------------------------------------------------------------------------ 1052 CALL DiffuseConvectiveGenCompose( & 1053 MASS, STIFF, FORCE, LOAD, & 1054 HeatCapacity, C0, C1*HeatCapacity(1:n), HeatConductivity, & 1055 PhaseSpatial, LocalEnthalpy_h, Enthalpy, U, V, W, & 1056 MU(1,1:n),MU(2,1:n),MU(3,1:n), Viscosity, Density, Pressure, & 1057 dPressureDt, PressureCoeff, CompressibilityModel /= Incompressible, & 1058 Stabilize, Element, n, ElementNodes ) 1059!------------------------------------------------------------------------------ 1060 END IF 1061!------------------------------------------------------------------------------ 1062 1063 IF ( HeaterControlLocal .AND. .NOT. TransientHeaterControl) THEN 1064 1065 IF ( TransientAssembly .AND. .NOT. ConstantBulk ) THEN 1066 CALL Default1stOrderTime( MASS, STIFF, FORCE ) 1067 END IF 1068 1069 CALL UpdateGlobalEquations( Solver % Matrix, STIFF, & 1070 ForceHeater, FORCE, n, 1, EnthalpyPerm(Element % NodeIndexes) ) 1071 ELSE 1072 Bubbles = UseBubbles .AND. .NOT.Stabilize .AND. & 1073 ( ConvectionFlag == 'computed' .OR. ConvectionFlag == 'constant' ) 1074 1075!------------------------------------------------------------------------------ 1076! If time dependent simulation add mass matrix to stiff matrix 1077!------------------------------------------------------------------------------ 1078 TimeForce = 0.0_dp 1079 IF ( TransientAssembly ) THEN 1080 IF ( ConstantBulk ) THEN 1081 CALL DefaultUpdateMass( MASS ) 1082 ELSE 1083 CALL Default1stOrderTime( MASS,STIFF,FORCE ) 1084 END IF 1085 ELSE IF ( Solver % NOFEigenValues>0 ) THEN 1086 CALL DefaultUpdateDamp(MASS) 1087 END IF 1088!------------------------------------------------------------------------------ 1089! Update global matrices from local matrices 1090!------------------------------------------------------------------------------ 1091 IF ( Bubbles ) THEN 1092 CALL Condensate( N, STIFF, FORCE, TimeForce ) 1093 END IF 1094 1095 CALL DefaultUpdateEquations( STIFF, FORCE ) 1096 END IF 1097!------------------------------------------------------------------------------ 1098 END DO ! Bulk elements 1099!------------------------------------------------------------------------------ 1100 1101 CALL DefaultFinishBulkAssembly() 1102 1103 11041000 CONTINUE 1105 1106!------------------------------------------------------------------------------ 1107! Neumann & Newton boundary conditions 1108!------------------------------------------------------------------------------ 1109 DO t=1, Solver % Mesh % NumberOfBoundaryElements 1110 Element => GetBoundaryElement(t) 1111 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 1112 1113 n = GetElementNOFNodes() 1114 1115 ! Check that the dimension of element is suitable for fluxes 1116 IF( .NOT. PossibleFluxElement(Element) ) CYCLE 1117 1118 BC => GetBC() 1119 IF ( .NOT. ASSOCIATED(BC) ) CYCLE 1120 1121 ! This checks whether there are any Dirichlet conditions on the 1122 ! smart heater boundary. If there are the r.h.s. must be zero as 1123 ! there can possibly not be any effect on Enthalpy_h. 1124 !----------------------------------------------------------------- 1125 IF ( HeaterControlLocal .AND. .NOT. TransientHeaterControl) THEN 1126 IF( ListCheckPresent(BC, Varname) ) THEN 1127 nd = GetElementDOFs(Indexes) 1128 ForceHeater(EnthalpyPerm(Indexes(1:nd))) = 0.0_dp 1129 END IF 1130 END IF 1131 1132 HeatFluxBC = GetLogical( BC, 'Enthalpy Heat Flux BC', Found ) 1133 IF ( Found .AND. .NOT. HeatFluxBC ) CYCLE 1134 1135 HeatGapBC = ListGetLogical( BC, 'Heat Gap', Found ) 1136 CALL AddHeatFluxBC() 1137 1138 IF ( HeatGapBC ) THEN 1139 CALL FindGapIndexes( Element, Indexes, n ) 1140 SaveIndexes(1:n) = Element % NodeIndexes 1141 Element % NodeIndexes = Indexes(1:n) 1142 CALL AddHeatFluxBC() 1143 Element % NodeIndexes = SaveIndexes(1:n) 1144 END IF 1145 1146 END DO ! Neumann & Newton BCs 1147!------------------------------------------------------------------------------ 1148 1149 IF ( TransientSimulation .AND. ConstantBulk ) CALL AddGlobalTime() 1150 1151 CALL DefaultFinishAssembly() 1152 CALL Info( 'HeatSolve', 'Assembly done', Level=4 ) 1153 1154 CALL DefaultDirichletBCs() 1155 1156!------------------------------------------------------------------------------ 1157 1158!------------------------------------------------------------------------------ 1159! Solve the system and check for convergence 1160!------------------------------------------------------------------------------ 1161 at = CPUTime() - at 1162 st = CPUTime() 1163 1164 PrevNorm = Norm 1165 1166 IF(SmartHeaterControl .AND. NewtonLinearization .AND. SmartTolReached) THEN 1167 1168 IF(.NOT. TransientHeaterControl) THEN 1169 1170 CALL ListAddLogical(SolverParams, & 1171 'Skip Compute Nonlinear Change',.TRUE.) 1172 1173 Relax = GetCReal( SolverParams, & 1174 'Nonlinear System Relaxation Factor', Found ) 1175 1176 IF ( Found .AND. Relax /= 1.0d0 ) THEN 1177 CALL ListAddConstReal( Solver % Values, & 1178 'Nonlinear System Relaxation Factor', 1.0d0 ) 1179 ELSE 1180 Relax = 1.0d0 1181 END IF 1182 1183 CALL SolveSystem( Solver % Matrix, ParMatrix, & 1184 ForceHeater, XX, Norm, 1, Solver ) 1185 1186 CALL SolveSystem( Solver % Matrix, ParMatrix, & 1187 Solver % Matrix % RHS, YY, Norm, 1, Solver ) 1188 1189 CALL ListAddLogical(SolverParams,'Skip Compute Nonlinear Change',.FALSE.) 1190 ELSE 1191 CALL SolveSystem( Solver % Matrix, ParMatrix, & 1192 Solver % Matrix % RHS, Enthalpy_h, Norm, 1, Solver ) 1193 YY = Enthalpy_h 1194 END IF 1195 1196 IF(.NOT. SmartHeaterAverage) THEN 1197 xave = XX(EnthalpyPerm(SmartHeaterNode)) 1198 yave = YY(EnthalpyPerm(SmartHeaterNode)) 1199 ELSE 1200 xave = 0.0d0 1201 yave = 0.0d0 1202 j = 0 1203 1204 DO k = Model % Mesh % NumberOfBulkElements + 1, & 1205 Model % Mesh % NumberOfBulkElements + Model % Mesh % NumberOfBoundaryElements 1206 1207 Element => Model % Mesh % Elements(k) 1208 IF ( Element % BoundaryInfo % Constraint == SmartHeaterBC ) THEN 1209 l = Element % TYPE % NumberOfNodes 1210 j = j + l 1211 xave = xave + SUM( XX(EnthalpyPerm(Element % NodeIndexes)) ) 1212 yave = yave + SUM( YY(EnthalpyPerm(Element % NodeIndexes)) ) 1213 END IF 1214 END DO 1215 xave = xave / j 1216 yave = yave / j 1217 CALL ListAddConstReal(Model % Simulation,'res: Smart Heater Enthalpy_h',yave) 1218 END IF 1219 1220 IF(.NOT. TransientHeaterControl) THEN 1221 IF ( ASSOCIATED(Solver % Variable % NonlinValues) ) THEN 1222 Solver % Variable % NonlinValues = Enthalpy_h 1223 END IF 1224 1225 PowerScaling = (MeltPoint - yave) / xave 1226 Enthalpy_h = YY + PowerScaling * XX 1227 1228 ! The change is computed separately for the controlled Enthalpy_h field 1229 !----------------------------------------------------------------------- 1230 CALL ComputeChange(Solver,.FALSE.,LocalNodes,Enthalpy_h) 1231 Norm = Solver % Variable % Norm 1232 1233 END IF 1234 1235 IF(dt > PowerTimeScale) THEN 1236 IF ( Relax /= 1.0d0 ) THEN 1237 CALL ListAddConstReal( Solver % Values, & 1238 'Nonlinear System Relaxation Factor', Relax ) 1239 END IF 1240 END IF 1241 ELSE 1242 Norm = DefaultSolve() 1243 END IF 1244 1245 1246 IF( SmartHeaterControl .OR. IntegralHeaterControl) THEN 1247 1248 CALL ListAddConstReal(Model % Simulation,'res: Heater Power Scaling',PowerScaling) 1249 1250 CALL Info( 'HeatSolve', 'Heater Control Information', Level=4 ) 1251 DO i=1,Model % NumberOfBodyForces 1252 IF( .NOT. (SmartHeaters(i) .OR. IntegralHeaters(i))) CYCLE 1253 IF( SmartHeaters(i) ) HeaterScaling(i) = PowerScaling 1254 1255 WRITE( Message, '(A,T35,I15)' ) 'Heater for body: ', i 1256 CALL Info( 'HeatSolve', Message, Level=4 ) 1257 IF(SmartHeaters(i)) WRITE( Message, '(A,T35,A)' ) 'Heater type:','Smart heater' 1258 IF(IntegralHeaters(i)) WRITE( Message, '(A,T35,A)' ) 'Heater type:','Integral heater' 1259 CALL Info( 'HeatSolve', Message, Level=4 ) 1260 1261 WRITE( Message,'(A,T35,ES15.4)' ) 'Heater Volume (m^3): ', HeaterArea(i) 1262 CALL Info( 'HeatSolve', Message, Level=4 ) 1263 s = HeaterSource(i) * HeaterScaling(i) 1264 WRITE( Message,'(A,T35,ES15.4)' ) 'Heater Power (W): ', s 1265 CALL Info( 'HeatSolve', Message, Level=4 ) 1266 1267 WRITE( Message,'(A,T35,ES15.4)' ) 'Heater scaling: ', HeaterScaling(i) 1268 CALL Info( 'HeatSolve', Message, Level=4 ) 1269 WRITE( Message, '(A,T35,ES15.4)' ) 'Heater Power Density (W/kg): ', s/(HeaterDensity(i) * HeaterArea(i)) 1270 CALL Info( 'HeatSolve', Message, Level=4 ) 1271 1272 IF( SmartHeaters(i)) CALL ListAddConstReal(Model % Simulation,'res: Heater Power Density',& 1273 s/(HeaterDensity(i) * HeaterArea(i))) 1274 END DO 1275 END IF 1276 1277 1278 st = CPUTIme()-st 1279 totat = totat + at 1280 totst = totst + st 1281 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat 1282 CALL Info( 'HeatSolve', Message, Level=4 ) 1283 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve: (s)', st, totst 1284 CALL Info( 'HeatSolve', Message, Level=4 ) 1285!------------------------------------------------------------------------------ 1286! If modelling phase change (and if requested by the user), check if any 1287! node has jumped over the phase change interval, and if so, reduce 1288! timestep and or relaxation and recompute. 1289!------------------------------------------------------------------------------ 1290 IF (PhaseChange .AND. CheckLatentHeatRelease ) THEN 1291!------------------------------------------------------------------------------ 1292 IF ( CheckLatentHeat() ) THEN 1293 Enthalpy_h(1:LocalNodes) = PrevSolution 1294 Norm = PrevNorm 1295 1296 IF ( TransientSimulation ) THEN 1297 dt = dt / 2 1298 Solver % dt = dt 1299 WRITE( Message, * ) & 1300 'Latent heat release check: reducing timestep to: ',dt 1301 CALL Info( 'HeatSolve', Message, Level=4 ) 1302 ELSE 1303 Relax = Relax / 2 1304 CALL ListAddConstReal( Solver % Values, & 1305 'Nonlinear System Relaxation Factor', Relax ) 1306 WRITE( Message, * ) & 1307 'Latent heat release check: reducing relaxation to: ',Relax 1308 CALL Info( 'HeatSolve', Message, Level=4 ) 1309 END IF 1310 1311 CYCLE 1312 END IF 1313 IF ( .NOT.TransientSimulation ) PrevSolution=Enthalpy_h(1:LocalNodes) 1314 END IF 1315!------------------------------------------------------------------------------ 1316 1317 RelativeChange = Solver % Variable % NonlinChange 1318 1319 WRITE( Message, * ) 'Result Norm : ',Norm 1320 CALL Info( 'HeatSolve', Message, Level=4 ) 1321 WRITE( Message, * ) 'Relative Change : ',RelativeChange 1322 CALL Info( 'HeatSolve', Message, Level=4 ) 1323 1324 IF ( RelativeChange < NewtonTol .OR. iter >= NewtonIter ) & 1325 NewtonLinearization = .TRUE. 1326 IF ( RelativeChange < NonlinearTol .AND. & 1327 (.NOT. SmartHeaterControl .OR. SmartTolReached)) EXIT 1328 1329 IF(SmartHeaterControl) THEN 1330 IF ( RelativeChange < SmartTol ) THEN 1331 SmartTolReached = .TRUE. 1332 YY = Enthalpy_h 1333 END IF 1334 END IF 1335 1336!------------------------------------------------------------------------------ 1337 END DO ! of the nonlinear iteration 1338!------------------------------------------------------------------------------ 1339 IF(TransientHeaterControl) THEN 1340 PowerRelax = GetCReal(Solver % Values,'Smart Heater Relaxation Factor', GotIt) 1341 IF(.NOT. GotIt) PowerRelax = 1.0_dp 1342 PowerSensitivity = ListGetConstReal(Solver % Values,'Smart Heater Power Sensivity',GotIt,UnFoundFatal=UnFoundFatal) 1343 !Previous default value: PowerSensitivity = 4.0_dp 1344 PowerScaling = PowerScaling * (1 + PowerSensitivity * PowerRelax * (MeltPoint/yave - 1.0d0) ) 1345 1346 IF( ListGetLogical( Solver % Values,'Smart Heater Transient Speedup',GotIt ) ) THEN 1347 Enthalpy_h = Enthalpy_h * (1 + PowerRelax * (MeltPoint/yave - 1.0d0) ) 1348 END IF 1349 YY = Enthalpy_h 1350 END IF 1351 1352 1353!------------------------------------------------------------------------------ 1354! Compute temperature homologous and water content 1355!------------------------------------------------------------------------------ 1356 1357 WaterName = GetString(SolverParams , 'Exported Variable 2', Found ) 1358 1359 IF (.NOT.Found) & 1360 CALL FATAL('EnthalpySolver','No value > Exported Variable 2 (water content) < found in Solver') 1361 1362 TempName = GetString(SolverParams , 'Exported Variable 3', Found ) 1363 1364 IF (.NOT.Found) & 1365 CALL FATAL('EnthalpySolver','No value > Exported Variable 3 (temperature) < found in Solver') 1366 1367 WaterVar => VariableGet(Solver % Mesh % Variables, WaterName) 1368 TemphomoVar => VariableGet(Solver % Mesh % Variables,TempName) 1369 1370 L_heat = GetConstReal(Model % Constants, "L_heat",GotIt) 1371 IF (.NOT.GotIt) THEN 1372 CALL WARN('EnthalpySolver', 'No Keyword >L_heat< defined in model constants. Using >334000Jkg-1< as default.') 1373 L_heat = 334000.0 1374 END IF 1375 1376do i=1,Model % NumberOfNodes 1377 1378hi = Enthalpy_h (EnthalpyPerm (i) ) 1379hm = PhaseChangeEnthValues (PhaseChangeEnthPerm(i)) 1380 1381 if (hi<hm) then 1382 WaterVar % values ( WaterVar % perm (i) ) = 0.0 1383 TemphomoVar % values ( TemphomoVar % perm (i) ) = & 1384 & (-B_cap+(B_cap**2+A_cap*2*(A_cap*0.5*Tref**2+B_cap*Tref+hi))**0.5 ) / A_cap - 273.16 ! Use CP(T)=B_cap+A_cap*T 1385 else 1386 TemphomoVar % values ( TemphomoVar % perm (i) ) = & 1387 & (-B_cap+(B_cap**2+A_cap*2*(A_cap*0.5*Tref**2+B_cap*Tref+hm))**0.5 ) / A_cap - 273.16 ! Use CP(T)=B_cap+A_cap*T 1388 WaterVar % values ( WaterVar % perm (i) ) = (hi-hm)/L_heat 1389 endif 1390 1391 1392 1393enddo 1394 1395!------------------------------------------------------------------------------ 1396! Compute cumulative time done by now and time remaining 1397!------------------------------------------------------------------------------ 1398 IF ( .NOT. TransientSimulation ) EXIT 1399 CumulativeTime = CumulativeTime + dt 1400 dt = Timestep - CumulativeTime 1401 1402 END DO ! time interval 1403 Solver % dt = Timestep 1404 1405!------------------------------------------------------------------------------ 1406 CALL ListAddConstReal( Solver % Values, & 1407 'Nonlinear System Relaxation Factor', SaveRelax ) 1408!------------------------------------------------------------------------------ 1409 1410 DEALLOCATE( PrevSolution ) 1411 1412 IF ( ListGetLogical( Solver % Values, 'Adaptive Mesh Refinement', Found ) ) & 1413 CALL RefineMesh( Model,Solver,Enthalpy_h,EnthalpyPerm, & 1414 HeatInsideResidual, HeatEdgeResidual, HeatBoundaryResidual ) 1415 1416CONTAINS 1417 1418 1419 1420!------------------------------------------------------------------------------ 1421 SUBROUTINE AddHeatFluxBC() 1422!------------------------------------------------------------------------------ 1423 CALL GetElementNodes( ElementNodes ) 1424 1425 HeatTransferCoeff = 0.0D0 1426 LOAD = 0.0D0 1427!------------------------------------------------------------------------------ 1428! BC: -k@T/@n = \epsilon\sigma(T^4 - Text^4) 1429!------------------------------------------------------------------------------ 1430 RadiationFlag = GetString( BC, 'Radiation', Found ) 1431 1432 IF ( Found .AND. RadiationFlag(1:4) /= 'none' ) THEN 1433 NodalEmissivity(1:n) = GetReal(BC, 'Emissivity', Found) 1434 IF(.NOT. Found) THEN 1435 NodalEmissivity(1:n) = GetParentMatProp( 'Emissivity' ) 1436 END IF 1437 Emissivity = SUM( NodalEmissivity(1:n) ) / n 1438 1439!------------------------------------------------------------------------------ 1440 IF ( RadiationFlag(1:9) == 'idealized' ) THEN 1441 AText(1:n) = GetReal( BC, 'Radiation External Enthalpy_h',Found ) 1442 IF(.NOT. Found) AText(1:n) = GetReal( BC, 'External Enthalpy_h' ) 1443 ELSE 1444 CALL DiffuseGrayRadiation( Model, Solver, Element, & 1445 Enthalpy_h, EnthalpyPerm, ForceVector, VisibleFraction, Text) 1446 1447 IF( GetLogical( BC, 'Radiation Boundary Open', Found) ) THEN 1448 AText(1:n) = GetReal( BC, 'Radiation External Enthalpy_h',Found ) 1449 IF(.NOT. Found) AText(1:n) = GetReal( BC, 'External Enthalpy_h' ) 1450 IF( VisibleFraction >= 1.0_dp ) THEN 1451 Atext(1:n) = Text 1452 ELSE 1453 Atext(1:n) = ( (1 - VisibleFraction) * Atext(1:n)**4 + & 1454 VisibleFraction * Text**4 ) ** 0.25_dp 1455 END IF 1456 ELSE 1457 AText(1:n) = Text 1458 END IF 1459 END IF 1460!------------------------------------------------------------------------------ 1461! Add our own contribution to surface Enthalpy_h (and external 1462! if using linear type iteration or idealized radiation) 1463!------------------------------------------------------------------------------ 1464 DO j=1,n 1465 k = EnthalpyPerm(Element % NodeIndexes(j)) 1466 Text = AText(j) 1467 1468 IF ( .NOT. HeatGapBC .AND. NewtonLinearization ) THEN 1469 HeatTransferCoeff(j) = Emissivity * 4*Enthalpy_h(k)**3 * & 1470 StefanBoltzmann 1471 LOAD(j) = Emissivity*(3*Enthalpy_h(k)**4+Text**4) * & 1472 StefanBoltzmann 1473 ELSE 1474 HeatTransferCoeff(j) = Emissivity * (Enthalpy_h(k)**3 + & 1475 Enthalpy_h(k)**2*Text+Enthalpy_h(k)*Text**2 + Text**3) * & 1476 StefanBoltzmann 1477 LOAD(j) = HeatTransferCoeff(j) * Text 1478 END IF 1479 END DO 1480 END IF ! of radition 1481!------------------------------------------------------------------------------ 1482 1483 Work(1:n) = GetReal( BC, 'Heat Transfer Coefficient',Found ) 1484 IF ( Found ) THEN 1485 AText(1:n) = GetReal( BC, 'External Enthalpy_h',Found ) 1486 DO j=1,n 1487!------------------------------------------------------------------------------ 1488! BC: -k@T/@n = \alpha(T - Text) 1489!------------------------------------------------------------------------------ 1490 k = EnthalpyPerm(Element % NodeIndexes(j)) 1491 LOAD(j) = LOAD(j) + Work(j) * AText(j) 1492 HeatTransferCoeff(j) = HeatTransferCoeff(j) + Work(j) 1493 END DO 1494 END IF 1495 1496!------------------------------------------------------------------------------ 1497! BC: -k@T/@n = (rho*L)*v.n 1498! Heating related to pulling is possible only in ss cases where pull velocity 1499! is desrcibed. 1500!------------------------------------------------------------------------------ 1501 1502 IF( GetLogical( BC, 'Phase Change',Found ) ) THEN 1503 PhaseVelocity(1,1:n) = GetReal( BC,'Phase Velocity 1', Found ) 1504 PhaseVelocity(2,1:n) = GetReal( BC,'Phase Velocity 2', Found ) 1505 PhaseVelocity(3,1:n) = GetReal( BC,'Phase Velocity 3', Found ) 1506 1507 ! Ensure that the latent heat and density come from the same side 1508 LatentHeat(1:n) = GetParentMatProp( 'Latent Heat', & 1509 UElement = Element, UParent = Parent ) 1510 IF(.NOT. ASSOCIATED(Parent) ) THEN 1511 CALL Warn('HeatSolve','Parent not associated') 1512 ELSE 1513 k = GetInteger(Model % Bodies(Parent % BodyId) % Values,'Material') 1514 Density(1:n) = GetReal( Model % Materials(k) % Values, 'Enthalpy Density' ) 1515 END IF 1516 1517 ! This could be rather put as a new type of BC into the assembly routine and 1518 ! then the Normal could be taken at the proper Gaussian integration points. 1519 Normal = NormalVector( Element, ElementNodes, 0.0_dp, 0.0_dp, .TRUE. ) 1520 1521 DO i=1,n 1522 LOAD(i) = LOAD(i) + & 1523 LatentHeat(i) * Density(i) * SUM( Normal(1:3) * PhaseVelocity(1:3,i)) 1524 END DO 1525 END IF 1526 1527!------------------------------------------------------------------------------ 1528! BC: -k@T/@n = g 1529!------------------------------------------------------------------------------ 1530 LOAD(1:n) = LOAD(1:n) + GetReal( BC, 'Enthalpy Heat Flux', Found ) 1531 WaterDiffusivity(1:n) = GetReal( Material, 'Enthalpy Water Diffusivity', Found ) 1532 1533 InfBC = ListGetLogical( BC,'Infinity BC '//TRIM(VarName),GotIt) 1534 IF( InfBC ) THEN 1535 AText(1:n) = GetReal( BC,'Infinity BC '//TRIM(VarName)//' Offset',GotIt) 1536 ! currently only isotropic heat conductivity supported 1537 HeatConductivityIso(1:n) = GetParentMatProp('Enthalpy Heat Diffusivity',Element,GotIt) 1538 1539 DO k=1,n 1540 IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(k)))<& 1541 &Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(k)))) then 1542 HeatConductivityIso(k) = WaterDiffusivity(k) 1543 ENDIF 1544 ENDDO 1545 1546 1547 IF(.NOT. GotIt) THEN 1548 CALL Fatal( 'HeatSolver','Could not find > Heat Conductivity < for parent!' ) 1549 END IF 1550 END IF 1551 1552!------------------------------------------------------------------------------ 1553! Get element matrix and rhs due to boundary conditions ... 1554!------------------------------------------------------------------------------ 1555 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1556 CALL DiffuseConvectiveBoundary( STIFF,FORCE, & 1557 LOAD,HeatTransferCoeff,InfBC,HeatConductivityIso,AText(1:n),& 1558 Element,n,ElementNodes ) 1559 ELSE 1560 IF( InfBC ) THEN 1561 CALL Fatal('HeatSolver','Infinity BC not implemented only for cartersian case!') 1562 END IF 1563 CALL DiffuseConvectiveGenBoundary(STIFF,FORCE,& 1564 LOAD,HeatTransferCoeff,Element,n,ElementNodes ) 1565 END IF 1566 1567!------------------------------------------------------------------------------ 1568! Update global matrices from local matrices 1569!------------------------------------------------------------------------------ 1570 IF ( TransientAssembly .AND. .NOT. ConstantBulk ) THEN 1571 MASS = 0.d0 1572 CALL Default1stOrderTime( MASS, STIFF, FORCE ) 1573 END IF 1574 1575 IF ( HeatGapBC ) & 1576 CALL AddHeatGap( Solver, Element, STIFF, EnthalpyPerm) 1577 1578 CALL DefaultUpdateEquations( STIFF, FORCE ) 1579!------------------------------------------------------------------------------ 1580 END SUBROUTINE AddHeatFluxBC 1581!------------------------------------------------------------------------------ 1582 1583 1584!------------------------------------------------------------------------------ 1585 SUBROUTINE AddGlobalTime() 1586!------------------------------------------------------------------------------ 1587 INTEGER :: i,j,k,l,n 1588 REAL(KIND=dp) :: FORCE(1) 1589 REAL(KIND=dp), POINTER :: SaveValues(:) => NULL() 1590 SAVE STIFF, MASS, X 1591 REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:),MASS(:,:), X(:,:) 1592 1593 IF ( .NOT.ASSOCIATED(Solver % Variable % Values, SaveValues) ) THEN 1594 IF ( ALLOCATED(STIFF) ) DEALLOCATE( STIFF,MASS,X ) 1595 n = 0 1596 DO i=1,Solver % Matrix % NumberOfRows 1597 n = MAX( n,Solver % Matrix % Rows(i+1)-Solver % Matrix % Rows(i) ) 1598 END DO 1599 k = SIZE(Solver % Variable % PrevValues,2) 1600 ALLOCATE( STIFF(1,n),MASS(1,n),X(n,k) ) 1601 1602 SaveValues => Solver % Variable % Values 1603 END IF 1604 1605 DO i=1,Solver % Matrix % NumberOFRows 1606 n = 0 1607 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 1608 n=n+1 1609 STIFF(1,n) = Solver % Matrix % Values(j) 1610 MASS(1,n) = Solver % Matrix % MassValues(j) 1611 X(n,:) = Solver % Variable % PrevValues(Solver % Matrix % Cols(j),:) 1612 END DO 1613 FORCE(1) = Solver % Matrix % RHS(i) 1614 Solver % Matrix % Force(i,1) = FORCE(1) 1615 k = MIN( Solver % DoneTime, Solver % Order ) 1616 CALL BDFLocal( n, dt, MASS, STIFF, FORCE, X, k ) 1617 1618 n = 0 1619 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 1620 n=n+1 1621 Solver % Matrix % Values(j) = STIFF(1,n) 1622 END DO 1623 Solver % Matrix % RHS(i) = FORCE(1) 1624 END DO 1625!------------------------------------------------------------------------------ 1626 END SUBROUTINE AddGlobalTime 1627!------------------------------------------------------------------------------ 1628 1629!------------------------------------------------------------------------------ 1630 SUBROUTINE DiffuseGrayRadiation( Model, Solver, Element, & 1631 Enthalpy_h, EnthalpyPerm, ForceVector,AngleFraction, Text) 1632!------------------------------------------------------------------------------ 1633 TYPE(Model_t) :: Model 1634 TYPE(Solver_t) :: Solver 1635 TYPE(Element_t), POINTER :: Element 1636 INTEGER :: EnthalpyPerm(:) 1637 REAL(KIND=dp) :: Enthalpy_h(:), ForceVector(:) 1638 REAL(KIND=dp) :: AngleFraction, Text 1639!------------------------------------------------------------------------------ 1640 REAL(KIND=dp) :: Area, Asum 1641 INTEGER :: i,j,k,l,m,ImplicitFactors 1642 INTEGER, POINTER :: ElementList(:) 1643!------------------------------------------------------------------------------ 1644! If linear iteration compute radiation load 1645!------------------------------------------------------------------------------ 1646 1647 1648 Asum = 0.0d0 1649 IF ( .NOT. NewtonLinearization ) THEN 1650 Text = ComputeRadiationLoad( Model, Solver % Mesh, Element, & 1651 Enthalpy_h, EnthalpyPerm, Emissivity, AngleFraction) 1652 ELSE ! Full Newton-Raphson solver 1653!------------------------------------------------------------------------------ 1654! Go trough surfaces (j) this surface (i) is getting 1655! radiated from. 1656!------------------------------------------------------------------------------ 1657 1658 Area = ElementArea( Solver % Mesh, Element, n ) 1659 ElementList => Element % BoundaryInfo % GebhardtFactors % Elements 1660 1661 DO j=1,Element % BoundaryInfo % GebhardtFactors % NumberOfFactors 1662 1663 RadiationElement => Solver % Mesh % Elements( ElementList(j) ) 1664 1665 Text = ComputeRadiationCoeff(Model,Solver % Mesh,Element,j) / ( Area ) 1666 Asum = Asum + Text 1667!------------------------------------------------------------------------------ 1668! Gebhardt factors are given elementwise at the center 1669! of the element, so take average of nodal Enthalpy_hs 1670! (or integrate over surface j) 1671!------------------------------------------------------------------------------ 1672 1673 k = RadiationElement % TYPE % NumberOfNodes 1674 ImplicitFactors = Element % BoundaryInfo % GebhardtFactors % NumberOfImplicitFactors 1675 IF(ImplicitFactors == 0) & 1676 ImplicitFactors = Element % BoundaryInfo % GebhardtFactors % NumberOfFactors 1677 1678 IF(j <= ImplicitFactors) THEN 1679 1680 S = (SUM( Enthalpy_h( EnthalpyPerm( RadiationElement % & 1681 NodeIndexes))**4 )/k )**(1.0d0/4.0d0) 1682!------------------------------------------------------------------------------ 1683! Linearization of the G_jiT^4_j term 1684!------------------------------------------------------------------------------ 1685 HeatTransferCoeff(1:n) = -4 * Text * S**3 * StefanBoltzmann 1686 LOAD(1:n) = -3 * Text * S**4 * StefanBoltzmann 1687!------------------------------------------------------------------------------ 1688! Integrate the contribution of surface j over surface i 1689! and add to global matrix 1690!------------------------------------------------------------------------------ 1691 CALL IntegOverA( STIFF, FORCE, LOAD, & 1692 HeatTransferCoeff, Element, n, k, ElementNodes ) 1693 1694 IF ( TransientAssembly ) THEN 1695 MASS = 0.d0 1696 CALL Add1stOrderTime( MASS, STIFF, & 1697 FORCE,dt,n,1,EnthalpyPerm(Element % NodeIndexes),Solver ) 1698 END IF 1699 1700 DO m=1,n 1701 k1 = EnthalpyPerm( Element % NodeIndexes(m) ) 1702 DO l=1,k 1703 k2 = EnthalpyPerm( RadiationElement % NodeIndexes(l) ) 1704 CALL AddToMatrixElement( StiffMatrix,k1, & 1705 k2,STIFF(m,l) ) 1706 END DO 1707 ForceVector(k1) = ForceVector(k1) + FORCE(m) 1708 END DO 1709 1710 ELSE 1711 1712 S = (SUM( Enthalpy_h( EnthalpyPerm( RadiationElement % & 1713 NodeIndexes))**4 )/k ) 1714 1715 HeatTransferCoeff(1:n) = 0.0d0 1716 LOAD(1:n) = Text * S * StefanBoltzmann 1717 1718 CALL IntegOverA( STIFF, FORCE, LOAD, & 1719 HeatTransferCoeff, Element, n, k, ElementNodes ) 1720 1721 DO m=1,n 1722 k1 = EnthalpyPerm( Element % NodeIndexes(m) ) 1723 ForceVector(k1) = ForceVector(k1) + FORCE(m) 1724 END DO 1725 1726 END IF 1727 1728 END DO 1729 1730!------------------------------------------------------------------------------ 1731! We have already added all external Enthalpy_h contributions 1732! to the matrix for the Newton type iteration 1733!------------------------------------------------------------------------------ 1734 AngleFraction = Asum / Emissivity 1735 Text = 0.0 1736 1737 END IF ! of newton-raphson 1738 1739 END SUBROUTINE DiffuseGrayRadiation 1740!------------------------------------------------------------------------------ 1741 1742 1743!------------------------------------------------------------------------------ 1744 SUBROUTINE EffectiveHeatCapacity() 1745 LOGICAL :: Found, Specific 1746 REAL(KIND=dp), ALLOCATABLE :: dT(:) 1747 1748!------------------------------------------------------------------------------ 1749! See if Enthalpy_h gradient indside the element is large enough 1750! to use the c_p = SQRT( (dH/dx)^2 / (dT/dx)^2 ), otherwise 1751! use c_p = dH/dT, or if in time dependent simulation, use 1752! c_p = (dH/dt) / (dT/dt), if requested. 1753!------------------------------------------------------------------------------ 1754 1755 SELECT CASE(PhaseModel) 1756!------------------------------------------------------------------------------ 1757 CASE( 'spatial 1' ) 1758 PhaseChangeModel = PHASE_SPATIAL_1 1759!------------------------------------------------------------------------------ 1760 1761 CASE( 'spatial 2' ) 1762!------------------------------------------------------------------------------ 1763! Check if local variation of Enthalpy_h is large enough to actually use the 1764! Spatial 2 model. Should perhaps be scaled to element size (or actually 1765! compute the gradient, but this will do for now...). 1766!------------------------------------------------------------------------------ 1767 s = MAXVAL(LocalEnthalpy_h(1:n))-MINVAL(LocalEnthalpy_h(1:n)) 1768 IF ( s < AEPS ) THEN 1769 PhaseChangeModel = PHASE_SPATIAL_1 1770 ELSE 1771 PhaseChangeModel = PHASE_SPATIAL_2 1772 END IF 1773 1774!------------------------------------------------------------------------------ 1775! Note that here HeatCapacity is miused for saving dT. 1776!------------------------------------------------------------------------------ 1777 CASE('temporal') 1778 IF ( TransientSimulation ) THEN 1779 ALLOCATE( dT(n) ) 1780 dT(1:n) = Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) - & 1781 PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes)) 1782 1783 IF ( ANY(ABS(dT(1:n)) < AEPS) ) THEN 1784 PhaseChangeModel = PHASE_SPATIAL_1 1785 ELSE 1786 PhaseChangeModel = PHASE_TEMPORAL 1787 END IF 1788 ELSE 1789 PhaseChangeModel = PHASE_SPATIAL_1 1790 END IF 1791 1792!------------------------------------------------------------------------------ 1793 CASE DEFAULT 1794 PhaseChangeModel = PHASE_SPATIAL_1 1795 1796 END SELECT 1797!------------------------------------------------------------------------------ 1798 1799 PhaseSpatial = ( PhaseChangeModel == PHASE_SPATIAL_2 ) 1800 Specific = ListCheckPresent( Material,'Specific Enthalpy') 1801 1802!----------------------------------------------------------------------------- 1803 SELECT CASE( PhaseChangeModel ) 1804 1805!------------------------------------------------------------------------------ 1806! This phase change model is available only for some type of real entries 1807! that have an implemented analytical derivation rule. 1808!----------------------------------------------------------------------------- 1809 CASE( PHASE_SPATIAL_1 ) 1810 HeatCapacity(1:n) = ListGetReal( Material, & 1811 'Effective Heat Capacity', n,Element % NodeIndexes, Found ) 1812 IF ( .NOT. Found ) THEN 1813 IF( Specific ) THEN 1814 HeatCapacity(1:n) = ListGetDerivValue( Material, & 1815 'Specific Enthalpy', n,Element % NodeIndexes ) 1816 HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n) 1817 ELSE 1818 HeatCapacity(1:n) = ListGetDerivValue( Material, & 1819 'Enthalpy', n,Element % NodeIndexes ) 1820 END IF 1821 END IF 1822 1823!--------------------------------------------------------------------------------------- 1824! Note that for the 'spatial 2' model the evaluation of c_p is done in each integration 1825! point and thus Enthalphy and PhaseSpatial flag are used instead of HeatCapacity directly. 1826!----------------------------------------------------------------------------------------- 1827 CASE( PHASE_SPATIAL_2 ) 1828 IF( Specific ) THEN 1829 Enthalpy(1:n) = ListGetReal(Material,'Specific Enthalpy',n,Element % NodeIndexes) 1830 Enthalpy(1:n) = Density(1:n) * Enthalpy(1:n) 1831 ELSE 1832 Enthalpy(1:n) = ListGetReal(Material,'Enthalpy',n,Element % NodeIndexes) 1833 END IF 1834 1835!------------------------------------------------------------------------------ 1836 CASE( PHASE_TEMPORAL ) 1837 ! When retrieving the value of enthalphy on the previous timestep 1838 ! the relevant entries of the Enthalpy_h solution in the global vector 1839 ! are tampered in order to make the ListGetReal command work as wanted. 1840 ! 1) Values at current Enthalpy_h 1841 !------------------------------------------------------------------------ 1842 IF( Specific ) THEN 1843 Work(1:n) = ListGetReal( Material,'Specific Enthalpy',n,Element % NodeIndexes ) 1844 ELSE 1845 Work(1:n) = ListGetReal( Material,'Enthalpy',n,Element % NodeIndexes ) 1846 END IF 1847 1848 ! 2) Values at previous Enthalpy_h 1849 Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) = & 1850 PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes)) 1851 1852 IF( Specific ) THEN 1853 Work(1:n) = Work(1:n) - ListGetReal( Material,'Specific Enthalpy', & 1854 n,Element % NodeIndexes ) 1855 HeatCapacity(1:n) = Density(1:n) * Work(1:n) / dT(1:n) 1856 ELSE 1857 Work(1:n) = Work(1:n) - ListGetReal( Material,'Enthalpy', & 1858 n,Element % NodeIndexes ) 1859 HeatCapacity(1:n) = Work(1:n) / dT(1:n) 1860 END IF 1861 1862 ! Revert to current Enthalpy_h 1863 Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) = & 1864 PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes)) + dT(1:n) 1865 1866!------------------------------------------------------------------------------ 1867 END SELECT 1868!------------------------------------------------------------------------------ 1869 END SUBROUTINE EffectiveHeatCapacity 1870!------------------------------------------------------------------------------ 1871 1872 1873 1874!------------------------------------------------------------------------------ 1875 FUNCTION CheckLatentHeat() RESULT(Failure) 1876!------------------------------------------------------------------------------ 1877 LOGICAL :: Failure, PhaseChange, CheckLatentHeatRelease 1878 INTEGER :: t, eq_id, body_id 1879 CHARACTER(LEN=MAX_NAME_LEN) :: PhaseModel 1880!------------------------------------------------------------------------------ 1881 1882 Failure = .FALSE. 1883!------------------------------------------------------------------------------ 1884 DO t=1,Solver % Mesh % NumberOfBulkElements 1885!------------------------------------------------------------------------------ 1886! Check if this element belongs to a body where Enthalpy_h 1887! has been calculated 1888!------------------------------------------------------------------------------ 1889 Element => Solver % Mesh % Elements(t) 1890 1891 NodeIndexes => Element % NodeIndexes 1892 IF ( ANY( EnthalpyPerm( NodeIndexes ) <= 0 ) ) CYCLE 1893 1894 body_id = Element % Bodyid 1895 eq_id = ListGetInteger( Model % Bodies(body_id) % Values, & 1896 'Equation', minv=1, maxv=Model % NumberOfEquations ) 1897 1898 PhaseModel = ListGetString( Model % Equations(eq_id) % Values, & 1899 'Phase Change Model',Found ) 1900 1901 PhaseChange = Found .AND. (PhaseModel(1:4) /= 'none') 1902 1903 IF ( PhaseChange ) THEN 1904 CheckLatentHeatRelease = ListGetLogical(Model % Equations(eq_id) % & 1905 Values, 'Check Latent Heat Release',Found ) 1906 END IF 1907 IF ( .NOT. ( PhaseChange .AND. CheckLatentHeatRelease ) ) CYCLE 1908 1909 n = Element % TYPE % NumberOfNodes 1910!------------------------------------------------------------------------------ 1911! Set the current element pointer in the model structure to 1912! reflect the element being processed 1913!------------------------------------------------------------------------------ 1914 Model % CurrentElement => Element 1915!------------------------------------------------------------------------------ 1916!------------------------------------------------------------------------------ 1917! Get element material parameters 1918!------------------------------------------------------------------------------ 1919 k = ListGetInteger( Model % Bodies(body_id) % Values,'Material', & 1920 minv=1, maxv=Model % NumberOfMaterials ) 1921 Material => Model % Materials(k) % Values 1922 1923 PhaseChangeIntervals => ListGetConstRealArray( Material, & 1924 'Phase Change Intervals' ) 1925 1926 DO k=1,n 1927 i = EnthalpyPerm( NodeIndexes(k) ) 1928 DO j=1,SIZE(PhaseChangeIntervals,2) 1929 IF ( ( Enthalpy_h(i) < PhaseChangeIntervals(1,j) .AND. & 1930 PrevSolution(i) > PhaseChangeIntervals(2,j) ).OR. & 1931 ( Enthalpy_h(i) > PhaseChangeIntervals(2,j) .AND. & 1932 PrevSolution(i) < PhaseChangeIntervals(1,j) ) ) THEN 1933 Failure = .TRUE. 1934 EXIT 1935 END IF 1936 END DO 1937 IF ( Failure ) EXIT 1938 END DO 1939 IF ( Failure ) EXIT 1940 END DO 1941!------------------------------------------------------------------------------ 1942 END FUNCTION CheckLatentHeat 1943!------------------------------------------------------------------------------ 1944 1945 1946 1947!------------------------------------------------------------------------------ 1948 SUBROUTINE IntegOverA( BoundaryMatrix, BoundaryVector, & 1949 LOAD, NodalAlpha, Element, n, m, Nodes ) 1950!------------------------------------------------------------------------------ 1951 REAL(KIND=dp) :: BoundaryMatrix(:,:),BoundaryVector(:), & 1952 LOAD(:),NodalAlpha(:) 1953 1954 TYPE(Nodes_t) :: Nodes 1955 TYPE(Element_t) :: Element 1956 1957 INTEGER :: n, m 1958 1959 REAL(KIND=dp) :: Basis(n) 1960 REAL(KIND=dp) :: dBasisdx(n,3),SqrtElementMetric 1961 1962 REAL(KIND=dp) :: u,v,w,s,x,y,z 1963 REAL(KIND=dp) :: Force,Alpha 1964 REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:) 1965 1966 INTEGER :: i,t,q,p,N_Integ 1967 1968 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1969 1970 LOGICAL :: stat 1971!------------------------------------------------------------------------------ 1972 1973 BoundaryVector = 0.0D0 1974 BoundaryMatrix = 0.0D0 1975!------------------------------------------------------------------------------ 1976! Integration stuff 1977!------------------------------------------------------------------------------ 1978 IntegStuff = GaussPoints( Element ) 1979 U_Integ => IntegStuff % u 1980 V_Integ => IntegStuff % v 1981 W_Integ => IntegStuff % w 1982 S_Integ => IntegStuff % s 1983 N_Integ = IntegStuff % n 1984 1985!------------------------------------------------------------------------------ 1986! Now we start integrating 1987!------------------------------------------------------------------------------ 1988 DO t=1,N_Integ 1989 u = U_Integ(t) 1990 v = V_Integ(t) 1991 w = W_Integ(t) 1992!------------------------------------------------------------------------------ 1993! Basis function values & derivatives at the integration point 1994!------------------------------------------------------------------------------ 1995 stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, & 1996 Basis,dBasisdx ) 1997 1998 s = SqrtElementMetric * S_Integ(t) 1999!------------------------------------------------------------------------------ 2000! Coordinatesystem dependent info 2001!------------------------------------------------------------------------------ 2002 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2003 x = SUM( Nodes % x(1:n)*Basis ) 2004 y = SUM( Nodes % y(1:n)*Basis ) 2005 z = SUM( Nodes % z(1:n)*Basis ) 2006 s = s * CoordinateSqrtMetric( x,y,z ) 2007 END IF 2008!------------------------------------------------------------------------------ 2009 Force = SUM( LOAD(1:n) * Basis ) 2010 Alpha = SUM( NodalAlpha(1:n) * Basis ) 2011 2012 DO p=1,N 2013 DO q=1,M 2014 BoundaryMatrix(p,q) = BoundaryMatrix(p,q) + & 2015 s * Alpha * Basis(p) / m 2016 END DO 2017 END DO 2018 2019 DO p=1,N 2020 BoundaryVector(p) = BoundaryVector(p) + s * Force * Basis(p) 2021 END DO 2022 END DO 2023 END SUBROUTINE IntegOverA 2024!------------------------------------------------------------------------------ 2025 2026 2027 2028!------------------------------------------------------------------------------ 2029 SUBROUTINE FindGapIndexes( Element, Indexes, n ) 2030!------------------------------------------------------------------------------ 2031 TYPE(Element_t) :: Element 2032 INTEGER :: n,Indexes(:) 2033!------------------------------------------------------------------------------ 2034 TYPE(Element_t), POINTER :: Parent,Left,Right 2035 INTEGER :: i,j,k,l 2036 REAL(KIND=dp) :: x0,y0,z0,x,y,z 2037!------------------------------------------------------------------------------ 2038 Left => Element % BoundaryInfo % Left 2039 Right => Element % BoundaryInfo % Right 2040 2041 IF ( .NOT.ASSOCIATED(Left) .OR. .NOT.ASSOCIATED(Right) ) RETURN 2042 2043 l = 0 2044 DO i=1,n 2045 Parent => Left 2046 k = Element % NodeIndexes(i) 2047 2048 IF ( ANY( Parent % NodeIndexes == k ) ) & 2049 Parent => Right 2050 2051 x0 = ElementNodes % x(i) 2052 y0 = ElementNodes % y(i) 2053 z0 = ElementNodes % z(i) 2054 DO j=1,Parent % TYPE % NumberOfNodes 2055 k = Parent % NodeIndexes(j) 2056 x = Solver % Mesh % Nodes % x(k) - x0 2057 y = Solver % Mesh % Nodes % y(k) - y0 2058 z = Solver % Mesh % Nodes % z(k) - z0 2059 IF ( x**2 + y**2 + z**2 < AEPS ) EXIT 2060 END DO 2061 Indexes(i) = k 2062 END DO 2063!------------------------------------------------------------------------------ 2064 END SUBROUTINE FindGapIndexes 2065!------------------------------------------------------------------------------ 2066 2067 2068!------------------------------------------------------------------------------ 2069 SUBROUTINE AddHeatGap( Solver, Element, STIFF, EnthalpyPerm ) 2070!------------------------------------------------------------------------------ 2071 TYPE(Solver_t) :: Solver 2072 REAL(KIND=dp) :: STIFF(:,:) 2073 INTEGER :: EnthalpyPerm(:) 2074 TYPE(Element_t) :: Element 2075!------------------------------------------------------------------------------ 2076 TYPE(Element_t), POINTER :: Parent,Left,Right 2077 INTEGER :: i,j,k,l, Ind(n) 2078 REAL(KIND=dp) :: x0,y0,z0,x,y,z 2079!------------------------------------------------------------------------------ 2080 CALL FindGapIndexes( Element, Ind, n ) 2081 DO i=1,n 2082 DO j=1,n 2083 k = EnthalpyPerm( Element % NodeIndexes(i) ) 2084 l = EnthalpyPerm( Ind(j) ) 2085 IF ( k > 0 .AND. l > 0 ) THEN 2086 CALL AddToMatrixElement( Solver % Matrix,k,l,-STIFF(i,j) ) 2087 END IF 2088 END DO 2089 END DO 2090!------------------------------------------------------------------------------ 2091 END SUBROUTINE AddHeatGap 2092!------------------------------------------------------------------------------ 2093 2094 2095!------------------------------------------------------------------------------ 2096 END SUBROUTINE EnthalpySolver 2097!------------------------------------------------------------------------------ 2098 2099 2100!------------------------------------------------------------------------------ 2101 FUNCTION HeatBoundaryResidual( Model, Edge, Mesh, Quant, Perm,Gnorm ) RESULT( Indicator ) 2102!------------------------------------------------------------------------------ 2103 USE DefUtils 2104 USE Radiation 2105 2106 IMPLICIT NONE 2107!------------------------------------------------------------------------------ 2108 TYPE(Model_t) :: Model 2109 INTEGER :: Perm(:) 2110 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 2111 TYPE( Mesh_t ), POINTER :: Mesh 2112 TYPE( Element_t ), POINTER :: Edge 2113!------------------------------------------------------------------------------ 2114 2115 TYPE(Nodes_t) :: Nodes, EdgeNodes 2116 TYPE(Element_t), POINTER :: Element, Bndry 2117 2118 INTEGER :: i,j,k,n,l,t,DIM,Pn,En 2119 LOGICAL :: stat, Found 2120 2121 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 2122 2123 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 2124 2125 REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:), ExtEnthalpy_h(:), & 2126 TransferCoeff(:), EdgeBasis(:), Basis(:), x(:), y(:), z(:), & 2127 dBasisdx(:,:), Enthalpy_h(:), Flux(:), NodalEmissivity(:) 2128 2129 REAL(KIND=dp) :: Conductivity, Emissivity, StefanBoltzmann 2130 2131 REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength, gx, gy, gz 2132 2133 REAL(KIND=dp) :: u, v, w, s, detJ 2134 2135 REAL(KIND=dp) :: Source, Residual, ResidualNorm, Area 2136 2137 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 2138 2139 LOGICAL :: First = .TRUE., Dirichlet 2140 SAVE Hwrk, First 2141!------------------------------------------------------------------------------ 2142 2143! Initialize: 2144! ----------- 2145 IF ( First ) THEN 2146 First = .FALSE. 2147 NULLIFY( Hwrk ) 2148 END IF 2149 2150 Indicator = 0.0d0 2151 Gnorm = 0.0d0 2152 2153 Metric = 0.0d0 2154 DO i=1,3 2155 Metric(i,i) = 1.0d0 2156 END DO 2157 2158 SELECT CASE( CurrentCoordinateSystem() ) 2159 CASE( AxisSymmetric, CylindricSymmetric ) 2160 dim = 3 2161 CASE DEFAULT 2162 dim = CoordinateSystemDimension() 2163 END SELECT 2164! 2165! --------------------------------------------- 2166 2167 Element => Edge % BoundaryInfo % Left 2168 2169 IF ( .NOT. ASSOCIATED( Element ) ) THEN 2170 Element => Edge % BoundaryInfo % Right 2171 ELSE IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) THEN 2172 Element => Edge % BoundaryInfo % Right 2173 END IF 2174 2175 IF ( .NOT. ASSOCIATED( Element ) ) RETURN 2176 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 2177 2178 En = Edge % TYPE % NumberOfNodes 2179 Pn = Element % TYPE % NumberOfNodes 2180 2181 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 2182 2183 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 2184 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 2185 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 2186 2187 ALLOCATE( Nodes % x(Pn), Nodes % y(Pn), Nodes % z(Pn) ) 2188 2189 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 2190 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 2191 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 2192 2193 ALLOCATE( Enthalpy_h(Pn), Basis(Pn), ExtEnthalpy_h(En), & 2194 TransferCoeff(En), x(En), y(En), z(En), EdgeBasis(En), & 2195 dBasisdx(Pn,3), NodalConductivity(En), Flux(En), & 2196 NodalEmissivity(En) ) 2197 2198 DO l = 1,En 2199 DO k = 1,Pn 2200 IF ( Edge % NodeIndexes(l) == Element % NodeIndexes(k) ) THEN 2201 x(l) = Element % TYPE % NodeU(k) 2202 y(l) = Element % TYPE % NodeV(k) 2203 z(l) = Element % TYPE % NodeW(k) 2204 EXIT 2205 END IF 2206 END DO 2207 END DO 2208! 2209! Integrate square of residual over boundary element: 2210! --------------------------------------------------- 2211 2212 Indicator = 0.0d0 2213 EdgeLength = 0.0d0 2214 ResidualNorm = 0.0d0 2215 2216 DO j=1,Model % NumberOfBCs 2217 IF ( Edge % BoundaryInfo % Constraint /= Model % BCs(j) % Tag ) CYCLE 2218 2219! IF ( .NOT. ListGetLogical( Model % BCs(j) % Values, & 2220! 'Enthalpy Heat Flux BC', Found ) ) CYCLE 2221 2222! 2223! Check if dirichlet BC given: 2224! ---------------------------- 2225 s = ListGetConstReal( Model % BCs(j) % Values,'Enthalpy_h',Dirichlet ) 2226 2227! Get various flux bc options: 2228! ---------------------------- 2229 2230! ...given flux: 2231! -------------- 2232 Flux(1:En) = ListGetReal( Model % BCs(j) % Values, & 2233 'Enthalpy Heat Flux', En, Edge % NodeIndexes, Found ) 2234 2235! ...convective heat transfer: 2236! ---------------------------- 2237 TransferCoeff(1:En) = ListGetReal( Model % BCs(j) % Values, & 2238 'Heat Transfer Coefficient', En, Edge % NodeIndexes, Found ) 2239 2240 ExtEnthalpy_h(1:En) = ListGetReal( Model % BCs(j) % Values, & 2241 'External Enthalpy_h', En, Edge % NodeIndexes, Found ) 2242 2243! ...black body radiation: 2244! ------------------------ 2245 Emissivity = 0.0d0 2246 StefanBoltzmann = 0.0d0 2247 2248 SELECT CASE(ListGetString(Model % BCs(j) % Values,'Radiation',Found)) 2249 !------------------ 2250 CASE( 'idealized' ) 2251 !------------------ 2252 2253 NodalEmissivity(1:En) = ListGetReal( Model % BCs(j) % Values, & 2254 'Emissivity', En, Edge % NodeIndexes, Found) 2255 IF(.NOT. Found) THEN 2256 NodalEmissivity(1:En) = GetParentMatProp( 'Emissivity', Edge) 2257 END IF 2258 Emissivity = SUM( NodalEmissivity(1:En)) / En 2259 2260 StefanBoltzMann = & 2261 ListGetConstReal( Model % Constants,'Stefan Boltzmann' ) 2262 2263 !--------------------- 2264 CASE( 'diffuse gray' ) 2265 !--------------------- 2266 2267 NodalEmissivity(1:En) = ListGetReal( Model % BCs(j) % Values, & 2268 'Emissivity', En, Edge % NodeIndexes, Found) 2269 IF(.NOT. Found) THEN 2270 NodalEmissivity(1:En) = GetParentMatProp( 'Emissivity', Edge) 2271 END IF 2272 Emissivity = SUM( NodalEmissivity(1:En)) / En 2273 2274 StefanBoltzMann = & 2275 ListGetConstReal( Model % Constants,'Stefan Boltzmann' ) 2276 2277 ExtEnthalpy_h(1:En) = ComputeRadiationLoad( Model, & 2278 Mesh, Edge, Quant, Perm, Emissivity ) 2279 END SELECT 2280 2281! get material parameters: 2282! ------------------------ 2283 k = ListGetInteger(Model % Bodies(Element % BodyId) % Values,'Material', & 2284 minv=1, maxv=Model % NumberOFMaterials) 2285 2286 CALL ListGetRealArray( Model % Materials(k) % Values, & 2287 'Enthalpy Heat Diffusivity', Hwrk, En, Edge % NodeIndexes ) 2288 2289 NodalConductivity( 1:En ) = Hwrk( 1,1,1:En ) 2290 2291! elementwise nodal solution: 2292! --------------------------- 2293 Enthalpy_h(1:Pn) = Quant( Perm(Element % NodeIndexes) ) 2294 2295! do the integration: 2296! ------------------- 2297 EdgeLength = 0.0d0 2298 ResidualNorm = 0.0d0 2299 2300 IntegStuff = GaussPoints( Edge ) 2301 2302 DO t=1,IntegStuff % n 2303 u = IntegStuff % u(t) 2304 v = IntegStuff % v(t) 2305 w = IntegStuff % w(t) 2306 2307 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 2308 EdgeBasis, dBasisdx ) 2309 2310 Normal = NormalVector( Edge, EdgeNodes, u, v, .TRUE. ) 2311 2312 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2313 s = IntegStuff % s(t) * detJ 2314 ELSE 2315 gx = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 2316 gy = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 2317 gz = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 2318 2319 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 2320 Symb, dSymb, gx, gy, gz ) 2321 2322 s = IntegStuff % s(t) * detJ * SqrtMetric 2323 END IF 2324 2325! 2326! Integration point in parent element local 2327! coordinates: 2328! ----------------------------------------- 2329 u = SUM( EdgeBasis(1:En) * x(1:En) ) 2330 v = SUM( EdgeBasis(1:En) * y(1:En) ) 2331 w = SUM( EdgeBasis(1:En) * z(1:En) ) 2332 2333 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 2334 Basis, dBasisdx ) 2335! 2336! Heat conductivity at the integration point: 2337! -------------------------------------------- 2338 Conductivity = SUM( NodalConductivity(1:En) * EdgeBasis(1:En) ) 2339! 2340! given flux at integration point: 2341! -------------------------------- 2342 Residual = -SUM( Flux(1:En) * EdgeBasis(1:En) ) 2343 2344! convective ...: 2345! ---------------- 2346 Residual = Residual + SUM(TransferCoeff(1:En) * EdgeBasis(1:En)) * & 2347 ( SUM( Enthalpy_h(1:Pn) * Basis(1:Pn) ) - & 2348 SUM( ExtEnthalpy_h(1:En) * EdgeBasis(1:En) ) ) 2349 2350! black body radiation...: 2351! ------------------------- 2352 Residual = Residual + & 2353 Emissivity * StefanBoltzmann * & 2354 ( SUM( Enthalpy_h(1:Pn) * Basis(1:Pn) ) ** 4 - & 2355 SUM( ExtEnthalpy_h(1:En) * EdgeBasis(1:En) ) ** 4 ) 2356 2357! flux given by the computed solution, and 2358! force norm for scaling the residual: 2359! ----------------------------------------- 2360 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2361 DO k=1,DIM 2362 Residual = Residual + Conductivity * & 2363 SUM( dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(k) 2364 2365 Gnorm = Gnorm + s * (Conductivity * & 2366 SUM(dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn)) * Normal(k))**2 2367 END DO 2368 ELSE 2369 DO k=1,DIM 2370 DO l=1,DIM 2371 Residual = Residual + Metric(k,l) * Conductivity * & 2372 SUM( dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(l) 2373 2374 Gnorm = Gnorm + s * (Metric(k,l) * Conductivity * & 2375 SUM(dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(l))**2 2376 END DO 2377 END DO 2378 END IF 2379 2380 EdgeLength = EdgeLength + s 2381 IF ( .NOT. Dirichlet ) THEN 2382 ResidualNorm = ResidualNorm + s * Residual ** 2 2383 END IF 2384 END DO 2385 EXIT 2386 END DO 2387 2388 IF ( CoordinateSystemDimension() == 3 ) THEN 2389 EdgeLength = SQRT(EdgeLength) 2390 END IF 2391 2392! Gnorm = EdgeLength * Gnorm 2393 Indicator = EdgeLength * ResidualNorm 2394 2395 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z) 2396 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 2397 2398 DEALLOCATE( Enthalpy_h, Basis, ExtEnthalpy_h, TransferCoeff, & 2399 x, y, z, EdgeBasis, dBasisdx, NodalConductivity, Flux, & 2400 NodalEmissivity ) 2401!------------------------------------------------------------------------------ 2402 END FUNCTION HeatBoundaryResidual 2403!------------------------------------------------------------------------------ 2404 2405 2406 2407!------------------------------------------------------------------------------ 2408 FUNCTION HeatEdgeResidual( Model, Edge, Mesh, Quant, Perm ) RESULT( Indicator ) 2409!------------------------------------------------------------------------------ 2410 USE DefUtils 2411 IMPLICIT NONE 2412 2413 TYPE(Model_t) :: Model 2414 INTEGER :: Perm(:) 2415 REAL(KIND=dp) :: Quant(:), Indicator(2) 2416 TYPE( Mesh_t ), POINTER :: Mesh 2417 TYPE( Element_t ), POINTER :: Edge 2418!------------------------------------------------------------------------------ 2419 2420 TYPE(Nodes_t) :: Nodes, EdgeNodes 2421 TYPE(Element_t), POINTER :: Element, Bndry 2422 2423 INTEGER :: i,j,k,l,n,t,DIM,En,Pn 2424 LOGICAL :: stat, Found 2425 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 2426 2427 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 2428 2429 REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:), x(:), y(:), z(:), & 2430 EdgeBasis(:), Basis(:), dBasisdx(:,:), Enthalpy_h(:) 2431 2432 REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength, Jump, Conductivity 2433 2434 REAL(KIND=dp) :: u, v, w, s, detJ 2435 2436 REAL(KIND=dp) :: Residual, ResidualNorm, Area 2437 2438 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 2439 2440 LOGICAL :: First = .TRUE. 2441 SAVE Hwrk, First 2442!------------------------------------------------------------------------------ 2443 2444 ! Initialize: 2445 ! ----------- 2446 IF ( First ) THEN 2447 First = .FALSE. 2448 NULLIFY( Hwrk ) 2449 END IF 2450 2451 SELECT CASE( CurrentCoordinateSystem() ) 2452 CASE( AxisSymmetric, CylindricSymmetric ) 2453 DIM = 3 2454 CASE DEFAULT 2455 DIM = CoordinateSystemDimension() 2456 END SELECT 2457 2458 Metric = 0.0d0 2459 DO i = 1,3 2460 Metric(i,i) = 1.0d0 2461 END DO 2462 2463 Grad = 0.0d0 2464! 2465! --------------------------------------------- 2466 2467 Element => Edge % BoundaryInfo % Left 2468 n = Element % TYPE % NumberOfNodes 2469 2470 Element => Edge % BoundaryInfo % Right 2471 n = MAX( n, Element % TYPE % NumberOfNodes ) 2472 2473 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 2474 2475 En = Edge % TYPE % NumberOfNodes 2476 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 2477 2478 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 2479 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 2480 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 2481 2482 ALLOCATE( NodalConductivity(En), EdgeBasis(En), Basis(n), & 2483 dBasisdx(n,3), x(En), y(En), z(En), Enthalpy_h(n) ) 2484 2485! Integrate square of jump over edge: 2486! ----------------------------------- 2487 ResidualNorm = 0.0d0 2488 EdgeLength = 0.0d0 2489 Indicator = 0.0d0 2490 2491 IntegStuff = GaussPoints( Edge ) 2492 2493 DO t=1,IntegStuff % n 2494 2495 u = IntegStuff % u(t) 2496 v = IntegStuff % v(t) 2497 w = IntegStuff % w(t) 2498 2499 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 2500 EdgeBasis, dBasisdx ) 2501 2502 Normal = NormalVector( Edge, EdgeNodes, u, v, .FALSE. ) 2503 2504 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2505 s = IntegStuff % s(t) * detJ 2506 ELSE 2507 u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 2508 v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 2509 w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 2510 2511 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 2512 Symb, dSymb, u, v, w ) 2513 s = IntegStuff % s(t) * detJ * SqrtMetric 2514 END IF 2515 2516 ! 2517 ! Compute flux over the edge as seen by elements 2518 ! on both sides of the edge: 2519 ! ---------------------------------------------- 2520 DO i = 1,2 2521 SELECT CASE(i) 2522 CASE(1) 2523 Element => Edge % BoundaryInfo % Left 2524 CASE(2) 2525 Element => Edge % BoundaryInfo % Right 2526 END SELECT 2527! 2528! Can this really happen (maybe it can...) ? 2529! ------------------------------------------- 2530 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) CYCLE 2531! 2532! Next, get the integration point in parent 2533! local coordinates: 2534! ----------------------------------------- 2535 Pn = Element % TYPE % NumberOfNodes 2536 2537 DO j = 1,En 2538 DO k = 1,Pn 2539 IF ( Edge % NodeIndexes(j) == Element % NodeIndexes(k) ) THEN 2540 x(j) = Element % TYPE % NodeU(k) 2541 y(j) = Element % TYPE % NodeV(k) 2542 z(j) = Element % TYPE % NodeW(k) 2543 EXIT 2544 END IF 2545 END DO 2546 END DO 2547 2548 u = SUM( EdgeBasis(1:En) * x(1:En) ) 2549 v = SUM( EdgeBasis(1:En) * y(1:En) ) 2550 w = SUM( EdgeBasis(1:En) * z(1:En) ) 2551! 2552! Get parent element basis & derivatives at the integration point: 2553! ----------------------------------------------------------------- 2554 Nodes % x(1:Pn) = Mesh % Nodes % x(Element % NodeIndexes) 2555 Nodes % y(1:Pn) = Mesh % Nodes % y(Element % NodeIndexes) 2556 Nodes % z(1:Pn) = Mesh % Nodes % z(Element % NodeIndexes) 2557 2558 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 2559 Basis, dBasisdx ) 2560! 2561! Material parameters: 2562! -------------------- 2563 k = ListGetInteger( Model % Bodies( & 2564 Element % BodyId) % Values, 'Material', & 2565 minv=1, maxv=Model % NumberOFMaterials ) 2566 2567 CALL ListGetRealArray( Model % Materials(k) % Values, & 2568 'Enthalpy Heat Diffusivity', Hwrk,En, Edge % NodeIndexes ) 2569 2570 NodalConductivity( 1:En ) = Hwrk( 1,1,1:En ) 2571 Conductivity = SUM( NodalConductivity(1:En) * EdgeBasis(1:En) ) 2572! 2573! Enthalpy_h at element nodal points: 2574! ------------------------------------ 2575 Enthalpy_h(1:Pn) = Quant( Perm(Element % NodeIndexes) ) 2576! 2577! Finally, the flux: 2578! ------------------ 2579 DO j=1,DIM 2580 Grad(j,i) = Conductivity * SUM( dBasisdx(1:Pn,j) * Enthalpy_h(1:Pn) ) 2581 END DO 2582 END DO 2583 2584! Compute squre of the flux jump: 2585! ------------------------------- 2586 EdgeLength = EdgeLength + s 2587 Jump = 0.0d0 2588 DO k=1,DIM 2589 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2590 Jump = Jump + (Grad(k,1) - Grad(k,2)) * Normal(k) 2591 ELSE 2592 DO l=1,DIM 2593 Jump = Jump + & 2594 Metric(k,l) * (Grad(k,1) - Grad(k,2)) * Normal(l) 2595 END DO 2596 END IF 2597 END DO 2598 ResidualNorm = ResidualNorm + s * Jump ** 2 2599 END DO 2600 2601 IF ( CoordinateSystemDimension() == 3 ) THEN 2602 EdgeLength = SQRT(EdgeLength) 2603 END IF 2604 Indicator = EdgeLength * ResidualNorm 2605 2606 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z, x, y, z) 2607 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 2608 DEALLOCATE( NodalConductivity, EdgeBasis, Basis, dBasisdx, Enthalpy_h) 2609 2610!------------------------------------------------------------------------------ 2611 END FUNCTION HeatEdgeResidual 2612!------------------------------------------------------------------------------ 2613 2614 2615!------------------------------------------------------------------------------ 2616 FUNCTION HeatInsideResidual( Model, Element, Mesh, & 2617 Quant, Perm, Fnorm ) RESULT( Indicator ) 2618!------------------------------------------------------------------------------ 2619 USE DefUtils 2620!------------------------------------------------------------------------------ 2621 IMPLICIT NONE 2622!------------------------------------------------------------------------------ 2623 TYPE(Model_t) :: Model 2624 INTEGER :: Perm(:) 2625 REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm 2626 TYPE( Mesh_t ), POINTER :: Mesh 2627 TYPE( Element_t ), POINTER :: Element 2628!------------------------------------------------------------------------------ 2629 2630 TYPE(Nodes_t) :: Nodes 2631 2632 INTEGER :: i,j,k,l,n,t,DIM 2633 2634 LOGICAL :: stat, Found, Compressible 2635 TYPE( Variable_t ), POINTER :: Var 2636 2637 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 2638 2639 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 2640 2641 REAL(KIND=dp), ALLOCATABLE :: NodalDensity(:) 2642 REAL(KIND=dp), ALLOCATABLE :: NodalCapacity(:) 2643 REAL(KIND=dp), ALLOCATABLE :: x(:), y(:), z(:) 2644 REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:) 2645 REAL(KIND=dp), ALLOCATABLE :: Velo(:,:), Pressure(:) 2646 REAL(KIND=dp), ALLOCATABLE :: NodalSource(:), Enthalpy_h(:), PrevTemp(:) 2647 REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:), ddBasisddx(:,:,:) 2648 2649 REAL(KIND=dp) :: u, v, w, s, detJ, Density, Capacity 2650 2651 REAL(KIND=dp) :: SpecificHeatRatio, ReferencePressure, dt 2652 REAL(KIND=dp) :: Source, Residual, ResidualNorm, Area, Conductivity 2653 2654 TYPE( ValueList_t ), POINTER :: Material 2655 2656 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 2657 2658 LOGICAL :: First = .TRUE. 2659 SAVE Hwrk, First 2660!------------------------------------------------------------------------------ 2661 2662! Initialize: 2663! ----------- 2664 Indicator = 0.0d0 2665 Fnorm = 0.0d0 2666! 2667! Check if this eq. computed in this element: 2668! ------------------------------------------- 2669 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 2670 2671 IF ( First ) THEN 2672 First = .FALSE. 2673 NULLIFY( Hwrk ) 2674 END IF 2675 2676 Metric = 0.0d0 2677 DO i=1,3 2678 Metric(i,i) = 1.0d0 2679 END DO 2680 2681 SELECT CASE( CurrentCoordinateSystem() ) 2682 CASE( AxisSymmetric, CylindricSymmetric ) 2683 DIM = 3 2684 CASE DEFAULT 2685 DIM = CoordinateSystemDimension() 2686 END SELECT 2687! 2688! Element nodal points: 2689! --------------------- 2690 n = Element % TYPE % NumberOfNodes 2691 2692 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 2693 2694 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 2695 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 2696 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 2697 2698 ALLOCATE( NodalDensity(n), NodalCapacity(n), NodalConductivity(n), & 2699 Velo(3,n), Pressure(n), NodalSource(n), Enthalpy_h(n), PrevTemp(n), & 2700 Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3) ) 2701! 2702! Elementwise nodal solution: 2703! --------------------------- 2704 Enthalpy_h(1:n) = Quant( Perm(Element % NodeIndexes) ) 2705! 2706! Check for time dep. 2707! ------------------- 2708 PrevTemp(1:n) = Enthalpy_h(1:n) 2709 dt = Model % Solver % dt 2710 IF ( ListGetString( Model % Simulation, 'Simulation Type') == 'transient' ) THEN 2711 Var => VariableGet( Model % Variables, 'Enthalpy_h', .TRUE. ) 2712 PrevTemp(1:n) = Var % PrevValues(Var % Perm(Element % NodeIndexes),1) 2713 END IF 2714! 2715! Material parameters: conductivity, heat capacity and density 2716! ------------------------------------------------------------- 2717 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', & 2718 minv=1, maxv=Model % NumberOfMaterials ) 2719 2720 Material => Model % Materials(k) % Values 2721 2722 CALL ListGetRealArray( Material, & 2723 'Enthalpy Heat Diffusivity', Hwrk,n, Element % NodeIndexes ) 2724 2725 NodalConductivity( 1:n ) = Hwrk( 1,1,1:n ) 2726 2727 NodalDensity(1:n) = ListGetReal( Material, & 2728 'Enthalpy Density', n, Element % NodeIndexes, Found ) 2729 2730 NodalCapacity(1:n) = 1.0 2731! 2732! Check for compressible flow equations: 2733! -------------------------------------- 2734 Compressible = .FALSE. 2735 2736 IF ( ListGetString( Material, 'Compressibility Model', Found ) == & 2737 'perfect gas equation 1' ) THEN 2738 2739 Compressible = .TRUE. 2740 2741 Pressure = 0.0d0 2742 Var => VariableGet( Mesh % Variables, 'Pressure', .TRUE. ) 2743 IF ( ASSOCIATED( Var ) ) THEN 2744 Pressure(1:n) = & 2745 Var % Values( Var % Perm(Element % NodeIndexes) ) 2746 END IF 2747 2748 ReferencePressure = ListGetConstReal( Material, & 2749 'Reference Pressure' ) 2750 2751 SpecificHeatRatio = ListGetConstReal( Material, & 2752 'Specific Heat Ratio' ) 2753 2754 NodalDensity(1:n) = (Pressure(1:n) + ReferencePressure) * SpecificHeatRatio / & 2755 ( (SpecificHeatRatio - 1) * NodalCapacity(1:n) * Enthalpy_h(1:n) ) 2756 END IF 2757! 2758! Get (possible) convection velocity at the nodes of the element: 2759! ---------------------------------------------------------------- 2760 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Equation', & 2761 minv=1, maxv=Model % NumberOFEquations ) 2762 2763 Velo = 0.0d0 2764 SELECT CASE( ListGetString( Model % Equations(k) % Values, & 2765 'Convection', Found ) ) 2766 2767 !----------------- 2768 CASE( 'constant' ) 2769 !----------------- 2770 2771 Velo(1,1:n) = ListGetReal( Material, & 2772 'Convection Velocity 1', n, Element % NodeIndexes, Found ) 2773 2774 Velo(2,1:n) = ListGetReal( Material, & 2775 'Convection Velocity 2', n, Element % NodeIndexes, Found ) 2776 2777 Velo(3,1:n) = ListGetReal( Material, & 2778 'Convection Velocity 3', n, Element % NodeIndexes, Found ) 2779 2780 !----------------- 2781 CASE( 'computed' ) 2782 !----------------- 2783 2784 Var => VariableGet( Mesh % Variables, 'Velocity 1', .TRUE. ) 2785 IF ( ASSOCIATED( Var ) ) THEN 2786 IF ( ALL( Var % Perm( Element % NodeIndexes ) > 0 ) ) THEN 2787 Velo(1,1:n) = Var % Values(Var % Perm(Element % NodeIndexes)) 2788 2789 Var => VariableGet( Mesh % Variables, 'Velocity 2', .TRUE. ) 2790 IF ( ASSOCIATED( Var ) ) & 2791 Velo(2,1:n) = Var % Values( & 2792 Var % Perm(Element % NodeIndexes ) ) 2793 2794 Var => VariableGet( Mesh % Variables, 'Velocity 3', .TRUE. ) 2795 IF ( ASSOCIATED( Var ) ) & 2796 Velo(3,1:n) = Var % Values( & 2797 Var % Perm( Element % NodeIndexes ) ) 2798 END IF 2799 END IF 2800 2801 END SELECT 2802 2803! 2804! Heat source: 2805! ------------ 2806! 2807 k = ListGetInteger( & 2808 Model % Bodies(Element % BodyId) % Values,'Body Force',Found, & 2809 1, Model % NumberOFBodyForces) 2810 2811 NodalSource = 0.0d0 2812 IF ( Found .AND. k > 0 ) THEN 2813 NodalSource(1:n) = ListGetReal( Model % BodyForces(k) % Values, & 2814 'Heat Source', n, Element % NodeIndexes, Found ) 2815 END IF 2816 2817! 2818! Integrate square of residual over element: 2819! ------------------------------------------ 2820 2821 ResidualNorm = 0.0d0 2822 Area = 0.0d0 2823 2824 IntegStuff = GaussPoints( Element ) 2825 2826 DO t=1,IntegStuff % n 2827 u = IntegStuff % u(t) 2828 v = IntegStuff % v(t) 2829 w = IntegStuff % w(t) 2830 2831 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 2832 Basis, dBasisdx, ddBasisddx, .TRUE., .FALSE. ) 2833 2834 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2835 s = IntegStuff % s(t) * detJ 2836 ELSE 2837 u = SUM( Basis(1:n) * Nodes % x(1:n) ) 2838 v = SUM( Basis(1:n) * Nodes % y(1:n) ) 2839 w = SUM( Basis(1:n) * Nodes % z(1:n) ) 2840 2841 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 2842 Symb, dSymb, u, v, w ) 2843 s = IntegStuff % s(t) * detJ * SqrtMetric 2844 END IF 2845 2846 Capacity = SUM( NodalCapacity(1:n) * Basis(1:n) ) 2847 Density = SUM( NodalDensity(1:n) * Basis(1:n) ) 2848 Conductivity = SUM( NodalConductivity(1:n) * Basis(1:n) ) 2849! 2850! Residual of the convection-diffusion (heat) equation: 2851! R = \rho * c_p * (@T/@t + u.grad(T)) - & 2852! div(C grad(T)) + p div(u) - h, 2853! --------------------------------------------------- 2854! 2855! or more generally: 2856! 2857! R = \rho * c_p * (@T/@t + u^j T_{,j}) - & 2858! g^{jk} (C T_{,j}}_{,k} + p div(u) - h 2859! --------------------------------------------------- 2860! 2861 Residual = -Density * SUM( NodalSource(1:n) * Basis(1:n) ) 2862 2863 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2864 DO j=1,DIM 2865! 2866! - grad(C).grad(T): 2867! -------------------- 2868! 2869 Residual = Residual - & 2870 SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) ) * & 2871 SUM( NodalConductivity(1:n) * dBasisdx(1:n,j) ) 2872 2873! 2874! - C div(grad(T)): 2875! ------------------- 2876! 2877 Residual = Residual - Conductivity * & 2878 SUM( Enthalpy_h(1:n) * ddBasisddx(1:n,j,j) ) 2879 END DO 2880 ELSE 2881 DO j=1,DIM 2882 DO k=1,DIM 2883! 2884! - g^{jk} C_{,k}T_{j}: 2885! --------------------- 2886! 2887 Residual = Residual - Metric(j,k) * & 2888 SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) ) * & 2889 SUM( NodalConductivity(1:n) * dBasisdx(1:n,k) ) 2890 2891! 2892! - g^{jk} C T_{,jk}: 2893! ------------------- 2894! 2895 Residual = Residual - Metric(j,k) * Conductivity * & 2896 SUM( Enthalpy_h(1:n) * ddBasisddx(1:n,j,k) ) 2897! 2898! + g^{jk} C {_jk^l} T_{,l}: 2899! --------------------------- 2900 DO l=1,DIM 2901 Residual = Residual + Metric(j,k) * Conductivity * & 2902 Symb(j,k,l) * SUM( Enthalpy_h(1:n) * dBasisdx(1:n,l) ) 2903 END DO 2904 END DO 2905 END DO 2906 END IF 2907 2908! + \rho * c_p * (@T/@t + u.grad(T)): 2909! ----------------------------------- 2910 Residual = Residual + Density * Capacity * & 2911 SUM((Enthalpy_h(1:n)-PrevTemp(1:n))*Basis(1:n)) / dt 2912 2913 DO j=1,DIM 2914 Residual = Residual + & 2915 Density * Capacity * SUM( Velo(j,1:n) * Basis(1:n) ) * & 2916 SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) ) 2917 END DO 2918 2919 2920 IF ( Compressible ) THEN 2921! 2922! + p div(u) or p u^j_{,j}: 2923! ------------------------- 2924! 2925 DO j=1,DIM 2926 Residual = Residual + & 2927 SUM( Pressure(1:n) * Basis(1:n) ) * & 2928 SUM( Velo(j,1:n) * dBasisdx(1:n,j) ) 2929 2930 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2931 DO k=1,DIM 2932 Residual = Residual + & 2933 SUM( Pressure(1:n) * Basis(1:n) ) * & 2934 Symb(j,k,j) * SUM( Velo(k,1:n) * Basis(1:n) ) 2935 END DO 2936 END IF 2937 END DO 2938 END IF 2939 2940! 2941! Compute also force norm for scaling the residual: 2942! ------------------------------------------------- 2943 DO i=1,DIM 2944 Fnorm = Fnorm + s * ( Density * & 2945 SUM( NodalSource(1:n) * Basis(1:n) ) ) ** 2 2946 END DO 2947 2948 Area = Area + s 2949 ResidualNorm = ResidualNorm + s * Residual ** 2 2950 END DO 2951 2952! Fnorm = Element % hk**2 * Fnorm 2953 Indicator = Element % hK**2 * ResidualNorm 2954 2955 DEALLOCATE( NodalDensity, NodalCapacity, NodalConductivity, & 2956 Velo, Pressure, NodalSource, Enthalpy_h, PrevTemp, Basis, & 2957 dBasisdx, ddBasisddx ) 2958!------------------------------------------------------------------------------ 2959 END FUNCTION HeatInsideResidual 2960!------------------------------------------------------------------------------ 2961 2962