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