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! * Authors: Juha Ruokolainen 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 08 Jun 1997 34! * 35! ****************************************************************************/ 36 37!> \ingroup Solvers 38!> \{ 39 40!------------------------------------------------------------------------------ 41!> Initialization of the main solver: AdvectionDiffusionSolver 42!------------------------------------------------------------------------------ 43 SUBROUTINE FlowSolver_init( Model,Solver,Timestep,TransientSimulation ) 44!------------------------------------------------------------------------------ 45 USE DefUtils 46 IMPLICIT NONE 47!------------------------------------------------------------------------------ 48 TYPE(Solver_t) :: Solver !< Linear & nonlinear equation solver options 49 TYPE(Model_t), TARGET :: Model !< All model information (mesh, materials, BCs, etc...) 50 REAL(KIND=dp) :: Timestep !< Timestep size for time dependent simulations 51 LOGICAL :: TransientSimulation !< Steady state or transient simulation 52!------------------------------------------------------------------------------ 53! Local variables 54!------------------------------------------------------------------------------ 55 TYPE(ValueList_t), POINTER :: Params 56 57 Params => GetSolverParams() 58 CALL ListAddInteger( Params,'Time Derivative Order',1 ) 59 60 END SUBROUTINE FlowSolver_init 61 62 63!------------------------------------------------------------------------------ 64!> Solver for the Navier-Stokes equation in various different coordinate systems. 65!------------------------------------------------------------------------------ 66 SUBROUTINE FlowSolver( Model,Solver,dt,TransientSimulation) 67!------------------------------------------------------------------------------ 68 USE NavierStokes 69 USE NavierStokesGeneral 70 USE NavierStokesCylindrical 71 USE Adaptive 72 USE DefUtils 73 USE FreeSurface 74 USE ElementDescription, ONLY: GetEdgeMap 75!------------------------------------------------------------------------------ 76 IMPLICIT NONE 77 78 TYPE(Model_t) :: Model 79 TYPE(Solver_t), TARGET :: Solver 80 81 REAL(KIND=dp) :: dt 82 LOGICAL :: TransientSimulation 83!------------------------------------------------------------------------------ 84! Local variables 85!------------------------------------------------------------------------------ 86 TYPE(Matrix_t),POINTER :: StiffMatrix 87 88 INTEGER :: i,j,k,n,nb,nd,t,iter,LocalNodes,istat,q,m 89 90 TYPE(ValueList_t),POINTER :: Material, BC, BodyForce, Equation 91 TYPE(Nodes_t) :: ElementNodes 92 TYPE(Element_t),POINTER :: Element 93 94 REAL(KIND=dp) :: RelativeChange,UNorm,Gravity(3),AngularVelocity(3), & 95 Tdiff,s,Relaxation,NewtonTol,NewtonUBound,NonlinearTol, & 96 ReferencePressure=0.0, SpecificHeatRatio, & 97 PseudoCompressibilityScale=1.0, FreeSTol, res 98 99 INTEGER :: NSDOFs,NewtonIter,NewtonMaxIter,NonlinearIter,FreeSIter 100 101 TYPE(Variable_t), POINTER :: DensitySol, TimeVar 102 TYPE(Variable_t), POINTER :: FlowSol, TempSol, MeshSol 103 104 INTEGER, POINTER :: FlowPerm(:),TempPerm(:), MeshPerm(:) 105 REAL(KIND=dp), POINTER :: FlowSolution(:), Temperature(:), & 106 gWork(:,:), ForceVector(:), LayerThickness(:), & 107 SurfaceRoughness(:),MeshVelocity(:) 108 109 REAL(KIND=dp), POINTER :: TempPrev(:) 110 REAL(KIND=DP), POINTER :: Pwrk(:,:,:) 111 112 LOGICAL :: Stabilize,NewtonLinearization = .FALSE., GotForceBC, GotIt, & 113 MBFlag, Convect = .TRUE., NormalTangential, RelaxBefore, & 114 divDiscretization, GradPDiscretization, ComputeFree=.FALSE., & 115 Transient, Rotating, AnyRotating, OutOfPlaneFlow=.FALSE.,& 116 RecheckNewton=.FALSE. 117 118! Which compressibility model is used 119 CHARACTER(LEN=MAX_NAME_LEN) :: CompressibilityFlag, StabilizeFlag, VarName 120 CHARACTER(LEN=MAX_NAME_LEN) :: LocalCoords, FlowModel 121 INTEGER :: CompressibilityModel, ModelCoords, ModelDim, NoActive 122 INTEGER :: body_id,bf_id,eq_id,DIM 123 INTEGER :: MidEdgeNodes(12), BrickFaceMap(6,4) 124 INTEGER, POINTER :: NodeIndexes(:), Indexes(:) 125 INTEGER, POINTER :: EdgeMap(:,:) 126 127 128 INTEGER, SAVE :: Timestep, SaveTimestep=-1 129 REAL(KIND=dp), ALLOCATABLE, SAVE :: pDensity0(:), pDensity1(:) 130! 131 LOGICAL :: AllocationsDone = .FALSE., FreeSurfaceFlag, & 132 PseudoPressureExists, PseudoCompressible, Bubbles, P2P1, & 133 Porous =.FALSE., PotentialForce=.FALSE., Hydrostatic=.FALSE., & 134 MagneticForce =.FALSE., UseLocalCoords, PseudoPressureUpdate, & 135 AllIncompressible 136 137 138 REAL(KIND=dp),ALLOCATABLE :: MASS(:,:),STIFF(:,:), LoadVector(:,:), & 139 Viscosity(:),FORCE(:), TimeForce(:), PrevDensity(:),Density(:), & 140 U(:),V(:),W(:),MU(:),MV(:),MW(:), Pressure(:),Alpha(:),Beta(:), & 141 ExtPressure(:),PrevPressure(:), HeatExpansionCoeff(:), & 142 ReferenceTemperature(:), Permeability(:),Mx(:),My(:),Mz(:), & 143 LocalTemperature(:), GasConstant(:), HeatCapacity(:), & 144 LocalTempPrev(:),SlipCoeff(:,:), PseudoCompressibility(:), & 145 PseudoPressure(:), Drag(:,:), PotentialField(:), & 146 PotentialCoefficient(:) 147 148 SAVE U,V,W,MASS,STIFF,LoadVector,Viscosity, TimeForce,FORCE,ElementNodes, & 149 Alpha,Beta,ExtPressure,Pressure,PrevPressure, PrevDensity,Density, & 150 AllocationsDone,LocalNodes, HeatExpansionCoeff,ReferenceTemperature, & 151 Permeability,Mx,My,Mz,LayerThickness, SlipCoeff, SurfaceRoughness, & 152 LocalTemperature, GasConstant, HeatCapacity, LocalTempPrev,MU,MV,MW, & 153 PseudoCompressibilityScale, PseudoCompressibility, PseudoPressure, & 154 PseudoPressureExists, Drag, PotentialField, PotentialCoefficient, & 155 ComputeFree, Indexes 156 157 REAL(KIND=dp) :: at,at0,at1,totat,st,totst 158!------------------------------------------------------------------------------ 159 160 INTERFACE 161 FUNCTION FlowBoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm ) RESULT(Indicator) 162 USE Types 163 TYPE(Element_t), POINTER :: Edge 164 TYPE(Model_t) :: Model 165 TYPE(Mesh_t), POINTER :: Mesh 166 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 167 INTEGER :: Perm(:) 168 END FUNCTION FlowBoundaryResidual 169 170 FUNCTION FlowEdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT(Indicator) 171 USE Types 172 TYPE(Element_t), POINTER :: Edge 173 TYPE(Model_t) :: Model 174 TYPE(Mesh_t), POINTER :: Mesh 175 REAL(KIND=dp) :: Quant(:), Indicator(2) 176 INTEGER :: Perm(:) 177 END FUNCTION FlowEdgeResidual 178 179 FUNCTION FlowInsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm ) RESULT(Indicator) 180 USE Types 181 TYPE(Element_t), POINTER :: Element 182 TYPE(Model_t) :: Model 183 TYPE(Mesh_t), POINTER :: Mesh 184 REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm 185 INTEGER :: Perm(:) 186 END FUNCTION FlowInsideResidual 187 END INTERFACE 188!------------------------------------------------------------------------------ 189 190!------------------------------------------------------------------------------ 191! Get variables needed for solving the system 192!------------------------------------------------------------------------------ 193 CALL Info('FlowSolver','Solving the Navier-Stokes equations',Level=6) 194 195 196 IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN 197 198 CALL DefaultStart() 199 200! Check for local coordinate system 201 202 LocalCoords = GetString( Solver % Values, 'Solver Coordinate System', & 203 UseLocalCoords ) 204 205 IF ( UseLocalCoords ) THEN 206 ModelCoords = Coordinates 207 ModelDim = Model % DIMENSION 208 SELECT CASE ( LocalCoords ) 209 CASE( 'cartesian 2d' ) 210 Coordinates = 1 211 Model % DIMENSION = 2 212 CALL Info( 'FlowSolve', 'Solver Coordinate System is cartesian 2d', LEVEL=31 ) 213 CASE( 'cartesian 3d' ) 214 Coordinates = 1 215 Model % DIMENSION = 3 216 CALL Info( 'FlowSolve', 'Solver Coordinate System is cartesian 3d', LEVEL=31 ) 217 CASE( 'axi symmetric' ) 218 Coordinates = 4 219 Model % DIMENSION = 2 220 CALL Info( 'FlowSolve', 'Solver Coordinate System is axi symmetric', LEVEL=31 ) 221 CASE( 'cylindric symmetric' ) 222 Coordinates = 3 223 Model % DIMENSION = 3 224 CALL Info( 'FlowSolve', 'Solver Coordinate System is cylindric symmetric', LEVEL=31 ) 225 CASE DEFAULT 226 CALL Warn( 'FlowSolve', 'Solver coordinate system not recognised, using original' ) 227 END SELECT 228 END IF 229 230 ! check for Flow model, one of 'full', 'no convection', 'stokes': 231 ! --------------------------------------------------------------- 232 Transient = TransientSimulation 233 Convect = .TRUE. 234 FlowModel = GetString( GetSolverParams(), 'Flow Model', Gotit ) 235 236 SELECT CASE(FlowModel) 237 CASE('no convection') 238 Convect = .FALSE. 239 CASE('stokes') 240 Convect = .FALSE. 241 Transient = .FALSE. 242 CASE DEFAULT 243 FlowModel = 'full' 244 END SELECT 245 246 DIM = CoordinateSystemDimension() 247 248 FlowSol => Solver % Variable 249 NSDOFs = FlowSol % DOFs 250 FlowPerm => FlowSol % Perm 251 FlowSolution => FlowSol % Values 252 253 VarName = GetVarName(FlowSol) 254 255 LocalNodes = COUNT( FlowPerm > 0 ) 256 IF ( LocalNodes <= 0 ) RETURN 257 258 TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' ) 259 IF ( ASSOCIATED( TempSol ) ) THEN 260 TempPerm => TempSol % Perm 261 Temperature => TempSol % Values 262 IF( Transient ) THEN 263 IF ( ASSOCIATED(TempSol % PrevValues) ) TempPrev => TempSol % PrevValues(:,1) 264 END IF 265 END IF 266 267 MeshSol => VariableGet( Solver % Mesh % Variables, 'Mesh Velocity') 268 NULLIFY( MeshVelocity ) 269 IF ( ASSOCIATED( MeshSol ) ) THEN 270 MeshPerm => MeshSol % Perm 271 MeshVelocity => MeshSol % Values 272 END IF 273 274 DensitySol => VariableGet( Solver % Mesh % Variables, 'Density' ) 275 276!------------------------------------------------------------------------------ 277 278 StiffMatrix => Solver % Matrix 279 ForceVector => StiffMatrix % RHS 280 UNorm = Solver % Variable % Norm 281 282!------------------------------------------------------------------------------ 283! Allocate some permanent storage, this is done first time only 284!------------------------------------------------------------------------------ 285 286 IF ( .NOT.AllocationsDone .OR. Solver % MeshChanged ) THEN 287 288 N = Solver % Mesh % MaxElementDOFs 289 290 IF( AllocationsDone ) THEN 291 DEALLOCATE( & 292 U, V, W, & 293 MU, MV, MW, & 294 Indexes, & 295 Pressure, & 296 PrevPressure, & 297 PseudoCompressibility, & 298 PrevDensity,Density, & 299 LayerThickness, & 300 SurfaceRoughness, & 301 Permeability, & 302 Mx,My,Mz, & 303 SlipCoeff, Drag, & 304 TimeForce,FORCE, Viscosity, & 305 MASS, STIFF, & 306 HeatExpansionCoeff, & 307 GasConstant, HeatCapacity, & 308 ReferenceTemperature, & 309 LocalTempPrev, LocalTemperature, & 310 PotentialField, PotentialCoefficient, & 311 LoadVector, Alpha, Beta, & 312 ExtPressure, STAT=istat ) 313 END IF 314 315 ALLOCATE( U(N), V(N), W(N), & 316 MU(N), MV(N), MW(N), & 317 Indexes( N ), & 318 Pressure( N ), & 319 PrevPressure( N ), & 320 PseudoCompressibility( N ), & 321 PrevDensity(N),Density( N ), & 322 LayerThickness(N), & 323 SurfaceRoughness(N), & 324 Permeability(N), & 325 Mx(N),My(N),Mz(N), & 326 SlipCoeff(3,N), Drag(3,N), & 327 TimeForce( 2*NSDOFs*N ), & 328 FORCE( 2*NSDOFs*N ), Viscosity( N ), & 329 MASS( 2*NSDOFs*N,2*NSDOFs*N ),& 330 STIFF( 2*NSDOFs*N,2*NSDOFs*N ),& 331 HeatExpansionCoeff(N), & 332 GasConstant( N ), HeatCapacity( N ), & 333 ReferenceTemperature(N), & 334 LocalTempPrev(N), LocalTemperature(N), & 335 PotentialField( N ), PotentialCoefficient( N ), & 336 LoadVector( 4,N ), Alpha( N ), Beta( N ), & 337 ExtPressure( N ), STAT=istat ) 338 339 Drag = 0.0d0 340 NULLIFY(Pwrk) 341 342 PseudoPressureExists = .FALSE. 343 AllIncompressible = .TRUE. 344 345 DO k=1,Model % NumberOfMaterials 346 Material => Model % Materials(k) % Values 347 CompressibilityFlag = ListGetString( Material, & 348 'Compressibility Model', GotIt) 349 IF (.NOT. gotIt ) CYCLE 350 IF( CompressibilityFlag == 'artificial compressible') THEN 351 PseudoPressureExists = .TRUE. 352 AllIncompressible = .FALSE. 353 ELSE IF ( CompressibilityFlag == 'incompressible' ) THEN 354 CONTINUE 355 ELSE 356 AllIncompressible = .FALSE. 357 END IF 358 END DO 359 360 IF ( PseudoPressureExists ) THEN 361 IF ( AllocationsDone ) THEN 362 DEALLOCATE( PseudoPressure ) 363 END IF 364 n = SIZE( FlowSolution ) / NSDOFs 365 ALLOCATE( PseudoPressure(n),STAT=istat ) 366 END IF 367 368!------------------------------------------------------------------------------ 369! This hack is needed cause of the fluctuating pressure levels 370!------------------------------------------------------------------------------ 371 IF( AllIncompressible ) THEN 372 CALL Info('FlowSolve','Enforcing relative pressure relaxation',Level=8) 373 CALL ListAddNewLogical( Solver % Values,'Relative Pressure Relaxation',.TRUE.) 374 END IF 375 376 IF ( istat /= 0 ) THEN 377 CALL Fatal( 'FlowSolve','Memory allocation error, Aborting.' ) 378 END IF 379 380!------------------------------------------------------------------------------ 381 382 AllocationsDone = .TRUE. 383 END IF 384!------------------------------------------------------------------------------ 385 386 387 TimeVar => VariableGet( Solver % Mesh % Variables, 'Timestep') 388 Timestep = NINT(Timevar % Values(1)) 389 IF ( SaveTimestep /= Timestep ) THEN 390 IF ( ALLOCATED(pDensity0) ) pDensity0 = pDensity1 391 SaveTimestep=Timestep 392 END IF 393 394!------------------------------------------------------------------------------ 395! Do some additional initialization, and go for it 396!------------------------------------------------------------------------------ 397 398 gWork => ListGetConstRealArray( Model % Constants,'Gravity',GotIt) 399 IF ( GotIt ) THEN 400 Gravity = gWork(1:3,1)*gWork(4,1) 401 ELSE 402 Gravity = 0.00D0 403 Gravity(2) = -9.81D0 404 END IF 405 406 AnyRotating = ListCheckPresentAnyBodyForce(Model,'Angular Velocity') .OR. & 407 ListCheckPresentAnyBodyForce(Model,'Angular Velocity 1') .OR. & 408 ListCheckPresentAnyBodyForce(Model,'Angular Velocity 2') .OR. & 409 ListCheckPresentAnyBodyForce(Model,'Angular Velocity 3') 410 411!------------------------------------------------------------------------------ 412 P2P1 = .FALSE. 413 Bubbles = ListGetLogical( Solver % Values,'Bubbles',GotIt ) 414 Stabilize = ListGetLogical( Solver % Values,'Stabilize',GotIt ) 415 416 StabilizeFlag = ListGetString( Solver % Values, & 417 'Stabilization Method', GotIt ) 418 IF ( .NOT. GotIt ) THEN 419 IF ( Stabilize ) THEN 420 StabilizeFlag = 'stabilized' 421 ELSE IF ( Bubbles ) THEN 422 StabilizeFlag = 'bubbles' 423 ELSE 424 StabilizeFlag = 'stabilized' 425 END IF 426 ELSE 427 IF (StabilizeFlag == 'p2/p1' .OR. StabilizeFlag == 'p2p1') THEN 428 P2P1 = .TRUE. 429 Bubbles = .FALSE. 430 Stabilize = .FALSE. 431 END IF 432 END IF 433 434 IF ( StabilizeFlag == 'bubbles' ) Bubbles = .TRUE. 435 436 DivDiscretization = ListGetLogical( Solver % Values, & 437 'Div Discretization', GotIt ) 438 439 GradPDiscretization = ListGetLogical( Solver % Values, & 440 'Gradp Discretization', GotIt ) 441 442 NonlinearTol = ListGetConstReal( Solver % Values, & 443 'Nonlinear System Convergence Tolerance',minv=0.0d0 ) 444 445 NewtonTol = ListGetConstReal( Solver % Values, & 446 'Nonlinear System Newton After Tolerance', minv=0.0d0 ) 447 448 !Option to switch back to picard if convergence exceeds certain tolerance 449 NewtonUBound = ListGetConstReal( Solver % Values, & 450 'Nonlinear System Newton Max Tolerance', GotIt ) 451 IF(GotIt) RecheckNewton = .TRUE. 452 453 NewtonIter = ListGetInteger( Solver % Values, & 454 'Nonlinear System Newton After Iterations', minv=0 ) 455 IF ( NewtonIter == 0 ) NewtonLinearization = .TRUE. 456 457 !Option to switch back to picard after NewtonMaxIter iterations 458 NewtonMaxIter = ListGetInteger( Solver % Values, & 459 'Nonlinear System Newton Max Iterations', GotIt ) 460 RecheckNewton = RecheckNewton .OR. GotIt 461 462 IF (GetLogical( GetSolverParams(), & 463 'Nonlinear System Reset Newton', GotIt)) NewtonLinearization=.FALSE. 464 465 NonlinearIter = ListGetInteger( Solver % Values, & 466 'Nonlinear System Max Iterations', minv=0 ) 467 468 IF ( .NOT. ListCheckPresent( Solver % Values, & 469 'Nonlinear System Norm Dofs' ) ) THEN 470 CALL ListAddInteger( Solver % Values, 'Nonlinear System Norm DOFs', NSDOFs-1 ) 471 END IF 472 473 FreeSTol = ListGetConstReal( Solver % Values, & 474 'Free Surface After Tolerance', GotIt, minv=0.0d0 ) 475 IF ( .NOT. GotIt ) FreeSTol = HUGE(1.0d0) 476 477 FreeSIter = ListGetInteger( Solver % Values, & 478 'Free Surface After Iterations', GotIt, minv=0 ) 479 IF ( .NOT. GotIt ) FreeSIter = 0 480 481!------------------------------------------------------------------------------ 482! Check if free surfaces present 483!------------------------------------------------------------------------------ 484 FreeSurfaceFlag = .FALSE. 485 DO i=1,Model % NumberOfBCs 486 FreeSurfaceFlag = FreeSurfaceFlag .OR. ListGetLogical( & 487 Model % BCs(i) % Values,'Free Surface', GotIt ) 488 IF ( FreeSurfaceFlag ) EXIT 489 END DO 490 491 CALL CheckCircleBoundary() 492!------------------------------------------------------------------------------ 493 494 totat = 0.0d0 495 totst = 0.0d0 496 497 498 ! Initialize the pressure to be used in artificial compressibility 499 IF(PseudoPressureExists) THEN 500 PseudoPressure = FlowSolution(NSDOFs:SIZE(FlowSolution):NSDOFs) 501 502 WRITE(Message,'(A,T25,E15.4)') 'PseudoPressure mean: ',& 503 SUM(PseudoPressure)/SIZE(PseudoPressure) 504 CALL Info('FlowSolve',Message,Level=5) 505 506 PseudoCompressibilityScale = ListGetConstReal( Model % Simulation, & 507 'Artificial Compressibility Scaling',GotIt) 508 509 IF(.NOT.GotIt) PseudoCompressibilityScale = 1.0 510 IF(Transient) THEN 511 PseudoCompressibilityScale = PseudoCompressibilityScale / dt 512 END IF 513 PseudoPressureUpdate = ListGetLogical( Model % Simulation, & 514 'Pseudo Pressure Update',GotIt) 515 IF (.NOT.GotIt) PseudoPressureUpdate = .FALSE. 516 END IF 517 518 DO iter=1,NonlinearIter 519 520 IF (PseudoPressureExists .AND. PseudoPressureUpdate) & 521 PseudoPressure = FlowSolution(NSDOFs:SIZE(FlowSolution):NSDOFs) 522 523 at = CPUTime() 524 at0 = RealTime() 525 at1 = RealTime() 526 527 CALL Info( 'FlowSolve', ' ', Level=4 ) 528 CALL Info( 'FlowSolve', ' ', Level=4 ) 529 CALL Info( 'FlowSolve', '-------------------------------------', Level=4 ) 530 WRITE( Message, * ) 'NAVIER-STOKES ITERATION', iter 531 CALL Info( 'FlowSolve',Message, Level=4 ) 532 CALL Info( 'FlowSolve','-------------------------------------', Level=4 ) 533 CALL Info( 'FlowSolve', ' ', Level=4 ) 534 CALL Info( 'FlowSolve','Starting Assembly...', Level=4 ) 535 536!------------------------------------------------------------------------------ 537 !CALL InitializeToZero( StiffMatrix, ForceVector ) 538 CALL DefaultInitialize() 539!------------------------------------------------------------------------------ 540 541 bf_id = -1 542 body_id = -1 543 544 CALL StartAdvanceOutput( 'FlowSolve', 'Assembly: ' ) 545 NoActive = GetNOFActive() 546 547 DO t = 1,NoActive 548 549 CALL AdvanceOutput( t, NoActive ) 550! 551 Element => GetActiveElement(t) 552 NodeIndexes => Element % NodeIndexes 553 554!------------------------------------------------------------------------------ 555 556 IF ( Element % BodyId /= body_id ) THEN 557 body_id = Element % BodyId 558 559 eq_id = ListGetInteger( Model % Bodies(body_id) % Values,'Equation', & 560 minv=1, maxv=Model % NumberOfEquations ) 561 Equation => Model % Equations(eq_id) % Values 562 563 bf_id = ListGetInteger( Model % Bodies(body_id) % Values, & 564 'Body Force', gotIt, 1, Model % NumberOfBodyForces ) 565 IF( bf_id > 0 ) THEN 566 BodyForce => Model % BodyForces(bf_id) % Values 567 END IF 568 569 IF ( FlowModel == 'full' ) THEN 570 Convect = ListGetLogical( Equation,'NS Convect', GotIt ) 571 IF ( .NOT. GotIt ) Convect = .TRUE. 572 ENDIF 573 574 k = ListGetInteger( Model % Bodies(body_id) % Values, 'Material', & 575 minv=1, maxv=Model % NumberOfMaterials ) 576 Material => Model % Materials(k) % Values 577 578!------------------------------------------------------------------------------ 579 CompressibilityFlag = ListGetString( Material, & 580 'Compressibility Model', GotIt) 581 IF ( .NOT.GotIt ) CompressibilityModel = Incompressible 582 PseudoCompressible = .FALSE. 583!------------------------------------------------------------------------------ 584 SELECT CASE( CompressibilityFlag ) 585!------------------------------------------------------------------------------ 586 CASE( 'incompressible' ) 587 CompressibilityModel = Incompressible 588 589 CASE( 'perfect gas', 'perfect gas equation 1' ) 590 CompressibilityModel = PerfectGas1 591 592 CASE( 'thermal' ) 593 CompressibilityModel = Thermal 594 595 CASE( 'user defined' ) 596 CompressibilityModel = UserDefined1 597 598 CASE( 'pressure dependent' ) 599 CompressibilityModel = UserDefined2 600 601 CASE( 'artificial compressible' ) 602 CompressibilityModel = Incompressible 603 PseudoCompressible = .TRUE. 604 605 CASE DEFAULT 606 CompressibilityModel = Incompressible 607!------------------------------------------------------------------------------ 608 END SELECT 609!------------------------------------------------------------------------------ 610 611 Gotit = .FALSE. 612 IF ( bf_id > 0 ) THEN 613 MagneticForce = ListGetLogical( BodyForce,'Lorentz Force', gotIt ) 614 Hydrostatic = ListGetLogical( BodyForce,'Hydrostatic Pressure',gotIt ) 615 END IF 616 IF ( .NOT. GotIt ) THEN 617 Hydrostatic = ListGetLogical( Equation,'Hydrostatic Pressure',gotIt ) 618 END IF 619!------------------------------------------------------------------------------ 620 621 Rotating = .FALSE. 622 IF( bf_id > 0 ) THEN 623 gWork => ListGetConstRealArray( BodyForce,'Angular Velocity',GotIt) 624 IF ( GotIt ) THEN 625 IF( Coordinates == Cartesian ) THEN 626 AngularVelocity = gWork(1:3,1) 627 Rotating = .TRUE. 628 ELSE 629 CALL Fatal('FlowSolve','Rotating coordinate implemented only for cartesian coodinates') 630 END IF 631 ELSE 632 AngularVelocity = 0.0_dp 633 END IF 634 END IF 635 END IF 636!------------------------------------------------------------------------------ 637 638 n = GetElementNOFNodes() 639 nb = GetElementNOFBDOFs() 640 nd = GetElementDOFs( Indexes ) 641 642 CALL GetElementNodes( ElementNodes ) 643 644 SELECT CASE( NSDOFs ) 645 CASE(3) 646 U(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-2) 647 V(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-1) 648 W(1:nd) = 0.0_dp 649 IF (bf_id > 0 ) THEN 650 W(1:n) = ListGetReal(BodyForce,'Out Of Plane Velocity',& 651 n, NodeIndexes(1:n),OutOfPlaneFlow) 652 IF (.NOT.OutOfPlaneFlow) & 653 W(1:n) = 0.0_dp 654 END IF 655 CASE(4) 656 U(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-3) 657 V(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-2) 658 W(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-1) 659 END SELECT 660 661 MU(1:nd) = 0.0d0 662 MV(1:nd) = 0.0d0 663 MW(1:nd) = 0.0d0 664 IF ( ASSOCIATED( MeshVelocity ) ) THEN 665 SELECT CASE( MeshSol % DOFs ) 666 CASE(2) 667 IF ( ALL( MeshPerm( Indexes(1:nd) ) > 0 ) ) THEN 668 MU(1:nd) = MeshVelocity(2*MeshPerm(Indexes(1:nd))-1) 669 MV(1:nd) = MeshVelocity(2*MeshPerm(Indexes(1:nd))-0) 670 END IF 671 672 CASE(3) 673 IF ( ALL( MeshPerm( NodeIndexes ) > 0 ) ) THEN 674 MU(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-2) 675 MV(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-1) 676 MW(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-0) 677 END IF 678 END SELECT 679 END IF 680 681 LocalTemperature = 0.0d0 682 LocalTempPrev = 0.0d0 683 IF ( ASSOCIATED( TempSol ) ) THEN 684 IF ( ALL( TempPerm(NodeIndexes) > 0 ) ) THEN 685 LocalTemperature(1:nd) = Temperature( TempPerm(Indexes(1:nd)) ) 686 IF ( Transient .AND. CompressibilityModel /= Incompressible) THEN 687 LocalTempPrev(1:nd) = TempPrev( TempPerm(Indexes(1:nd)) ) 688 END IF 689 END IF 690 END IF 691 ReferencePressure = 0.0d0 692 693 PrevDensity = 0.0d0 694 Density = 0.0d0 695!------------------------------------------------------------------------------ 696 SELECT CASE( CompressibilityModel ) 697!------------------------------------------------------------------------------ 698 CASE( Incompressible ) 699!------------------------------------------------------------------------------ 700 Pressure(1:nd) = FlowSolution( NSDOFs*FlowPerm(Indexes(1:nd)) ) 701 Density(1:n) = GetReal( Material, 'Density' ) 702 703 IF(PseudoCompressible) THEN 704 Pressure(1:n) = GetReal( Material,'Artificial Pressure', GotIt ) 705 IF(.NOT. GotIt) THEN 706 Pressure(1:nd) = PseudoPressure(FlowPerm(Indexes(1:nd))) 707 ELSE 708 Pressure(n+1:nd) = 0.0d0 709 END IF 710 PseudoCompressibility(1:n) = PseudoCompressibilityScale * & 711 GetReal(Material,'Artificial Compressibility', GotIt ) 712 IF(.NOT. gotIt) PseudoCompressibility(1:n) = 0.0d0 713 END IF 714 715!------------------------------------------------------------------------------ 716 CASE( PerfectGas1 ) 717 718 ! Use ReferenceTemperature in .sif file for fixed temperature 719 ! field. At the moment can not have both fixed T ideal gas and 720 ! Boussinesq force: 721 !------------------------------------------------------------- 722 IF ( .NOT. ASSOCIATED( TempSol ) ) THEN 723 LocalTemperature(1:n) = GetReal( Material, & 724 'Reference Temperature' ) 725 LocalTempPrev = LocalTemperature 726 END IF 727 728 HeatCapacity(1:n) = GetReal( Material, 'Heat Capacity', GotIt ) 729 730 731 ! Read Specific Heat Ratio: 732 !-------------------------- 733 SpecificHeatRatio = ListGetConstReal( Material, & 734 'Specific Heat Ratio', GotIt ) 735 IF ( .NOT.GotIt ) SpecificHeatRatio = 5.d0/3.d0 736 737 738 ! For an ideal gas, \gamma, c_p and R are really a constant 739 ! GasConstant is an array only since HeatCapacity formally is 740 !------------------------------------------------------------ 741 GasConstant(1:n) = ( SpecificHeatRatio - 1.d0 ) * & 742 HeatCapacity(1:n) / SpecificHeatRatio 743 744 745 ! For ideal gases take pressure deviation p_d as the 746 ! dependent variable: p = p_0 + p_d 747 ! Read p_0 748 !--------------------------------------------------- 749 ReferencePressure = ListGetConstReal( Material, & 750 'Reference Pressure', GotIt ) 751 IF ( .NOT.GotIt ) ReferencePressure = 0.0d0 752 753 Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))) 754 IF ( Transient ) THEN 755 PrevPressure(1:nd) = Solver % Variable % PrevValues( & 756 NSDOFs*FlowPerm(Indexes(1:nd)),1 ) 757 END IF 758 Density(1:n) = ( Pressure(1:n) + ReferencePressure ) / & 759 ( GasConstant(1:n) * LocalTemperature(1:n) ) 760 761 CASE( UserDefined1 ) 762 Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)) ) 763 IF ( ASSOCIATED( DensitySol ) ) THEN 764 Density(1:nd) = DensitySol % Values( DensitySol % Perm(Indexes(1:nd)) ) 765 IF ( Transient ) THEN 766 PrevDensity(1:nd) = DensitySol % PrevValues( & 767 DensitySol % Perm(Indexes(1:nd)),1) 768 END IF 769 ELSE 770 Density(1:n) = GetReal( Material,'Density' ) 771 IF ( Transient ) THEN 772 IF (.NOT.ALLOCATED(pDensity0)) THEN 773 ALLOCATE(pDensity0(LocalNodes), & 774 pDensity1(LocalNodes)) 775 END IF 776 777 IF ( Timestep==1 ) & 778 pDensity0(Indexes(1:n)) = Density(1:n) 779 pDensity1(Indexes(1:n)) = Density(1:n) 780 PrevDensity(1:n) = pDensity0(Indexes(1:n)) 781 END IF 782 END IF 783 784 CASE( UserDefined2 ) 785 Density(1:n) = GetReal( Material,'Density' ) 786 Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))) 787 788 CASE( Thermal ) 789 Pressure(1:n) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))) 790 791 HeatExpansionCoeff(1:n) = GetReal( Material, & 792 'Heat Expansion Coefficient' ) 793 794 ReferenceTemperature(1:n) = GetReal( Material, & 795 'Reference Temperature' ) 796 797 Density(1:n) = GetReal( Material,'Density' ) 798 799 IF( Transient ) THEN 800 PrevDensity(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n) * & 801 ( LocalTempPrev(1:n) - ReferenceTemperature(1:n) ) ) 802 END IF 803 Density(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n) * & 804 ( LocalTemperature(1:n) - ReferenceTemperature(1:n) ) ) 805!------------------------------------------------------------------------------ 806 END SELECT 807!------------------------------------------------------------------------------ 808 809!------------------------------------------------------------------------------ 810! Read in porous media defs 811!------------------------------------------------------------------------------ 812 Porous = ListGetLogical( Material,'Porous Media', GotIt) 813 IF(Porous) THEN 814 CALL GetRealArray( Material, Pwrk,'Porous Resistivity',GotIt) 815 816 IF( .NOT. GotIt ) THEN 817 Drag( 1,1:n) = GetReal( Material,'Porous Resistivity 1',GotIt ) 818 Drag( 2,1:n) = GetReal( Material,'Porous Resistivity 2',GotIt ) 819 IF( NSDOFs -1 > 2 ) THEN 820 Drag( 3,1:n) = GetReal( Material,'Porous Resistivity 3',GotIt ) 821 END IF 822 ELSE IF ( SIZE(Pwrk,1) == 1 ) THEN 823 DO i=1,NSDOFs-1 824 Drag( i,1:n ) = Pwrk( 1,1,1:n ) 825 END DO 826 ELSE 827 DO i=1,MIN(NSDOFs,SIZE(Pwrk,1)) 828 Drag(i,1:n) = Pwrk(i,1,1:n) 829 END DO 830 END IF 831 END IF 832 833!------------------------------------------------------------------------------ 834! Viscosity = Laminar viscosity 835!------------------------------------------------------------------------------ 836 Viscosity(1:n) = GetReal( Material,'Viscosity' ) 837 838!------------------------------------------------------------------------------ 839! Set body forces, if any 840!------------------------------------------------------------------------------ 841 LoadVector = 0.0D0 842 843 IF ( bf_id > 0 ) THEN 844 HeatExpansionCoeff = 0.0D0 845 ReferenceTemperature = 0.0D0 846 847!------------------------------------------------------------------------------ 848! Boussinesq body force & gravity 849!------------------------------------------------------------------------------ 850 IF ( ListGetLogical( BodyForce,'Boussinesq',gotIt) ) THEN 851 852 HeatExpansionCoeff(1:n) = GetReal( Material, & 853 'Heat Expansion Coefficient' ) 854 855 ReferenceTemperature(1:n) = GetReal( Material, & 856 'Reference Temperature' ) 857 858 DO i=1,n 859 k = TempPerm(NodeIndexes(i)) 860 IF ( k > 0 ) THEN 861 IF ( Hydrostatic ) THEN 862 Tdiff = 1 - HeatExpansionCoeff(i) * & 863 (Temperature(k) - ReferenceTemperature(i)) 864 865 IF ( Tdiff <= 0.0D0 ) THEN 866 CALL Warn( 'FlowSolve','Zero or negative density.' ) 867 END IF 868 ELSE 869 Tdiff = -HeatExpansionCoeff(i) * & 870 (Temperature(k) - ReferenceTemperature(i)) 871 END IF 872 873 LoadVector(1,i) = Gravity(1) * Tdiff 874 LoadVector(2,i) = Gravity(2) * Tdiff 875 IF ( NSDOFs > 3 ) THEN 876 LoadVector(3,i) = Gravity(3) * Tdiff 877 END IF 878 END IF 879 END DO 880 ELSE IF ( Hydrostatic ) THEN 881 LoadVector(1,1:n) = Gravity(1) 882 LoadVector(2,1:n) = Gravity(2) 883 IF ( NSDOFs > 3 ) LoadVector(3,1:n) = Gravity(3) 884 END IF 885!------------------------------------------------------------------------------ 886 LoadVector(1,1:n) = LoadVector(1,1:n) + ListGetReal( BodyForce, & 887 'Flow Bodyforce 1',n,NodeIndexes,gotIt ) 888 889 LoadVector(2,1:n) = LoadVector(2,1:n) + ListGetReal( BodyForce, & 890 'Flow Bodyforce 2',n,NodeIndexes,gotIt ) 891 892 IF ( NSDOFs > 3 ) THEN 893 LoadVector(3,1:n) = LoadVector(3,1:n) + ListGetReal( BodyForce, & 894 'Flow Bodyforce 3',n,NodeIndexes,gotIt ) 895 END IF 896 897!------------------------------------------------------------------------------ 898 899 PotentialForce = ListGetLogical( BodyForce,'Potential Force',gotIt) 900 IF(PotentialForce) THEN 901 PotentialField(1:n) = ListGetReal( BodyForce, & 902 'Potential Field',n,NodeIndexes) 903 PotentialCoefficient(1:n) = ListGetReal( BodyForce, & 904 'Potential Coefficient',n,NodeIndexes) 905 END IF 906 907 908 909!------------------------------------------------------------------------------ 910 END IF ! of body forces 911 912!------------------------------------------------------------------------------ 913! 914! NOTE: LoadVector is multiplied by density inside *Navier* routines 915! 916 IF ( Transient ) THEN 917 SELECT CASE( CompressibilityModel ) 918 CASE( PerfectGas1 ) 919 IF ( ASSOCIATED( TempSol ) ) THEN 920 DO i=1,n 921 k = TempPerm(NodeIndexes(i)) 922 IF ( k > 0 ) THEN 923 LoadVector(NSDOFs,i) = LoadVector(NSDOFs,i) + & 924 ( Temperature(k) - TempPrev(k) ) / dt 925 END IF 926 END DO 927 END IF 928 CASE( UserDefined1, Thermal ) 929 DO i=1,n 930 LoadVector(NSDOFs,i) = LoadVector(NSDOFs,i) - & 931 ( Density(i) - PrevDensity(i) ) / (Density(i)*dt) 932 END DO 933 END SELECT 934 END IF 935 936!------------------------------------------------------------------------------ 937! Get element local stiffness & mass matrices 938!------------------------------------------------------------------------------ 939 SELECT CASE(Coordinates) 940 CASE( Cartesian ) 941!------------------------------------------------------------------------------ 942 SELECT CASE( CompressibilityModel ) 943!------------------------------------------------------------------------------ 944 CASE( Incompressible,PerfectGas1,UserDefined1,UserDefined2,Thermal) 945!------------------------------------------------------------------------------ 946! Density needed for steady-state, also pressure for transient 947!------------------------------------------------------------------------------ 948 CALL NavierStokesCompose( MASS,STIFF,FORCE, LoadVector, & 949 Viscosity,Density,U,V,W,MU,MV,MW,ReferencePressure+Pressure(1:n), & 950 LocalTemperature, Convect, StabilizeFlag, CompressibilityModel, & 951 PseudoCompressible, PseudoCompressibility, GasConstant, Porous, & 952 Drag, PotentialForce, PotentialField, PotentialCoefficient, & 953 MagneticForce, Rotating, AngularVelocity, DivDiscretization, & 954 GradPDiscretization, NewtonLinearization, Element,n,ElementNodes) 955!------------------------------------------------------------------------------ 956 END SELECT 957!------------------------------------------------------------------------------ 958 959 CASE( Cylindric,CylindricSymmetric,AxisSymmetric ) 960! Same comments as Cartesian 961!------------------------------------------------------------------------------ 962 SELECT CASE( CompressibilityModel ) 963!------------------------------------------------------------------------------ 964 CASE( Incompressible,PerfectGas1) 965!------------------------------------------------------------------------------ 966 CALL NavierStokesCylindricalCompose( & 967 MASS,STIFF,FORCE, & 968 LoadVector, Viscosity,Density,U,V,W,MU,MV,MW, & 969 ReferencePressure+Pressure(1:n),LocalTemperature,& 970 Convect, StabilizeFlag, CompressibilityModel /= Incompressible, & 971 PseudoCompressible, PseudoCompressibility, GasConstant, Porous, Drag, & 972 PotentialForce, PotentialField, PotentialCoefficient, & 973 MagneticForce, divDiscretization, gradpDiscretization, NewtonLinearization, & 974 Element,n,ElementNodes ) 975!------------------------------------------------------------------------------ 976 CASE DEFAULT 977 CALL Fatal('FlowSolver','Missing compressibility model in cylindrical coordinates') 978 979 980 END SELECT 981!------------------------------------------------------------------------------ 982 983!------------------------------------------------------------------------------ 984 CASE DEFAULT 985!------------------------------------------------------------------------------ 986 987 SELECT CASE( CompressibilityModel ) 988!------------------------------------------------------------------------------ 989 CASE( Incompressible,PerfectGas1) 990 991 CALL NavierStokesGeneralCompose( & 992 MASS,STIFF,FORCE, & 993 LoadVector, Viscosity,Density,U,V,W,MU,MV,MW,Stabilize, & 994 NewtonLinearization,Element,n,ElementNodes ) 995 996 CASE DEFAULT 997 CALL Fatal('FlowSolver','Missing compressibility model in general coordinates') 998 999 END SELECT 1000 1001!------------------------------------------------------------------------------ 1002 END SELECT 1003!------------------------------------------------------------------------------ 1004! If time dependent simulation, add mass matrix to global 1005! matrix and global RHS vector 1006!------------------------------------------------------------------------------ 1007 IF ( CompressibilityModel /= Incompressible .AND. & 1008 StabilizeFlag == 'stabilized' ) THEN 1009 Bubbles = .TRUE. 1010 StabilizeFlag = 'bubbles' 1011 END IF 1012 IF ( Element % TYPE % BasisFunctionDegree <= 1 .AND. P2P1 ) THEN 1013 Bubbles = .TRUE. 1014 StabilizeFlag = 'bubbles' 1015 END IF 1016 1017 IF ( nb==0 .AND. Bubbles ) nb = n 1018 1019 TimeForce = 0.0_dp 1020 IF ( Transient ) THEN 1021!------------------------------------------------------------------------------ 1022! NOTE: the following will replace STIFF and FORCE 1023! with the combined information 1024!------------------------------------------------------------------------------ 1025 CALL Default1stOrderTime( MASS, STIFF, FORCE ) 1026 END IF 1027 1028 IF ( nb > 0 ) THEN 1029 CALL NSCondensate( nd, nb, NSDOFs-1, STIFF, FORCE, TimeForce ) 1030 END IF 1031 1032!------------------------------------------------------------------------------ 1033! Add local stiffness matrix and force vector to global matrix & vector 1034!------------------------------------------------------------------------------ 1035 CALL DefaultUpdateEquations( STIFF, FORCE ) 1036!------------------------------------------------------------------------------ 1037 END DO 1038!------------------------------------------------------------------------------ 1039 1040 CALL DefaultFinishBulkAssembly() 1041 1042 CALL Info( 'FlowSolve', 'Assembly done', Level=4 ) 1043 1044!------------------------------------------------------------------------------ 1045! Neumann & Newton boundary conditions 1046!------------------------------------------------------------------------------ 1047 NoActive = GetNOFBoundaryElements() 1048 1049 DO t = 1,NoActive 1050 1051 Element => GetBoundaryElement(t) 1052 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 1053 1054 n = GetElementNOFNodes() 1055 1056 CALL GetElementNodes( ElementNodes ) 1057 NodeIndexes => Element % NodeIndexes 1058 1059 BC => GetBC() 1060 IF ( .NOT. ASSOCIATED(BC) ) CYCLE 1061 1062!------------------------------------------------------------------------------ 1063 GotForceBC = GetLogical( BC, 'Flow Force BC',gotIt ) 1064 IF ( .NOT. gotIt ) GotForceBC = .TRUE. 1065 1066 IF ( GotForceBC ) THEN 1067 LoadVector = 0.0d0 1068 Alpha = 0.0d0 1069 ExtPressure = 0.0d0 1070 Beta = 0.0d0 1071 SlipCoeff = 0.0d0 1072 STIFF = 0.0d0 1073 FORCE = 0.0d0 1074 1075!------------------------------------------------------------------------------ 1076! (at the moment the following is done...) 1077! BC: \tau \cdot n = \alpha n + @\beta/@t + R_k u_k + F 1078!------------------------------------------------------------------------------ 1079!------------------------------------------------------------------------------ 1080! normal force BC: \tau\cdot n = \alpha n 1081!------------------------------------------------------------------------------ 1082 IF ( GetLogical( BC, 'Free Surface',gotIt) ) THEN 1083 Alpha(1:n) = GetReal( BC,'Surface Tension Coefficient', gotIt ) 1084 END IF 1085 1086 ExtPressure(1:n) = GetReal( BC, 'External Pressure', GotForceBC ) 1087 IF(.NOT. GotForceBC) ExtPressure(1:n) = GetReal( BC, 'Normal Pressure', GotForceBC ) 1088!------------------------------------------------------------------------------ 1089! tangential force BC: 1090! \tau\cdot n = @\beta/@t (tangential derivative of something) 1091!------------------------------------------------------------------------------ 1092 1093 IF ( ASSOCIATED( TempSol ) ) THEN 1094 Beta(1:n) = GetReal( BC, & 1095 'Surface Tension Expansion Coefficient',gotIt ) 1096 1097 IF ( gotIt ) THEN 1098 DO j=1,n 1099 k = TempPerm( NodeIndexes(j) ) 1100 IF ( k>0 ) Beta(j) = 1.0_dp - Beta(j) * Temperature(k) 1101 END DO 1102 Beta(1:n) = Beta(1:n) * GetReal(BC, 'Surface Tension Coefficient' ) 1103 ELSE 1104 Beta(1:n) = GetReal( BC,'Surface Tension Coefficient', gotIt ) 1105 END IF 1106 END IF 1107 1108!------------------------------------------------------------------------------ 1109! force in given direction BC: \tau\cdot n = F 1110!------------------------------------------------------------------------------ 1111 1112 LoadVector(1,1:n) = GetReal( BC, 'Pressure 1', GotIt ) 1113 LoadVector(2,1:n) = GetReal( BC, 'Pressure 2', GotIt ) 1114 LoadVector(3,1:n) = GetReal( BC, 'Pressure 3', GotIt ) 1115 LoadVector(4,1:n) = GetReal( BC, 'Mass Flux', GotIt ) 1116 1117!------------------------------------------------------------------------------ 1118! slip boundary condition BC: \tau\cdot n = R_k u_k 1119!------------------------------------------------------------------------------ 1120 1121 SlipCoeff = 0.0d0 1122 SlipCoeff(1,1:n) = GetReal( BC, 'Slip Coefficient 1',GotIt ) 1123 SlipCoeff(2,1:n) = GetReal( BC, 'Slip Coefficient 2',GotIt ) 1124 SlipCoeff(3,1:n) = GetReal( BC, 'Slip Coefficient 3',GotIt ) 1125 1126 NormalTangential = GetLogical( BC, & 1127 'Normal-Tangential Velocity', GotIt ) 1128 1129 IF (.NOT.GotIt) THEN 1130 NormalTangential = GetLogical( BC, & 1131 'Normal-Tangential '//GetVarName(Solver % Variable), GotIt ) 1132 END IF 1133!------------------------------------------------------------------------------ 1134 SELECT CASE( CurrentCoordinateSystem() ) 1135 CASE( Cartesian ) 1136 1137 CALL NavierStokesBoundary( STIFF, FORCE, & 1138 LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, NormalTangential, & 1139 Element, n, ElementNodes ) 1140 1141 CASE( Cylindric, CylindricSymmetric, AxisSymmetric ) 1142 1143 CALL NavierStokesCylindricalBoundary( STIFF, & 1144 FORCE, LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, & 1145 NormalTangential, Element, n, ElementNodes) 1146 1147 CASE DEFAULT 1148 1149 CALL NavierStokesGeneralBoundary( STIFF, & 1150 FORCE, LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, & 1151 Element, n, ElementNodes) 1152 1153 END SELECT 1154 1155!------------------------------------------------------------------------------ 1156 1157 IF ( GetLogical( BC, 'Wall Law',GotIt ) ) THEN 1158 !/* 1159 ! * TODO: note that the following is not really valid, the 1160 ! * pointer to the Material structure is from the remains 1161 ! * of the last of the bulk elements. 1162 ! */ 1163 Density(1:n) = GetParentMatProp( 'Density' ) 1164 Viscosity(1:n) = GetParentMatProp( 'Viscosity' ) 1165 1166 LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness' ) 1167 SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt ) 1168 1169 SELECT CASE( NSDOFs ) 1170 CASE(3) 1171 U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 ) 1172 V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 ) 1173 W(1:n) = 0.0d0 1174 1175 CASE(4) 1176 U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-3 ) 1177 V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 ) 1178 W(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 ) 1179 END SELECT 1180 1181 CALL NavierStokesWallLaw( STIFF,FORCE, & 1182 LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, & 1183 Element,n, ElementNodes ) 1184 ELSE IF ( GetLogical( BC, 'VMS Wall', GotIt ) ) THEN 1185 Density(1:n) = GetParentMatProp( 'Density' ) 1186 Viscosity(1:n) = GetParentMatProp( 'Viscosity' ) 1187 1188 LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness', GotIt ) 1189 SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt ) 1190 1191 SELECT CASE( NSDOFs ) 1192 CASE(3) 1193 U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 ) 1194 V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 ) 1195 W(1:n) = 0.0d0 1196 1197 CASE(4) 1198 U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-3 ) 1199 V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 ) 1200 W(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 ) 1201 END SELECT 1202 CALL VMSWalls( STIFF,FORCE, & 1203 LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, & 1204 Element,n, ElementNodes ) 1205 1206 END IF 1207!------------------------------------------------------------------------------ 1208 1209 IF ( Transient ) THEN 1210 MASS = 0.0d0 1211 CALL Default1stOrderTime( MASS, STIFF, FORCE ) 1212 END IF 1213 1214!------------------------------------------------------------------------------ 1215! Add local stiffness matrix and force vector to 1216! global matrix & vector 1217!-------------------------------------------------------------------------- 1218 CALL DefaultUpdateEquations( STIFF, FORCE ) 1219!------------------------------------------------------------------------------ 1220 END IF 1221 END DO 1222!------------------------------------------------------------------------------ 1223 ! 1224 ! IMPLEMENT NOSLIP WALL BC CODE: 1225 ! ------------------------------ 1226 DO i=1,Model % NumberOFBCs 1227 BC => Model % BCs(i) % Values 1228 IF ( GetLogical( BC, 'Noslip wall BC', gotit ) ) THEN 1229 IF ( VarName == 'flow solution' ) THEN 1230 CALL ListAddConstReal( BC, 'Velocity 1', 0.0_dp ) 1231 CALL ListAddConstReal( BC, 'Velocity 2', 0.0_dp ) 1232 IF ( NSDOFs>3 ) CALL ListAddConstReal( BC, 'Velocity 3', 0.0_dp ) 1233 ELSE 1234 DO j=1,NSDOFs-1 1235 CALL ListAddConstReal( BC, ComponentName( & 1236 Solver % Variable % Name, j), 0.0_dp ) 1237 END DO 1238 END IF 1239 END IF 1240 END DO 1241 1242 CALL DefaultFinishBoundaryAssembly() 1243 CALL DefaultFinishAssembly() 1244 1245!------------------------------------------------------------------------------ 1246! Dirichlet boundary conditions 1247!------------------------------------------------------------------------------ 1248 CALL DefaultDirichletBCs() 1249 CALL Info( 'FlowSolve', 'Dirichlet conditions done', Level=4 ) 1250 1251!------------------------------------------------------------------------------ 1252! Solve the system and check for convergence 1253!------------------------------------------------------------------------------ 1254 at = CPUTime() - at 1255 st = CPUTime() 1256 1257 Unorm = DefaultSolve() 1258 1259 st = CPUTIme()-st 1260 totat = totat + at 1261 totst = totst + st 1262 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat 1263 CALL Info( 'FlowSolve', Message, Level=4 ) 1264 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve: (s)', st, totst 1265 CALL Info( 'FlowSolve', Message, Level=4 ) 1266 1267 n = NSDOFs * LocalNodes 1268 1269 1270!------------------------------------------------------------------------------ 1271 1272 WRITE( Message, * ) 'Result Norm : ',Solver % Variable % Norm 1273 CALL Info( 'FlowSolve', Message, Level=4 ) 1274 WRITE( Message, * ) 'Relative Change : ',Solver % Variable % NonlinChange 1275 CALL Info( 'FlowSolve', Message, Level=4 ) 1276 1277 RelativeChange = Solver % Variable % NonlinChange 1278 IF ( RelativeChange < NewtonTol .OR. & 1279 iter > NewtonIter ) NewtonLinearization = .TRUE. 1280 1281 IF ( RecheckNewton .AND. NewtonLinearization .AND. (RelativeChange > NewtonUBound)) THEN 1282 NewtonLinearization = .FALSE. 1283 CALL Info('FlowSolve', 'Newton tolerance exceeded, switching back to picard', Level=6) 1284 END IF 1285 1286 IF ( RecheckNewton .AND. NewtonLinearization .AND. (iter >= NewtonMaxIter)) THEN 1287 NewtonLinearization = .FALSE. 1288 CALL Info('FlowSolve', 'Newton iteration limit exceeded, switching back to picard', Level=6) 1289 END IF 1290 1291 IF ( RecheckNewton .AND. (RelativeChange > NewtonUBound) .AND. NewtonLinearization ) THEN 1292 NewtonLinearization = .FALSE. 1293 CALL Info('FlowSolve', 'Newton tolerance exceeded, switching back to picard', Level=6) 1294 END IF 1295 1296 IF ( RelativeChange < NonLinearTol .AND. Iter<NonlinearIter ) EXIT 1297 1298!------------------------------------------------------------------------------ 1299! If free surfaces in model, this will move the nodal points 1300!------------------------------------------------------------------------------ 1301 IF ( FreeSurfaceFlag ) THEN 1302 1303 IF ( RelativeChange < FreeSTol .OR. & 1304 iter > FreeSIter ) ComputeFree = .TRUE. 1305 1306 IF ( ComputeFree ) THEN 1307 Relaxation = GetCReal( Solver % Values, & 1308 'Free Surface Relaxation Factor', GotIt ) 1309 1310 IF ( .NOT.GotIt ) Relaxation = 1.0d0 1311 1312 MBFlag = GetLogical(GetSolverParams(), 'Internal Move Boundary', GotIt) 1313 IF ( MBFlag .OR. .NOT. GotIt ) CALL MoveBoundary( Model, Relaxation ) 1314 END IF 1315 END IF 1316!------------------------------------------------------------------------------ 1317 END DO ! of nonlinear iteration 1318 1319 IF ( P2P1 ) THEN 1320 !---------------------------------------------------------------------------------------- 1321 ! Replace the zero pressure solution at the nodes which are not needed in the linear 1322 ! pressure approximation by the interpolated values for right visualization: 1323 !---------------------------------------------------------------------------------------- 1324 NoActive = GetNOFActive() 1325 1326 DO t=1,NoActive 1327 ! First the midedge nodes: 1328 Element => GetActiveElement(t) 1329 IF ( Element % TYPE % BasisFunctionDegree <= 1 ) CYCLE 1330 1331 nd = GetElementDOFs( Indexes ) 1332 k = GetElementFamily() 1333 EdgeMap => GetEdgeMap(k) 1334 SELECT CASE(k) 1335 CASE (3) 1336 MidEdgeNodes(1:3) = (/ 4, 5, 6 /) 1337 CASE (4) 1338 MidEdgeNodes(1:4) = (/ 5, 6, 7, 8 /) 1339 CASE (5) 1340 MidEdgeNodes(1:6) = (/ 5, 6, 7, 8, 9, 10 /) 1341 CASE (6) 1342 MidEdgeNodes(1:8) = (/ 6, 7, 8, 9, 10, 11, 12, 13 /) 1343 CASE (7) 1344 MidEdgeNodes(1:9) = (/ 7, 8, 9, 10, 11, 12, 13, 14, 15 /) 1345 CASE (8) 1346 MidEdgeNodes(1:12) = (/ 9, 10, 11, 12, 17, 18, 19, 20, 13, 14, 15, 16 /) 1347 END SELECT 1348 1349 DO q=1,SIZE(EdgeMap,1) 1350 m = (dim+1) * Solver % Variable % Perm(Indexes(MidEdgeNodes(q))) 1351 i = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,1))) 1352 j = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,2))) 1353 Solver % Variable % Values(m) = 0.5d0 * ( Solver % Variable % Values(i) + & 1354 Solver % Variable % Values(j) ) 1355 END DO 1356 1357 ! The pressure at the midface nodes for 409 elements: 1358 IF (k==4 .AND. nd==9) THEN 1359 res = 0.0d0 1360 DO q=1,4 1361 i = (dim+1) * Solver % Variable % Perm(Indexes(q)) 1362 res = res + Solver % Variable % Values(i) 1363 END DO 1364 m = (dim+1) * Solver % Variable % Perm(Indexes(9)) 1365 Solver % Variable % Values(m) = 0.25d0 * res 1366 END IF 1367 1368 ! The pressure at the midpoint and at the midface nodes for 827 elements: 1369 IF (k==8 .AND. nd==27) THEN 1370 BrickFaceMap(1,:) = (/ 1,2,6,5 /) 1371 BrickFaceMap(2,:) = (/ 2,3,7,6 /) 1372 BrickFaceMap(3,:) = (/ 4,3,7,8 /) 1373 BrickFaceMap(4,:) = (/ 1,4,8,5 /) 1374 BrickFaceMap(5,:) = (/ 1,2,3,4 /) 1375 BrickFaceMap(6,:) = (/ 5,6,7,8 /) 1376 DO j=1,6 1377 res = 0.0d0 1378 DO q=1,4 1379 i = (dim+1) * Solver % Variable % Perm(Indexes(BrickFaceMap(j,q))) 1380 res = res + Solver % Variable % Values(i) 1381 END DO 1382 m = (dim+1) * Solver % Variable % Perm(Indexes(20+j)) 1383 Solver % Variable % Values(m) = 0.25d0 * res 1384 END DO 1385 1386 res = 0.0d0 1387 DO q=1,8 1388 i = (dim+1) * Solver % Variable % Perm(Indexes(q)) 1389 res = res + Solver % Variable % Values(i) 1390 END DO 1391 m = (dim+1) * Solver % Variable % Perm(Indexes(27)) 1392 Solver % Variable % Values(m) = 0.125d0 * res 1393 END IF 1394 1395 END DO 1396 END IF 1397 1398 IF (ListGetLogical(Solver % Values,'Adaptive Mesh Refinement',GotIt)) & 1399 CALL RefineMesh( Model,Solver,FlowSolution,FlowPerm, & 1400 FlowInsideResidual, FlowEdgeResidual, FlowBoundaryResidual ) 1401 1402!------------------------------------------------------------------------------ 1403 CALL CheckCircleBoundary() 1404!------------------------------------------------------------------------------ 1405 1406 IF ( UseLocalCoords ) THEN 1407 Coordinates = ModelCoords 1408 Model % DIMENSION = ModelDim 1409 END IF 1410 1411CONTAINS 1412!------------------------------------------------------------------------------ 1413 1414 1415!------------------------------------------------------------------------------ 1416 SUBROUTINE CheckCircleBoundary() 1417!------------------------------------------------------------------------------ 1418 REAL(KIND=dp) :: x,y,phi,x0,y0,r 1419 LOGICAL :: GotIt 1420 INTEGER :: i,j,k,l 1421!------------------------------------------------------------------------------ 1422 1423 l = 0 1424 DO i=1,Model % NumberOfBCs 1425 IF ( .NOT.ListgetLogical( Model % BCs(i) % Values, & 1426 'Circle Boundary', GotIt ) ) CYCLE 1427 1428 x0 = ListGetConstReal( Model % BCs(i) % Values, 'Circle X', GotIt ) 1429 IF ( .NOT. GotIt ) x0 = 0.0d0 1430 1431 y0 = ListGetConstReal( Model % BCs(i) % Values, 'Circle Y', GotIt ) 1432 IF ( .NOT. GotIt ) y0 = 0.0d0 1433 1434 R = ListGetConstReal( Model % BCs(i) % Values, 'Circle R', GotIt ) 1435 IF ( .NOT. GotIt ) R = 1.0d0 1436 1437 DO j=Solver % Mesh % NumberOfBulkElements+1, & 1438 Solver % Mesh % NumberOfBulkElements+ & 1439 Solver % Mesh % NumberOfBoundaryElements 1440 Element => Solver % Mesh % Elements(j) 1441 IF ( Element % BoundaryInfo % Constraint & 1442 /= Model % BCs(i) % Tag ) CYCLE 1443 1444 n = Element % TYPE % NumberOfNodes 1445 NodeIndexes => Element % NodeIndexes 1446 DO k=1,n 1447 x = Solver % Mesh % Nodes % x(NodeIndexes(k)) - x0 1448 y = Solver % Mesh % Nodes % y(NodeIndexes(k)) - y0 1449 1450 phi = ATAN2( y,x ) 1451 x = R * COS( phi ) 1452 y = R * SIN( phi ) 1453 1454 Solver % Mesh % Nodes % x(NodeIndexes(k)) = x + x0 1455 Solver % Mesh % Nodes % y(NodeIndexes(k)) = y + y0 1456 END DO 1457 l = l + 1 1458 END DO 1459 END DO 1460 1461 IF ( l > 0 ) THEN 1462 WRITE( Message, * ) 'Elements on Circle', l 1463 CALL Info( 'FlowSolve', Message, Level=6 ) 1464 END IF 1465!------------------------------------------------------------------------------ 1466 END SUBROUTINE CheckCircleBoundary 1467!------------------------------------------------------------------------------ 1468 1469!------------------------------------------------------------------------------ 1470 END SUBROUTINE FlowSolver 1471!------------------------------------------------------------------------------ 1472 1473!------------------------------------------------------------------------------ 1474!> Compute the residual of the Navier-Stokes equation for the boundary elements. 1475!------------------------------------------------------------------------------ 1476 FUNCTION FlowBoundaryResidual( Model, Edge, Mesh, & 1477 Quant, Perm, Gnorm ) RESULT( Indicator ) 1478!------------------------------------------------------------------------------ 1479 USE DefUtils 1480 IMPLICIT NONE 1481!------------------------------------------------------------------------------ 1482 TYPE(Model_t) :: Model 1483 INTEGER :: Perm(:) 1484 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 1485 TYPE( Mesh_t ), POINTER :: Mesh 1486 TYPE( Element_t ), POINTER :: Edge 1487!------------------------------------------------------------------------------ 1488 1489 TYPE(Nodes_t) :: Nodes, EdgeNodes 1490 TYPE(Element_t), POINTER :: Element 1491 1492 INTEGER :: i,k,n,l,t,bc,DIM,DOFs,Pn,En 1493 LOGICAL :: stat, GotIt, Compressible 1494 1495 REAL(KIND=dp) :: Grad(3,3), Grad1(3,3), Stress(3,3), Normal(3), ForceSolved(3), & 1496 EdgeLength,u, v, w, s, detJ 1497 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 1498 1499 REAL(KIND=dp), ALLOCATABLE :: EdgeBasis(:), dEdgeBasisdx(:,:), Basis(:),dBasisdx(:,:) 1500 REAL(KIND=dp), ALLOCATABLE :: x(:), y(:), z(:), ExtPressure(:) 1501 REAL(KIND=dp), ALLOCATABLE :: Temperature(:), Tension(:), SlipCoeff(:,:) 1502 REAL(KIND=dp), ALLOCATABLE :: Velocity(:,:), Pressure(:), Force(:,:), NodalViscosity(:) 1503 1504 REAL(KIND=dp) :: Residual(3), ResidualNorm, Viscosity, Slip, Dir(3) 1505 1506 TYPE(Variable_t), POINTER :: TempSol 1507 TYPE(ValueList_t), POINTER :: Material 1508 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1509!------------------------------------------------------------------------------ 1510 1511! Initialize: 1512! ----------- 1513 Indicator = 0.0d0 1514 Gnorm = 0.0d0 1515 1516 Metric = 0.0d0 1517 DO i=1,3 1518 Metric(i,i) = 1.0d0 1519 END DO 1520 1521 SELECT CASE( CurrentCoordinateSystem() ) 1522 CASE( AxisSymmetric, CylindricSymmetric ) 1523 DIM = 3 1524 CASE DEFAULT 1525 DIM = CoordinateSystemDimension() 1526 END SELECT 1527 1528 DOFs = DIM + 1 1529 IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs-1 1530! 1531! -------------------------------------------------- 1532 Element => Edge % BoundaryInfo % Left 1533 1534 IF ( .NOT. ASSOCIATED( Element ) ) THEN 1535 1536 Element => Edge % BoundaryInfo % Right 1537 1538 ELSE IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) THEN 1539 1540 Element => Edge % BoundaryInfo % Right 1541 1542 END IF 1543 1544 IF ( .NOT. ASSOCIATED( Element ) ) RETURN 1545 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 1546 1547 En = Edge % TYPE % NumberOfNodes 1548 Pn = Element % TYPE % NumberOfNodes 1549 1550 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 1551 1552 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 1553 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 1554 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 1555 1556 ALLOCATE( Nodes % x(Pn), Nodes % y(Pn), Nodes % z(Pn) ) 1557 1558 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 1559 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 1560 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 1561 1562 ALLOCATE( EdgeBasis(En), dEdgeBasisdx(En,3), Basis(Pn), dBasisdx(Pn,3), & 1563 x(En), y(En), z(En), ExtPressure(En), Temperature(Pn), Tension(En), & 1564 SlipCoeff(3,En), Velocity(3,Pn), Pressure(Pn), Force(3,En), NodalViscosity(En) ) 1565 1566 DO l = 1,En 1567 DO k = 1,Pn 1568 IF ( Edge % NodeIndexes(l) == Element % NodeIndexes(k) ) THEN 1569 x(l) = Element % TYPE % NodeU(k) 1570 y(l) = Element % TYPE % NodeV(k) 1571 z(l) = Element % TYPE % NodeW(k) 1572 EXIT 1573 END IF 1574 END DO 1575 END DO 1576! 1577! Integrate square of residual over boundary element: 1578! --------------------------------------------------- 1579 1580 Indicator = 0.0d0 1581 EdgeLength = 0.0d0 1582 ResidualNorm = 0.0d0 1583 1584 DO bc=1,Model % NumberOfBCs 1585 IF ( Edge % BoundaryInfo % Constraint /= Model % BCs(bc) % Tag ) CYCLE 1586 1587! IF ( .NOT. ListGetLogical( Model % BCs(bc) % Values, & 1588! 'Flow Force BC', gotIt ) ) CYCLE 1589! 1590! Get material parameters: 1591! ------------------------ 1592 1593 k = ListGetInteger(Model % Bodies(Element % BodyId) % Values,'Material', & 1594 minv=1, maxv=Model % NumberOfMaterials ) 1595 Material => Model % Materials(k) % Values 1596 1597 NodalViscosity(1:En) = ListGetReal( Material, & 1598 'Viscosity', En, Edge % NodeIndexes, GotIt ) 1599 1600 Compressible = .FALSE. 1601 IF ( ListGetString( Material, 'Compressibility Model', GotIt ) == & 1602 'perfect gas equation 1' ) Compressible = .TRUE. 1603 1604 1605! Given traction: 1606! --------------- 1607 Force = 0.0d0 1608 1609 Force(1,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1610 'Pressure 1', En, Edge % NodeIndexes, GotIt ) 1611 1612 Force(2,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1613 'Pressure 2', En, Edge % NodeIndexes, GotIt ) 1614 1615 Force(3,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1616 'Pressure 3', En, Edge % NodeIndexes, GotIt ) 1617 1618! 1619! Force in normal direction: 1620! --------------------------- 1621 ExtPressure(1:En) = ListGetReal( Model % BCs(bc) % Values, & 1622 'External Pressure', En, Edge % NodeIndexes, GotIt ) 1623 1624! 1625! Slip BC condition: 1626! ------------------ 1627 SlipCoeff = 0.0d0 1628 SlipCoeff(1,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1629 'Slip Coefficient 1',En,Edge % NodeIndexes,GotIt ) 1630 1631 SlipCoeff(2,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1632 'Slip Coefficient 2',En,Edge % NodeIndexes,GotIt ) 1633 1634 SlipCoeff(3,1:En) = ListGetReal( Model % BCs(bc) % Values, & 1635 'Slip Coefficient 3',En,Edge % NodeIndexes,GotIt ) 1636 1637! 1638! Surface tension induced by temperature gradient (or otherwise): 1639! --------------------------------------------------------------- 1640 TempSol => VariableGet( Mesh % Variables, 'Temperature', .TRUE. ) 1641 1642 IF ( ASSOCIATED( TempSol ) ) THEN 1643 Tension(1:En) = ListGetReal( Model % BCs(bc) % Values, & 1644 'Surface Tension Expansion Coefficient',En,Edge % NodeIndexes,gotIt ) 1645 1646 IF ( gotIt ) THEN 1647 DO n=1,En 1648 k = TempSol % Perm( Edge % NodeIndexes(n) ) 1649 IF (k>0) Tension(n) = 1.0d0 - Tension(n) * TempSol % Values(k) 1650 END DO 1651 1652 Tension(1:En) = Tension(1:En) * ListGetReal( & 1653 Model % BCs(bc) % Values,'Surface Tension Coefficient', & 1654 En, Edge % NodeIndexes ) 1655 ELSE 1656 Tension(1:En) = ListGetReal( & 1657 Model % BCs(bc) % Values,'Surface Tension Coefficient', & 1658 En, Edge % NodeIndexes,gotIt ) 1659 END IF 1660 ELSE 1661 Tension(1:En) = ListGetReal( & 1662 Model % BCs(bc) % Values,'Surface Tension Coefficient', & 1663 En, Edge % NodeIndexes,gotIt ) 1664 END IF 1665 1666! 1667! If dirichlet BC for velocity in any direction given, 1668! nullify force in that directon: 1669! ------------------------------------------------------------------ 1670 Dir = 1 1671 s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 1', GotIt ) 1672 IF ( GotIt ) Dir(1) = 0 1673 1674 s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 2', GotIt ) 1675 IF ( GotIt ) Dir(2) = 0 1676 1677 s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 3', GotIt ) 1678 IF ( GotIt ) Dir(3) = 0 1679 1680! 1681! Elementwise nodal solution: 1682! --------------------------- 1683 Velocity = 0.0d0 1684 DO k=1,DOFs-1 1685 Velocity(k,1:Pn) = Quant(DOFs*Perm(Element % NodeIndexes)-DOFs + k) 1686 END DO 1687 Pressure(1:Pn) = Quant( DOFs*Perm(Element % NodeIndexes) ) 1688 1689! do the integration: 1690! ------------------- 1691 EdgeLength = 0.0d0 1692 ResidualNorm = 0.0d0 1693 1694 IntegStuff = GaussPoints( Edge ) 1695 1696 DO t=1,IntegStuff % n 1697 u = IntegStuff % u(t) 1698 v = IntegStuff % v(t) 1699 w = IntegStuff % w(t) 1700 1701 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 1702 EdgeBasis, dEdgeBasisdx ) 1703 1704 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1705 s = IntegStuff % s(t) * detJ 1706 ELSE 1707 u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 1708 v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 1709 w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 1710 1711 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 1712 Symb, dSymb, u, v, w ) 1713 1714 s = IntegStuff % s(t) * detJ * SqrtMetric 1715 END IF 1716 1717 Normal = NormalVector( Edge, EdgeNodes, u, v, .TRUE. ) 1718 1719 u = SUM( EdgeBasis(1:En) * x(1:En) ) 1720 v = SUM( EdgeBasis(1:En) * y(1:En) ) 1721 w = SUM( EdgeBasis(1:En) * z(1:En) ) 1722 1723 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 1724 Basis, dBasisdx ) 1725 1726 Viscosity = SUM( NodalViscosity(1:En) * EdgeBasis(1:En) ) 1727 1728 Residual = 0.0d0 1729! 1730! Given force at the integration point: 1731! ------------------------------------- 1732 Residual = Residual + MATMUL( Force(:,1:En), EdgeBasis(1:En) ) - & 1733 SUM( ExtPressure(1:En) * EdgeBasis(1:En) ) * Normal 1734 1735! 1736! Slip velocity BC: 1737! ----------------- 1738 DO i=1,DIM 1739 Slip = SUM( SlipCoeff(i,1:En) * EdgeBasis(1:En) ) 1740 Residual(i) = Residual(i) - & 1741 Slip * SUM( Velocity(i,1:Pn) * Basis(1:Pn) ) 1742 END DO 1743 1744! 1745! Tangential tension force: 1746! ------------------------- 1747 DO i=1,DIM 1748 Residual(i) = Residual(i) + & 1749 SUM( dEdgeBasisdx(1:En,i) * Tension(1:En) ) 1750 END DO 1751 1752! 1753! Force given by the computed solution: 1754! ------------------------------------- 1755! 1756! Stress tensor on the boundary: 1757! ------------------------------ 1758 Grad = MATMUL( Velocity(:,1:Pn), dBasisdx(1:Pn,:) ) 1759 1760 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 1761 Grad1 = Grad 1762 DO i=1,DIM 1763 DO k=1,DIM 1764 DO l=1,DIM 1765 Grad1(i,k) = Grad1(i,k) - & 1766 Symb(k,l,i) * SUM ( Velocity(l,1:Pn) * Basis(1:Pn) ) 1767 END DO 1768 END DO 1769 END DO 1770 1771 Grad = 0.0d0 1772 DO i=1,DIM 1773 DO k=1,DIM 1774 DO l=1,DIM 1775 Grad(i,k) = Grad(i,k) + Metric(k,l) * Grad1(i,l) 1776 END DO 1777 END DO 1778 END DO 1779 END IF 1780 1781 Stress = Viscosity * ( Grad + TRANSPOSE(Grad) ) 1782 Stress = Stress - Metric * SUM( Pressure(1:Pn) * Basis(1:Pn) ) 1783 1784 IF ( Compressible ) THEN 1785 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1786 DO i=1,DIM 1787 DO k=1,DIM 1788 Stress(i,i) = Stress(i,i) - & 1789 (2.0d0/3.0d0) * Viscosity * Grad(k,k) 1790 END DO 1791 END DO 1792 ELSE 1793 DO i=1,DIM 1794 DO k=1,DIM 1795 DO l=1,DIM 1796 Stress(i,k) = Stress(i,k) - & 1797 Metric(i,k) * (2.0d0/3.0d0) * Viscosity * Grad(l,l) 1798 END DO 1799 END DO 1800 END DO 1801 END IF 1802 END IF 1803 1804 ForceSolved = MATMUL(Stress,Normal) 1805 Residual = Residual - ForceSolved * Dir 1806 1807 EdgeLength = EdgeLength + s 1808 1809 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1810 Gnorm = Gnorm + s * SUM( ForceSolved**2 ) 1811 ResidualNorm = ResidualNorm + s * SUM( Residual(1:DIM) ** 2 ) 1812 ELSE 1813 CALL InvertMatrix( Metric,3 ) 1814 DO i=1,DIM 1815 DO k=1,DIM 1816 ResidualNorm = ResidualNorm + & 1817 s * Metric(i,k) * Residual(i) * Residual(k) 1818 Gnorm = GNorm + s * Metric(i,k) * & 1819 ForceSolved(i) * ForceSolved(k) 1820 END DO 1821 END DO 1822 END IF 1823 END DO 1824 EXIT 1825 END DO 1826 1827 IF ( CoordinateSystemDimension() == 3 ) EdgeLength = SQRT(EdgeLength) 1828 Indicator = EdgeLength * ResidualNorm 1829 1830 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z) 1831 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 1832 1833 DEALLOCATE( EdgeBasis, dEdgeBasisdx, Basis, dBasisdx, x, y, z, & 1834 ExtPressure, Temperature, Tension,SlipCoeff, Velocity, Pressure, & 1835 Force, NodalViscosity ) 1836!------------------------------------------------------------------------------ 1837 END FUNCTION FlowBoundaryResidual 1838!------------------------------------------------------------------------------ 1839 1840 1841!------------------------------------------------------------------------------ 1842!> Compute the residual of the Navier-Stokes equation for the edge elements. 1843!------------------------------------------------------------------------------ 1844 FUNCTION FlowEdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT( Indicator ) 1845!------------------------------------------------------------------------------ 1846 USE DefUtils 1847 IMPLICIT NONE 1848 1849 TYPE(Model_t) :: Model 1850 INTEGER :: Perm(:) 1851 REAL(KIND=dp) :: Quant(:), Indicator(2) 1852 TYPE( Mesh_t ), POINTER :: Mesh 1853 TYPE( Element_t ), POINTER :: Edge 1854!------------------------------------------------------------------------------ 1855 1856 TYPE(Nodes_t) :: Nodes, EdgeNodes 1857 TYPE(Element_t), POINTER :: Element 1858 1859 INTEGER :: i,j,k,l,n,t,DIM,DOFs,En,Pn 1860 LOGICAL :: stat, GotIt 1861 1862 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 1863 REAL(KIND=dp) :: Stress(3,3,2), Jump(3), Viscosity 1864 1865 REAL(KIND=dp) :: Grad(3,3), Grad1(3,3), Normal(3) 1866 1867 REAL(KIND=dp), ALLOCATABLE :: NodalViscosity(:), x(:), y(:), z(:), & 1868 EdgeBasis(:), Basis(:), dBasisdx(:,:) 1869 REAL(KIND=dp), ALLOCATABLE :: Velocity(:,:), Pressure(:) 1870 1871 REAL(KIND=dp) :: ResidualNorm, EdgeLength, u, v, w, s, detJ 1872 1873 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1874!------------------------------------------------------------------------------ 1875 1876! Initialize: 1877! ----------- 1878 SELECT CASE( CurrentCoordinateSystem() ) 1879 CASE( AxisSymmetric, CylindricSymmetric ) 1880 DIM = 3 1881 CASE DEFAULT 1882 DIM = CoordinateSystemDimension() 1883 END SELECT 1884 1885 DOFs = DIM + 1 1886 IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs - 1 1887 1888 Metric = 0.0d0 1889 DO i = 1,3 1890 Metric(i,i) = 1.0d0 1891 END DO 1892 1893 Grad = 0.0d0 1894! 1895! --------------------------------------------- 1896 1897 Element => Edge % BoundaryInfo % Left 1898 n = Element % TYPE % NumberOfNodes 1899 1900 Element => Edge % BoundaryInfo % Right 1901 n = MAX( n, Element % TYPE % NumberOfNodes ) 1902 1903 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 1904 1905 En = Edge % TYPE % NumberOfNodes 1906 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 1907 1908 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 1909 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 1910 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 1911 1912 ALLOCATE( NodalViscosity(En), x(En), y(En), z(En), EdgeBasis(En), & 1913 Basis(n), dBasisdx(n,3), Velocity(3,n), Pressure(n) ) 1914 1915! Integrate square of jump over edge: 1916! ------------------------------------ 1917 ResidualNorm = 0.0d0 1918 EdgeLength = 0.0d0 1919 Indicator = 0.0d0 1920 1921 IntegStuff = GaussPoints( Edge ) 1922 1923 DO t=1,IntegStuff % n 1924 1925 u = IntegStuff % u(t) 1926 v = IntegStuff % v(t) 1927 w = IntegStuff % w(t) 1928 1929 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 1930 EdgeBasis, dBasisdx ) 1931 1932 Normal = NormalVector( Edge, EdgeNodes, u, v, .FALSE. ) 1933 1934 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1935 s = IntegStuff % s(t) * detJ 1936 ELSE 1937 u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 1938 v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 1939 w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 1940 1941 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 1942 Symb, dSymb, u, v, w ) 1943 s = IntegStuff % s(t) * detJ * SqrtMetric 1944 END IF 1945 1946 Stress = 0.0d0 1947 DO i = 1,2 1948 IF ( i==1 ) THEN 1949 Element => Edge % BoundaryInfo % Left 1950 ELSE 1951 Element => Edge % BoundaryInfo % Right 1952 END IF 1953 1954 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) CYCLE 1955 1956 Pn = Element % TYPE % NumberOfNodes 1957 Nodes % x(1:Pn) = Mesh % Nodes % x(Element % NodeIndexes) 1958 Nodes % y(1:Pn) = Mesh % Nodes % y(Element % NodeIndexes) 1959 Nodes % z(1:Pn) = Mesh % Nodes % z(Element % NodeIndexes) 1960 1961 DO j = 1,En 1962 DO k = 1,Pn 1963 IF ( Edge % NodeIndexes(j) == Element % NodeIndexes(k) ) THEN 1964 x(j) = Element % TYPE % NodeU(k) 1965 y(j) = Element % TYPE % NodeV(k) 1966 z(j) = Element % TYPE % NodeW(k) 1967 EXIT 1968 END IF 1969 END DO 1970 END DO 1971 1972 u = SUM( EdgeBasis(1:En) * x(1:En) ) 1973 v = SUM( EdgeBasis(1:En) * y(1:En) ) 1974 w = SUM( EdgeBasis(1:En) * z(1:En) ) 1975 1976 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 1977 Basis, dBasisdx ) 1978 1979 k = ListGetInteger( Model % Bodies( Element % BodyId) % Values, 'Material', & 1980 minv=1, maxv=Model % NumberOfMaterials ) 1981 1982 NodalViscosity(1:En) = ListGetReal( & 1983 Model % Materials(k) % Values, 'Viscosity', & 1984 En, Edge % NodeIndexes, GotIt ) 1985 1986 Viscosity = SUM( NodalViscosity(1:En) * EdgeBasis(1:En) ) 1987! 1988! Elementwise nodal solution: 1989! --------------------------- 1990 Velocity = 0.0d0 1991 DO k=1,DOFs-1 1992 Velocity(k,1:Pn) = Quant(DOFs*Perm(Element % NodeIndexes)-DOFs+k) 1993 END DO 1994 Pressure(1:Pn) = Quant( DOFs*Perm(Element % NodeIndexes) ) 1995! 1996! Stress tensor on the edge: 1997! -------------------------- 1998 Grad = MATMUL( Velocity(:,1:Pn), dBasisdx(1:Pn,:) ) 1999 2000 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2001 Grad1 = Grad 2002 DO j=1,DIM 2003 DO k=1,DIM 2004 DO l=1,DIM 2005 Grad1(j,k) = Grad1(j,k) - & 2006 Symb(k,l,j) * SUM ( Velocity(l,1:Pn) * Basis(1:Pn) ) 2007 END DO 2008 END DO 2009 END DO 2010 2011 Grad = 0.0d0 2012 DO j=1,DIM 2013 DO k=1,DIM 2014 DO l=1,DIM 2015 Grad(j,k) = Grad(j,k) + Metric(k,l) * Grad1(j,l) 2016 END DO 2017 END DO 2018 END DO 2019 END IF 2020 2021 Stress(:,:,i) = Viscosity * ( Grad + TRANSPOSE(Grad) ) 2022 2023 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2024 DO j=1,DIM 2025 Stress(j,j,i) = Stress(j,j,i) - SUM( Pressure(1:Pn) * Basis(1:Pn)) 2026 DO k=1,DIM 2027 Stress(j,j,i) = Stress(j,j,i) - (2.0d0/3.0d0)*Viscosity*Grad(k,k) 2028 END DO 2029 END DO 2030 ELSE 2031 DO j=1,DIM 2032 DO k=1,DIM 2033 Stress(j,k,i) = Stress(j,k,i) - & 2034 Metric(j,k) * SUM( Pressure(1:Pn) * Basis(1:Pn) ) 2035 2036 DO l=1,DIM 2037 Stress(j,k,i) = Stress(j,k,i) - & 2038 Metric(j,k) * (2.0d0/3.0d0) * Viscosity * Grad(l,l) 2039 END DO 2040 END DO 2041 END DO 2042 END IF 2043 2044 END DO 2045 2046 EdgeLength = EdgeLength + s 2047 2048 Jump = MATMUL( ( Stress(:,:,1) - Stress(:,:,2)), Normal ) 2049 2050 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2051 ResidualNorm = ResidualNorm + s * SUM( Jump(1:DIM) ** 2 ) 2052 ELSE 2053 CALL InvertMatrix( Metric,3 ) 2054 DO i=1,DIM 2055 DO j=1,DIM 2056 ResidualNorm = ResidualNorm + s*Metric(i,j)*Jump(i)*Jump(j) 2057 END DO 2058 END DO 2059 END IF 2060 END DO 2061 2062 Indicator = EdgeLength * ResidualNorm 2063 2064 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z) 2065 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 2066 2067 DEALLOCATE( NodalViscosity, x, y, z, EdgeBasis, & 2068 Basis, dBasisdx, Velocity, Pressure ) 2069!------------------------------------------------------------------------------ 2070 END FUNCTION FlowEdgeResidual 2071!------------------------------------------------------------------------------ 2072 2073 2074!------------------------------------------------------------------------------ 2075!> Compute the residual of the Navier-Stokes equation for the bulk elements. 2076!------------------------------------------------------------------------------ 2077 FUNCTION FlowInsideResidual( Model, Element, & 2078 Mesh, Quant, Perm, Fnorm ) RESULT( Indicator ) 2079!------------------------------------------------------------------------------ 2080 USE DefUtils 2081!------------------------------------------------------------------------------ 2082 IMPLICIT NONE 2083!------------------------------------------------------------------------------ 2084 TYPE(Model_t) :: Model 2085 INTEGER :: Perm(:) 2086 REAL(KIND=dp) :: Quant(:), Indicator(2), FNorm 2087 TYPE( Mesh_t ), POINTER :: Mesh 2088 TYPE( Element_t ), POINTER :: Element 2089!------------------------------------------------------------------------------ 2090 2091 TYPE(Nodes_t) :: Nodes 2092 2093 INTEGER :: i,j,k,l,m,n,t,DIM,DOFs 2094 2095 LOGICAL :: stat, GotIt, Compressible, Convect 2096 2097 TYPE( Variable_t ), POINTER :: Var 2098 2099 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 2100 REAL(KIND=dp) :: Density, Viscosity,u, v, w, s, detJ 2101 REAL(KIND=dp) :: Residual(4), ResidualNorm, Area, ReferencePressure, dt 2102 2103 REAL(KIND=dp), ALLOCATABLE :: NodalViscosity(:), NodalDensity(:), & 2104 Basis(:), dBasisdx(:,:), ddBasisddx(:,:,:) 2105 REAL(KIND=dp),ALLOCATABLE :: Velocity(:,:), Pressure(:) 2106 REAL(KIND=dp),ALLOCATABLE :: PrevVelo(:,:), PrevPres(:) 2107 REAL(KIND=dp),ALLOCATABLE :: Temperature(:), NodalForce(:,:) 2108 REAL(KIND=dp),ALLOCATABLE :: HeatCapacity(:), ReferenceTemperature(:), & 2109 HeatExpansionCoeff(:) 2110 2111 REAL(KIND=dp) :: SpecificHeatRatio 2112 2113 REAL(KIND=dp), POINTER :: Gravity(:,:) 2114 2115 TYPE(ValueList_t), POINTER :: Material 2116 2117 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 2118!------------------------------------------------------------------------------ 2119 2120! Initialize: 2121! ----------- 2122 Indicator = 0.0d0 2123 FNorm = 0.0d0 2124 2125 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 2126 2127 Metric = 0.0d0 2128 DO i=1,3 2129 Metric(i,i) = 1.0d0 2130 END DO 2131 2132 SELECT CASE( CurrentCoordinateSystem() ) 2133 CASE( AxisSymmetric, CylindricSymmetric ) 2134 DIM = 3 2135 CASE DEFAULT 2136 DIM = CoordinateSystemDimension() 2137 END SELECT 2138 2139 DOFs = DIM + 1 2140 IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs-1 2141! 2142! Element nodal points: 2143! --------------------- 2144 n = Element % TYPE % NumberOfNodes 2145 2146 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 2147 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 2148 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 2149 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 2150 2151 ALLOCATE( NodalViscosity(n), NodalDensity(n), Basis(n), dBasisdx(n,3), & 2152 ddBasisddx(n,3,3), Velocity(3,n), Pressure(n), PrevVelo(3,n), & 2153 PrevPres(n), Temperature(n), NodalForce(4,n), HeatCapacity(n), & 2154 ReferenceTemperature(n), HeatExpansionCoeff(n) ) 2155 2156! 2157! Material parameters: density, viscosity, etc. 2158! ---------------------------------------------- 2159 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', & 2160 minv=1, maxv=Model % NumberOfMaterials ) 2161 2162 Material => Model % Materials(k) % Values 2163 2164 NodalDensity(1:n) = ListGetReal( & 2165 Material, 'Density', n, Element % NodeIndexes, GotIt ) 2166 2167 NodalViscosity(1:n) = ListGetReal( & 2168 Material, 'Viscosity', n, Element % NodeIndexes, GotIt ) 2169 2170 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values,'Equation', & 2171 minv=1, maxv=Model % NumberOfEquations ) 2172 2173 Convect = ListGetLogical( Model % Equations(k) % Values, & 2174 'NS Convect', GotIt ) 2175 IF ( .NOT. GotIt ) Convect = .TRUE. 2176! 2177! Elementwise nodal solution: 2178! --------------------------- 2179 Velocity = 0.0d0 2180 DO k=1,DOFs-1 2181 Velocity(k,1:n) = Quant( DOFs*Perm(Element % NodeIndexes)-DOFs+k ) 2182 END DO 2183 Pressure(1:n) = Quant( DOFs*Perm(Element % NodeIndexes) ) 2184 2185! 2186! Check for time dep. 2187! ------------------- 2188 PrevPres(1:n) = Pressure(1:n) 2189 PrevVelo(1:3,1:n) = Velocity(1:3,1:n) 2190 2191 dt = Model % Solver % dt 2192 2193 IF ( ListGetString( Model % Simulation, 'Simulation Type') == 'transient' ) THEN 2194 Var => VariableGet( Model % Variables, 'Flow Solution', .TRUE. ) 2195 2196 PrevVelo = 0.0d0 2197 DO k=1,DOFs-1 2198 PrevVelo(k,1:n) = & 2199 Var % PrevValues(DOFs*Var % Perm(Element % NodeIndexes)-DOFs+k,1) 2200 END DO 2201 PrevPres(1:n)=Var % PrevValues(DOFs*Var % Perm(Element % NodeIndexes),1) 2202 END IF 2203 2204 2205! 2206! Check for compressible flow equations: 2207! -------------------------------------- 2208 Compressible = .FALSE. 2209 2210 IF ( ListGetString( Material, 'Compressibility Model', GotIt ) == & 2211 'perfect gas equation 1' ) THEN 2212 2213 Compressible = .TRUE. 2214 2215 Var => VariableGet( Mesh % Variables, 'Temperature', .TRUE. ) 2216 IF ( ASSOCIATED( Var ) ) THEN 2217 Temperature(1:n) = & 2218 Var % Values( Var % Perm(Element % NodeIndexes) ) 2219 ELSE 2220 Temperature(1:n) = ListGetReal( Material, & 2221 'Reference Temperature',n,Element % NodeIndexes ) 2222 END IF 2223 2224 SpecificHeatRatio = ListGetConstReal( Material, & 2225 'Specific Heat Ratio' ) 2226 2227 ReferencePressure = ListGetConstReal( Material, & 2228 'Reference Pressure' ) 2229 2230 HeatCapacity(1:n) = ListGetReal( Material, & 2231 'Heat Capacity',n,Element % NodeIndexes ) 2232 2233 NodalDensity(1:n) = (Pressure(1:n) + ReferencePressure) * SpecificHeatRatio / & 2234 ( (SpecificHeatRatio - 1) * HeatCapacity(1:n) * Temperature(1:n) ) 2235 END IF 2236! 2237! Body Forces: 2238! ------------ 2239! 2240 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, & 2241 'Body Force', GotIt, 1, Model % NumberOfBodyForces ) 2242 2243 NodalForce = 0.0d0 2244 2245 IF ( GotIt .AND. k > 0 ) THEN 2246! 2247! Boussinesq approximation of heat expansion for 2248! incompressible flow equations: 2249! 2250! Density for the force term equals to 2251! 2252! \rho = rho_0 (1-\beta(T-T_0)), 2253! 2254! where \beta is the heat expansion coefficient, 2255! T temperature and \rho_0 and T_0 correspond to 2256! stress free state. Otherwise density is assumed 2257! constant. 2258! ---------------------------------------------- 2259 IF (ListGetLogical(Model % BodyForces(k) % Values,'Boussinesq',GotIt)) THEN 2260 2261 Var => VariableGet( Mesh % Variables, 'Temperature', .TRUE. ) 2262 IF ( ASSOCIATED( Var ) ) THEN 2263 Temperature(1:n) = & 2264 Var % Values( Var % Perm(Element % NodeIndexes) ) 2265 2266 HeatExpansionCoeff(1:n) = ListGetReal( Material, & 2267 'Heat Expansion Coefficient',n,Element % NodeIndexes ) 2268 2269 ReferenceTemperature(1:n) = ListGetReal( Material, & 2270 'Reference Temperature',n,Element % NodeIndexes ) 2271 2272 Gravity => ListGetConstRealArray( Model % Constants, & 2273 'Gravity' ) 2274 2275 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values,'Equation', & 2276 minv=1, maxv=Model % NumberOfEquations ) 2277 2278 IF ( ListGetLogical( Model % Equations(k) % Values, & 2279 'Hydrostatic Pressure', GotIt) ) THEN 2280 DO i=1,DIM 2281 NodalForce(i,1:n) = ( 1 - HeatExpansionCoeff(1:n) * & 2282 ( Temperature(1:n) - ReferenceTemperature(1:n) ) ) * & 2283 Gravity(i,1) * Gravity(4,1) 2284 END DO 2285 ELSE 2286 DO i=1,DIM 2287 NodalForce(i,1:n) = ( -HeatExpansionCoeff(1:n) * & 2288 ( Temperature(1:n) - ReferenceTemperature(1:n) ) ) * & 2289 Gravity(i,1) * Gravity(4,1) 2290 END DO 2291 END IF 2292 END IF 2293 END IF 2294 2295! 2296! Given external force: 2297! --------------------- 2298 NodalForce(1,1:n) = NodalForce(1,1:n) + ListGetReal( & 2299 Model % BodyForces(k) % Values, 'Flow BodyForce 1', & 2300 n, Element % NodeIndexes, GotIt ) 2301 2302 NodalForce(2,1:n) = NodalForce(2,1:n) + ListGetReal( & 2303 Model % BodyForces(k) % Values, 'Flow BodyForce 2', & 2304 n, Element % NodeIndexes, GotIt ) 2305 2306 NodalForce(3,1:n) = NodalForce(3,1:n) + ListGetReal( & 2307 Model % BodyForces(k) % Values, 'Flow BodyForce 3', & 2308 n, Element % NodeIndexes, GotIt ) 2309 END IF 2310! 2311! Integrate square of residual over element: 2312! ------------------------------------------ 2313 ResidualNorm = 0.0d0 2314 Area = 0.0d0 2315 2316 IntegStuff = GaussPoints( Element ) 2317 2318 DO t=1,IntegStuff % n 2319 u = IntegStuff % u(t) 2320 v = IntegStuff % v(t) 2321 w = IntegStuff % w(t) 2322 2323 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 2324 Basis, dBasisdx, ddBasisddx, .TRUE. ) 2325 2326 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2327 s = IntegStuff % s(t) * detJ 2328 ELSE 2329 u = SUM( Basis(1:n) * Nodes % x(1:n) ) 2330 v = SUM( Basis(1:n) * Nodes % y(1:n) ) 2331 w = SUM( Basis(1:n) * Nodes % z(1:n) ) 2332 2333 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 2334 Symb, dSymb, u, v, w ) 2335 s = IntegStuff % s(t) * detJ * SqrtMetric 2336 END IF 2337 2338 Density = SUM( NodalDensity(1:n) * Basis(1:n) ) 2339 Viscosity = SUM( NodalViscosity(1:n) * Basis(1:n) ) 2340! 2341! Residual of the navier-stokes equations: 2342! 2343! or more generally: 2344! 2345! ---------------------------------------------------------- 2346! 2347 Residual = 0.0d0 2348 DO i=1,DIM 2349! 2350! given force: 2351! ------------- 2352 Residual(i) = -Density * SUM( NodalForce(i,1:n) * Basis(1:n) ) 2353 2354 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2355! + grad(p): 2356! ---------- 2357 Residual(i) = Residual(i) + SUM( Pressure(1:n) * dBasisdx(1:n,i) ) 2358 2359 DO j=1,DIM 2360! 2361! - 2 ( \mu \epsilon^{ij} )_{,j}: 2362! ------------------------------- 2363 Residual(i) = Residual(i) - Viscosity * & 2364 SUM( Velocity(i,1:n) * ddBasisddx(1:n,j,j) ) 2365 2366 Residual(i) = Residual(i) - & 2367 SUM( NodalViscosity(1:n) * dBasisdx(1:n,j) ) * & 2368 SUM( Velocity(i,1:n) * dBasisdx(1:n,j) ) 2369 2370 Residual(i) = Residual(i) - Viscosity * & 2371 SUM( Velocity(j,1:n) * ddBasisddx(1:n,i,j) ) 2372 2373 Residual(i) = Residual(i) - & 2374 SUM( NodalViscosity(1:n) * dBasisdx(1:n,j) ) * & 2375 SUM( Velocity(j,1:n) * dBasisdx(1:n,i) ) 2376 2377 IF ( Compressible ) THEN 2378! 2379! + (2/3) grad(\mu div(u)): 2380! ------------------------- 2381 Residual(i) = Residual(i) + & 2382 Viscosity * ( 2.0d0 / 3.0d0 ) * & 2383 SUM( Velocity(j,1:n) * ddBasisddx(1:n,j,i) ) 2384 2385 Residual(i) = Residual(i) + & 2386 SUM( NodalViscosity(1:n) * dBasisdx(1:n,i) ) * & 2387 SUM( Velocity(j,1:n) * dBasisdx(1:n,j) ) 2388 2389 END IF 2390 END DO 2391 2392 IF ( Convect ) THEN 2393! 2394! + \rho * (@u/@t + u.grad(u)): 2395! ----------------------------- 2396 Residual(i) = Residual(i) + Density * & 2397 SUM((Velocity(i,1:n)-PrevVelo(i,1:n))*Basis(1:n)) / dt 2398 2399 DO j=1,DIM 2400 Residual(i) = Residual(i) + & 2401 Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * & 2402 SUM( Velocity(i,1:n) * dBasisdx(1:n,j) ) 2403 END DO 2404 END IF 2405 ELSE 2406! + g^{ij}p_{,j}: 2407! --------------- 2408 DO j=1,DIM 2409 Residual(i) = Residual(i) + Metric(i,j) * & 2410 SUM( Pressure(1:n) * dBasisdx(1:n,i) ) 2411 END DO 2412 2413! - g^{jk} (\mu u^i_{,k})_{,j}): 2414! ------------------------------ 2415 DO j=1,DIM 2416 DO k=1,DIM 2417 Residual(i) = Residual(i) - & 2418 Metric(j,k) * Viscosity * & 2419 SUM( Velocity(i,1:n) * ddBasisddx(1:n,j,k) ) 2420 2421 DO l=1,DIM 2422 Residual(i) = Residual(i) + & 2423 Metric(j,k) * Viscosity * Symb(j,k,l) * & 2424 SUM( Velocity(i,1:n) * dBasisdx(1:n,l) ) 2425 2426 Residual(i) = Residual(i) - & 2427 Metric(j,k) * Viscosity * Symb(l,j,i) * & 2428 SUM( Velocity(l,1:n) * dBasisdx(1:n,k) ) 2429 2430 Residual(i) = Residual(i) - & 2431 Metric(j,k) * Viscosity * Symb(l,k,i) * & 2432 SUM( Velocity(l,1:n) * dBasisdx(1:n,j) ) 2433 2434 Residual(i) = Residual(i) - & 2435 Metric(j,k) * Viscosity * dSymb(l,j,i,k) * & 2436 SUM( Velocity(l,1:n) * Basis(1:n) ) 2437 2438 DO m=1,DIM 2439 Residual(i) = Residual(i) - Metric(j,k) * Viscosity *& 2440 Symb(m,k,i) * Symb(l,j,m) * & 2441 SUM( Velocity(l,1:n) * Basis(1:n) ) 2442 2443 Residual(i) = Residual(i) + Metric(j,k) * Viscosity *& 2444 Symb(j,k,m) * Symb(l,m,i) * & 2445 SUM( Velocity(l,1:n) * Basis(1:n) ) 2446 END DO 2447 END DO 2448 END DO 2449 END DO 2450 2451! - g^{ik} (\mu u^j_{,k})_{,j}): 2452! ------------------------------ 2453 DO j=1,DIM 2454 DO k=1,DIM 2455 Residual(i) = Residual(i) - & 2456 Metric(i,k) * Viscosity * & 2457 SUM( Velocity(j,1:n) * ddBasisddx(1:n,j,k) ) 2458 2459 DO l=1,DIM 2460 Residual(i) = Residual(i) + & 2461 Metric(i,k) * Viscosity * Symb(j,k,l) * & 2462 SUM( Velocity(j,1:n) * dBasisdx(1:n,l) ) 2463 2464 Residual(i) = Residual(i) - & 2465 Metric(i,k) * Viscosity * Symb(l,j,j) * & 2466 SUM( Velocity(l,1:n) * dBasisdx(1:n,k) ) 2467 2468 Residual(i) = Residual(i) - & 2469 Metric(i,k) * Viscosity * Symb(l,k,j) * & 2470 SUM( Velocity(l,1:n) * dBasisdx(1:n,j) ) 2471 2472 Residual(i) = Residual(i) - & 2473 Metric(i,k) * Viscosity * dSymb(l,j,j,k) * & 2474 SUM( Velocity(l,1:n) * Basis(1:n) ) 2475 2476 DO m=1,DIM 2477 Residual(i) = Residual(i) - Metric(i,k) * Viscosity *& 2478 Symb(m,k,j) * Symb(l,j,m) * & 2479 SUM( Velocity(l,1:n) * Basis(1:n) ) 2480 2481 Residual(i) = Residual(i) + Metric(i,k) * Viscosity *& 2482 Symb(j,k,m) * Symb(l,m,j) * & 2483 SUM( Velocity(l,1:n) * Basis(1:n) ) 2484 END DO 2485 END DO 2486 END DO 2487 END DO 2488 2489 IF ( Convect ) THEN 2490! 2491! + \rho * (@u/@t + u^j u^i_{,j}): 2492! -------------------------------- 2493 Residual(i) = Residual(i) + Density * & 2494 SUM((Velocity(i,1:n)-PrevVelo(i,1:n))*Basis(1:n)) / dt 2495 2496 DO j=1,DIM 2497 Residual(i) = Residual(i) + & 2498 Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * & 2499 SUM( Velocity(i,1:n) * dBasisdx(1:n,j) ) 2500 2501 DO k=1,DIM 2502 Residual(i) = Residual(i) + & 2503 Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * & 2504 Symb(j,k,i) * SUM( Velocity(k,1:n) * Basis(1:n) ) 2505 END DO 2506 END DO 2507 END IF 2508 END IF 2509 END DO 2510 2511! 2512! Continuity equation: 2513! -------------------- 2514 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2515! 2516! + \rho * div(u): 2517! ---------------- 2518 DO j=1,DIM 2519 Residual(DIM+1) = Residual(DIM+1) + & 2520 Density * SUM( Velocity(j,1:n) * dBasisdx(1:n,j) ) 2521 END DO 2522 2523 IF ( Compressible ) THEN 2524! 2525! + u.grad(\rho): 2526! ---------------- 2527 DO j=1,DIM 2528 Residual(DIM+1) = Residual(DIM+1) + & 2529 SUM( Velocity(j,1:n) * Basis(1:n) ) * & 2530 SUM( NodalDensity(1:n) * dBasisdx(1:n,j) ) 2531 END DO 2532 END IF 2533 ELSE 2534! 2535! + \rho * u^j_{,j}: 2536! ------------------ 2537 DO j=1,DIM 2538 Residual(DIM+1) = Residual(DIM+1) + & 2539 Density * SUM( Velocity(j,1:n) * dBasisdx(1:n,j) ) 2540 2541 DO k=1,DIM 2542 Residual(DIM+1) = Residual(DIM+1) + Density * & 2543 Symb(k,j,j) * SUM( Velocity(k,1:n) * Basis(1:n) ) 2544 END DO 2545 END DO 2546 2547 IF ( Compressible ) THEN 2548! 2549! + u^j \rho_{,j}: 2550! ---------------- 2551 DO j=1,DIM 2552 Residual(DIM+1) = Residual(DIM+1) + & 2553 SUM( Velocity(j,1:n) * Basis(1:n) ) * & 2554 SUM( NodalDensity(1:n) * dBasisdx(1:n,j) ) 2555 END DO 2556 END IF 2557 END IF 2558 2559 DO i=1,DIM 2560 FNorm = FNorm + s * (Density * SUM(NodalForce(i,1:n)*Basis(1:n))**2) 2561 END DO 2562 Area = Area + s 2563 2564 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 2565 ResidualNorm = ResidualNorm + & 2566 s * (Element % hK**2 * SUM(Residual(1:dim)**2) + Residual(dim+1)**2 ) 2567 ELSE 2568 CALL InvertMatrix( Metric,3 ) 2569 DO i=1,dim 2570 DO j=1,dim 2571 ResidualNorm = ResidualNorm + & 2572 s * Element % hK **2 * Metric(i,j) * Residual(i) * Residual(j) 2573 END DO 2574 END DO 2575 ResidualNorm = ResidualNorm + s * Residual(dim+1)**2 2576 END IF 2577 END DO 2578 2579! FNorm = Area * FNorm 2580 Indicator = ResidualNorm 2581 2582 DEALLOCATE( NodalViscosity, NodalDensity, Basis, dBasisdx, & 2583 ddBasisddx, Velocity, Pressure, PrevVelo, PrevPres, Temperature, & 2584 NodalForce, HeatCapacity, ReferenceTemperature, HeatExpansionCoeff,& 2585 Nodes % x, Nodes % y, Nodes % z ) 2586!------------------------------------------------------------------------------ 2587 END FUNCTION FlowInsideResidual 2588!------------------------------------------------------------------------------ 2589 2590!> \} 2591