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! * Modified by: Jussi Heikonen, Ville Savolainen, Peter Raback 35! * Modification date: 2.9.2006 36! * 37! *****************************************************************************/ 38 39 40!------------------------------------------------------------------------------ 41!> Initialization for the primary solver: PhaseChangeSolve 42!> \ingroup Solvers 43!------------------------------------------------------------------------------ 44SUBROUTINE PhaseChangeSolve_Init( Model,Solver,dt,TransientSimulation) 45!------------------------------------------------------------------------------ 46 USE DefUtils 47 IMPLICIT NONE 48!------------------------------------------------------------------------------ 49 TYPE(Model_t) :: Model 50 TYPE(Solver_t), TARGET :: Solver 51 LOGICAL :: TransientSimulation 52 REAL(KIND=dp) :: dt 53!------------------------------------------------------------------------------ 54 LOGICAL :: Found 55 TYPE(ValueList_t), POINTER :: Params 56 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName 57 58 Params => GetSolverParams() 59 VariableName = GetString(Params,'Variable') 60 CALL ListAddString( Params,NextFreeKeyword('Exported Variable ',Params), & 61 '-nooutput '//TRIM(ComponentName(VariableName))//'Move' ) 62 63 IF (ListGetLogical(Params,'Use Average Velocity',Found)) & 64 CALL ListAddString( Params,NextFreeKeyword('Exported Variable ',Params), & 65 '-nooutput '//TRIM(ComponentName(VariableName))//'MoveAve' ) 66!------------------------------------------------------------------------------ 67END SUBROUTINE PhaseChangeSolve_Init 68!------------------------------------------------------------------------------ 69 70 71!------------------------------------------------------------------------------ 72!> Solve the free surface in the phase change problem using Lagrangian techniques. 73!> \deprecated This had been replaced by separate versions for transient and steady state phase change. 74!> \ingroup Solvers 75!------------------------------------------------------------------------------ 76SUBROUTINE PhaseChangeSolve( Model,Solver,dt,TransientSimulation ) 77!------------------------------------------------------------------------------ 78 USE DefUtils 79 80 IMPLICIT NONE 81!------------------------------------------------------------------------------ 82 TYPE(Model_t) :: Model 83 TYPE(Solver_t), TARGET :: Solver 84 LOGICAL :: TransientSimulation 85 REAL(KIND=dp) :: dt 86!------------------------------------------------------------------------------ 87! Local variables 88!------------------------------------------------------------------------------ 89 TYPE(Element_t), POINTER :: CurrentElement, Parent, Element 90 TYPE(Variable_t), POINTER :: SurfSol, TempSol, HelpSol, HelpSol2 91 TYPE(Nodes_t) :: Nodes, PNodes 92 TYPE(GaussIntegrationPoints_t) :: IntegStuff 93 TYPE(Matrix_t),POINTER :: StiffMatrix 94 TYPE(ValueList_t), POINTER :: Material 95 TYPE(Solver_t), POINTER :: PSolver 96 97 REAL(KIND=dp) :: Normal(3), u, v, w, UPull(3), PrevUpull(3), & 98 Density, Update, MaxUpdate, MaxTempDiff, Relax, LocalRelax, AverageRelax, & 99 AverageRelax2, surf, xx, yy, r, detJ, Temp, MeltPoint, NonlinearTol, & 100 Temp1,Temp2,Temppi,dxmin,dymin, xmin, xmax, d, HeatCond, s, tave, prevtave, cvol, clim, ccum=1, & 101 tabs, prevtabs, volabs, prevvolabs, CoordMin(3), CoordMax(3), RelativeChange, & 102 dTdz, ttemp, NewtonAfterTol, area, volume, prevvolume, Eps, & 103 Norm, maxds, maxds0, ds, FluxCorrect = 1.0, dpos, & 104 pos0=0.0, prevpos0, ThermalInertiaFactor, Coeff, OrigNorm, & 105 SteadyDt, AveragingDt, Trip_Temp, MeanSurface, PullRelax, & 106 SpeedUp, PrevNorm 107 REAL(KIND=dp), POINTER :: Surface(:), Temperature(:), ForceVector(:), & 108 x(:), y(:), z(:), Basis(:), dBasisdx(:,:), NodalTemp(:), & 109 Conductivity(:), LatentHeat(:), TempDiff(:), Taverage(:), Taverage2(:), & 110 Normals(:), Weights(:), SurfaceMove(:), SurfaceMoveAve(:) 111 REAL (KIND=dp), ALLOCATABLE :: PrevTemp(:), IsoSurf(:,:) 112 113 REAL(KIND=dp), ALLOCATABLE :: & 114 LocalStiffMatrix(:,:), LocalForceVector(:), LocalMassMatrix(:,:) 115 116 INTEGER :: i,j,k,t,n,nn,pn,DIM,kl,kr,l, bc, Trip_node, NoBNodes, NonlinearIter, istat, & 117 NElems,ElementCode,Next,Vertex,ii,imin,NewtonAfterIter,Node, iter, LiquidInd, Visited = -1, & 118 SubroutineVisited = 0, NormalDirection, TangentDirection, CoordMini(3), CoordMaxi(3), & 119 Axis_node 120 INTEGER, POINTER :: NodeIndexes(:),TempPerm(:),SurfPerm(:),NormalsPerm(:) 121 122 LOGICAL :: Stat, FirstTime = .TRUE., Newton = .FALSE., UseTaverage, Debug, & 123 PullControl, PullVelocitySet = .FALSE., IsoSurfAllocated, AllocationsDone = .FALSE., & 124 TransientAlgo = .FALSE., SteadyAlgo = .FALSE., UseLoads, & 125 AverageNormal, SurfaceVelocitySet = .FALSE., TriplePointFixed 126 LOGICAL, POINTER :: NodeDone(:) 127 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName, str 128 129 SAVE FirstTime, Trip_node, NoBNodes, SubroutineVisited, prevpos0, & 130 PrevTemp, Newton, ccum, NormalDirection, TangentDirection, MeltPoint, & 131 ForceVector, FluxCorrect, NodeDone, Eps, PullControl, & 132 Visited, Nodes, PNodes, NodalTemp, Conductivity, LatentHeat, & 133 TempDiff, AllocationsDone, LocalStiffMatrix, LocalForceVector, LocalMassMatrix, & 134 x, y, z, Basis, dBasisdx, norm, Axis_node, PullVelocitySet, & 135 Normals, Weights, NormalsPerm, AverageNormal, SurfaceMove, SurfaceMoveAve, & 136 SurfaceVelocitySet, CoordMax, CoordMin, CoordMaxi, CoordMini, UPull 137 138 !------------------------------------------------------------------------------ 139 ! Decide which kind of algorithm to use for the current timestep size 140 !------------------------------------------------------------------------------ 141 142 TransientAlgo = TransientSimulation 143 IF(TransientAlgo) THEN 144 SteadyDt = ListGetConstReal(Solver % Values,'Steady Transition Timestep',Stat) 145 IF(Stat) TransientAlgo = (dt < SteadyDt) 146 END IF 147 SteadyAlgo = .NOT. TransientAlgo 148 149 CALL Info('PhaseChangeSolve', '--------------------------------------------') 150 IF(SteadyAlgo) CALL Info('PhaseChangeSolve', 'Using steady algorithm to find the isotherm') 151 IF(TransientAlgo) CALL Info('PhaseChangeSolve','Using transient algorithm for surface update') 152 CALL Info('PhaseChangeSolve', '--------------------------------------------') 153 154 SubroutineVisited = SubroutineVisited + 1 155 156 DIM = CoordinateSystemDimension() 157 IF(DIM /= 2) THEN 158 CALL Fatal('PhaseChangeSolve','Implemented only in 2D') 159 END IF 160 161!------------------------------------------------------------------------------ 162! Get variables needed for solution 163!------------------------------------------------------------------------------ 164 165 PSolver => Solver 166 SurfSol => Solver % Variable 167 Surface => SurfSol % Values 168 SurfPerm => SurfSol % Perm 169 IF(.NOT. ASSOCIATED (Surface) .OR. ALL(SurfPerm <= 0) ) THEN 170 CALL Fatal('PhaseChangeSolve','Surface field needed for Phase Change') 171 END IF 172 173 VariableName = ListGetString( Solver % Values, 'Phase Change Variable', Stat ) 174 IF(Stat) THEN 175 TempSol => VariableGet( Solver % Mesh % Variables, TRIM(VariableName) ) 176 ELSE 177 TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' ) 178 END IF 179 TempPerm => TempSol % Perm 180 Temperature => TempSol % Values 181 IF(.NOT. ASSOCIATED (Temperature) .OR. ALL(TempPerm <= 0) ) THEN 182 CALL Fatal('PhaseChangeSolve','Temperature field needed for Phase Change') 183 END IF 184 185 Relax = GetCReal( Solver % Values, & 186 'Nonlinear System Relaxation Factor', stat ) 187 IF ( .NOT. stat ) Relax = 1.0d0 188 NonlinearIter = ListGetInteger( Solver % Values, & 189 'Nonlinear System Max Iterations', stat ) 190 IF ( .NOT. stat ) NonlinearIter = 1 191 NonlinearTol = ListGetConstReal( Solver % Values, & 192 'Nonlinear System Convergence Tolerance', stat ) 193 194 PullControl = ListGetLogical( Solver % Values,'Pull Rate Control',stat) 195 TriplePointFixed = ListGetLogical( Solver % Values,'Triple Point Fixed',stat) 196 197!--------------------------------------------------------------------------------- 198! The first time the main axis of the free surface is determined 199! and some permanent vectors related to the surface are allocated. 200!--------------------------------------------------------------------------------- 201 202 IF(FirstTime) THEN 203 UPull = 0.0 204 NoBNodes = 0 205 CoordMax = -HUGE(CoordMax) 206 CoordMin = HUGE(CoordMin) 207 208 DO k=1, Model % Mesh % NumberOfNodes 209 IF( SurfPerm(k) <= 0) CYCLE 210 NoBnodes = NoBnodes + 1 211 212 DO j=1,DIM 213 IF(j==1) xx = Model % Mesh % Nodes % x(k) 214 IF(j==2) xx = Model % Mesh % Nodes % y(k) 215 IF(j==3) xx = Model % Mesh % Nodes % z(k) 216 IF(xx > CoordMax(j)) THEN 217 CoordMax(j) = xx 218 CoordMaxi(j) = k 219 END IF 220 IF(xx < CoordMin(j)) THEN 221 CoordMin(j) = xx 222 CoordMini(j) = k 223 END IF 224 END DO 225 END DO 226 227 ! Direction of minimum change 228 j = 1 229 DO i=1,DIM 230 IF(CoordMax(i)-CoordMin(i) < CoordMax(j)-CoordMin(j)) THEN 231 j = i 232 END IF 233 END DO 234 NormalDirection = j 235 236 ! Direction of maximum change 237 j = 1 238 DO i=1,DIM 239 IF(CoordMax(i)-CoordMin(i) > CoordMax(j)-CoordMin(j)) THEN 240 j = i 241 END IF 242 END DO 243 TangentDirection = j 244 245 ! Find triple point: 246 !------------------- 247 Trip_node = CoordMaxi(TangentDirection) 248 Axis_node = CoordMini(TangentDirection) 249 250 n = Solver % Mesh % MaxElementNodes 251 252 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n), & 253 PNodes % x(n), PNodes % y(n), PNodes % z(n), & 254 x(n), y(n), z(n), Basis(n), dBasisdx(n,3), NodalTemp(n), & 255 Conductivity(n), LatentHeat(n), TempDiff(n), & 256 LocalStiffMatrix(n,n), LocalForceVector(n), LocalMassMatrix(n,n), & 257 PrevTemp(NobNodes), & 258 STAT=istat) 259 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error 1.' ) 260 261 Nodes % x = 0.0d0 262 Nodes % y = 0.0d0 263 Nodes % z = 0.0d0 264 PrevTemp = 0.0d0 265 266 IF( SteadyAlgo ) THEN 267 ALLOCATE( NodeDone( NobNodes ), STAT=istat) 268 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error 3.' ) 269 END IF 270 271 VariableName = ListGetString( Solver % Values, 'Normal Variable', Stat ) 272 IF(Stat) THEN 273 HelpSol => VariableGet( Solver % Mesh % Variables, TRIM(VariableName), ThisOnly=.TRUE. ) 274 ELSE 275 HelpSol => VariableGet( Solver % Mesh % Variables, 'Normals',ThisOnly=.TRUE. ) 276 END IF 277 IF(ASSOCIATED(HelpSol)) THEN 278 Normals => HelpSol % Values 279 NormalsPerm => HelpSol % Perm 280 AverageNormal = .TRUE. 281 ELSE 282 AverageNormal = .FALSE. 283 END IF 284 285 HelpSol => VariableGet( Solver % Mesh % Variables, & 286 TRIM(ComponentName(Solver % Variable))//'Move', stat) 287 IF(ASSOCIATED(HelpSol)) THEN 288 SurfaceMove => HelpSol % Values 289 ELSE 290 ALLOCATE( SurfaceMove(SIZE(Surface)), STAT=istat) 291 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' ) 292 SurfaceMove = 0.0d0 293 CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, & 294 PSolver,TRIM(ComponentName(Solver % Variable))//'Move',1, & 295 SurfaceMove,SurfPerm,Output=.FALSE.) 296 END IF 297 298 IF(ListGetLogical(Solver % Values,'Use Average Velocity',Stat)) THEN 299 HelpSol => VariableGet( Solver % Mesh % Variables, & 300 TRIM(ComponentName(Solver % Variable))//'MoveAve', stat) 301 IF(ASSOCIATED(HelpSol)) THEN 302 SurfaceMoveAve => HelpSol % Values 303 ELSE 304 ALLOCATE( SurfaceMoveAve(SIZE(Surface)), STAT=istat) 305 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' ) 306 SurfaceMoveAve = 0.0d0 307 CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, & 308 PSolver,TRIM(ComponentName(Solver % Variable))//'MoveAve', & 309 1,SurfaceMoveAve,SurfPerm,Output=.FALSE.) 310 END IF 311 END IF 312 313 AllocationsDone = .TRUE. 314 END IF 315 316 ! Find triple point temperature 317 !------------------------------- 318 319 Trip_Temp = Temperature( TempPerm(Trip_node) ) 320 321 322 i = ListGetInteger( Solver % Values,'Passive Steps',Stat) 323 IF( i >= SubroutineVisited) GOTO 200 324 325 ! Find the constant material parameters 326 !-------------------------------------- 327 328 CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (1)) 329 n = CurrentElement % TYPE % NumberOfNodes 330 NodeIndexes => CurrentElement % NodeIndexes 331 k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material') 332 Material => Model % Materials(k) % Values 333 MeltPoint = ListGetConstReal( Material,'Melting Point' ) 334 Density = ListGetConstReal( Material, 'Density' ) 335 336 ! The first velocity should always be set 337 IF ( .NOT. PullVelocitySet ) THEN 338 Upull(1) = ListGetConstReal( Material, 'Convection Velocity 1', Stat ) 339 Upull(2) = ListGetConstReal( Material, 'Convection Velocity 2', Stat ) 340 Upull(3) = 0.0 341 PullVelocitySet = .TRUE. 342 END IF 343 344!-------------------------------------------------------------------- 345! The transient algorithm 346! In the transient algorhitm a heat flux over the interface is computed and 347! it is assumed to be used solely in the melting of the solid into liquid. 348! This melting speed gives an estimate for the melting speed that may be 349! improved by iteration. 350!-------------------------------------------------------------------- 351 352 IF ( TransientAlgo ) THEN 353 StiffMatrix => Solver % Matrix 354 ForceVector => Solver % Matrix % RHS 355 356 UseLoads = ListGetLogical( Solver % Values,'Use Heat Load',Stat) 357 IF(UseLoads) THEN 358 HelpSol => VariableGet( Solver % Mesh % Variables, 'Nodal Heat Load' ) 359 END IF 360 361 PrevUpull = Upull 362 PrevTemp = SurfaceMove 363 364 ! nonlinear iteration could only be associated if normal is computed from the solution 365 DO iter = 1, NonlinearIter 366 367 ! First solve the velocity field 368 CALL InitializeToZero( StiffMatrix, ForceVector ) 369 370 IF(UseLoads) THEN 371 DO t=1,Solver % Mesh % NumberOfNodes 372 i = SurfPerm(t) 373 IF( i <= 0) CYCLE 374 ForceVector(i) = HelpSol % Values(HelpSol % Perm(t)) 375 END DO 376 END IF 377 378 DO t = 1, Solver % NumberOfActiveElements 379 CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (t)) 380 n = CurrentElement % TYPE % NumberOfNodes 381 NodeIndexes => CurrentElement % NodeIndexes 382 383 k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material') 384 Material => Model % Materials(k) % Values 385 386 LatentHeat(1:n) = ListGetReal( Material, 'Latent Heat', n, NodeIndexes ) 387 Conductivity(1:n) = ListGetReal( Material,'Heat Conductivity', n, NodeIndexes ) 388 TempDiff(1:n) = Temperature( TempPerm(NodeIndexes) ) - MeltPoint 389 390 CALL VelocityLocalMatrix( LocalStiffMatrix, LocalMassMatrix, LocalForceVector,& 391 CurrentElement, n, Nodes ) 392 393 CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, & 394 ForceVector, LocalForceVector, n, 1, SurfPerm(NodeIndexes) ) 395 END DO 396 397 CALL FinishAssembly( Solver, ForceVector ) 398 399 ! No Dirihtlet conditions here since 400 ! One should not really try to force the phase change at some point, 401 ! rather use feedback to tune the pull velocity 402 403 CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, SurfaceMove, Norm, 1, Solver ) 404 RelativeChange = Solver % Variable % NonlinChange 405 406 WRITE( Message, * ) 'Result Norm : ',Norm 407 CALL Info( 'PhaseChangeSolve', Message, Level=4 ) 408 409 WRITE( Message, * ) 'Relative Change : ',RelativeChange 410 CALL Info( 'PhaseChangeSolve', Message, Level=4 ) 411 412 IF ( RelativeChange < NonLinearTol ) EXIT 413 END DO 414 415 IF(ListGetLogical(Solver % Values,'Use Average Velocity',Stat)) THEN 416 IF(SurfaceVelocitySet) THEN 417 LocalRelax = GetCReal(Solver % Values,'Velocity Averaging Factor') 418 SurfaceMoveAve = LocalRelax * SurfaceMove + (1-LocalRelax) * SurfaceMoveAve 419 ELSE 420 SurfaceMoveAve = SurfaceMove 421 SurfaceVelocitySet = .TRUE. 422 END IF 423 LocalRelax = GetCReal(Solver % Values,'Velocity Relaxation Factor',Stat) 424 IF(.NOT. Stat) LocalRelax = 1.0d0 425 SurfaceMove = LocalRelax * SurfaceMove + (1-LocalRelax) * SurfaceMoveAve 426 ELSE 427 LocalRelax = GetCReal(Solver % Values,'Velocity Relaxation Factor',Stat) 428 IF(.NOT. Stat) LocalRelax = 1.0d0 429 SurfaceMove = LocalRelax * SurfaceMove + (1-LocalRelax) * PrevTemp 430 END IF 431 432 IF(PullControl .OR. TriplePointFixed) THEN 433 Upull(NormalDirection) = -SurfaceMove(SurfPerm(Trip_node)) 434 IF(PullControl) THEN 435 WRITE(Message,*) 'Pull velocity: ', Upull(NormalDirection) 436 CALL Info('PhaseChangeSolve',Message) 437 END IF 438 END IF 439 440 441 ! Then solve the corresponding update in displacement field 442 SpeedUp = ListGetConstReal( solver % Values,'Transient SpeedUp',Stat) 443 IF(.NOT. Stat) SpeedUp = 1.0d0 444 !----------------------------------------------------------------------------------------- 445 ! If no smoothing is applied directly to the free surface then its update may be computed 446 ! directly from the velocity field. 447 !----------------------------------------------------------------------------------------- 448 IF( ListGetConstReal(Solver % Values,'Surface Smoothing Factor') < 1.0d-20) THEN 449 Surface = Surface + SpeedUp * dt * ( SurfaceMove + Upull(NormalDirection)) 450 ELSE 451 CALL InitializeToZero( StiffMatrix, ForceVector ) 452 453 DO t = 1, Solver % NumberOfActiveElements 454 CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (t)) 455 n = CurrentElement % TYPE % NumberOfNodes 456 NodeIndexes => CurrentElement % NodeIndexes 457 458 CALL SurfaceLocalMatrix( LocalStiffMatrix, LocalMassMatrix, LocalForceVector,& 459 CurrentElement, n, Nodes ) 460 461 CALL Add1stOrderTime( LocalMassMatrix, LocalStiffMatrix, & 462 LocalForceVector, dt, n, 1, SurfPerm(NodeIndexes), Solver ) 463 464 CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, & 465 ForceVector, LocalForceVector, n, 1, SurfPerm(NodeIndexes) ) 466 END DO 467 468 CALL FinishAssembly( Solver, ForceVector ) 469 CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, Surface, Norm, 1, Solver ) 470 END IF 471 472 IF(PullControl .OR. TriplePointFixed) THEN 473 IF(Visited /= Solver % DoneTime) THEN 474 IF(Solver % DoneTime == 1) THEN 475 prevpos0 = 0.0 476 ELSE 477 prevpos0 = pos0 478 END IF 479 Visited = Solver % DoneTime 480 END IF 481 482 pos0 = prevpos0 + UPull(NormalDirection) * dt 483 484 IF(PullControl) THEN 485 ! This sets the maximum position of the crystal side that is still 486 ! aligned with the pull direction. 487 CALL FindPullBoundary() 488 END IF 489 490 END IF 491 END IF 492 493!-------------------------------------------------------------------- 494! Transient algorithm end 495!-------------------------------------------------------------------- 496 497 498!---------------------------------------------------------------------------- 499! The Steady State the simulation is based on a geometric determination of the 500! isotherm. The solution may be accelerated using local or global Newton 501! type of iteration. It is activated only after the solution is quite accurate. 502!----------------------------------------------------------------------------- 503 504 IF( SteadyAlgo ) THEN 505 506 UseTaverage = ListGetLogical( Solver % Values,'Average Temperature', stat ) 507 AveragingDt = ListGetConstReal(Solver % Values,'Steady Averaging Timestep',Stat) 508 509 IF(UseTaverage) THEN 510 HelpSol => VariableGet( Solver % Mesh % Variables, 'Taverage' ) 511 IF(.NOT. ASSOCIATED (HelpSol)) THEN 512 ALLOCATE( Taverage( Model % NumberOfNodes ), STAT=istat ) 513 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' ) 514 CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, & 515 PSolver, 'Taverage', 1, Taverage, TempPerm) 516 HelpSol => VariableGet( Solver % Mesh % Variables, 'Taverage' ) 517 END IF 518 519 NULLIFY(HelpSol2) 520 AverageRelax2 = GetCReal(Solver % Values,'Temperature Averaging Factor',stat) 521 IF(stat) THEN 522 HelpSol2 => VariableGet( Solver % Mesh % Variables, 'Taverage Slow' ) 523 IF(.NOT. ASSOCIATED (HelpSol2)) THEN 524 ALLOCATE( Taverage2( Model % NumberOfNodes ), STAT=istat ) 525 IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' ) 526 CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, & 527 PSolver, 'Taverage Slow', 1, Taverage2, TempPerm) 528 HelpSol2 => VariableGet( Solver % Mesh % Variables, 'Taverage Slow' ) 529 END IF 530 END IF 531 532 IF(dt >= AveragingDt) THEN 533 HelpSol % Values = TempSol % Values 534 IF(ASSOCIATED(HelpSol2)) HelpSol2 % Values = TempSol % Values 535 ELSE 536 AverageRelax = GetCReal( Solver % Values,'Temperature Relaxation Factor') 537 IF(.NOT. ASSOCIATED(HelpSol2)) THEN 538 HelpSol % Values = (1.0d0-AverageRelax) * HelpSol % Values + AverageRelax * TempSol % Values 539 ELSE 540 HelpSol2 % Values = (1.0d0-AverageRelax2) * HelpSol2 % Values + AverageRelax2 * TempSol % Values 541 HelpSol % Values = (1.0d0-AverageRelax) * HelpSol2 % Values + AverageRelax * TempSol % Values 542 END IF 543 PRINT *,'temperature set to average temperature' 544 Temperature => HelpSol % Values 545 END IF 546 END IF 547 548 NewtonAfterIter = ListGetInteger( Solver % Values, & 549 'Nonlinear System Newton After Iterations', stat ) 550 IF ( stat .AND. SubroutineVisited > NewtonAfterIter ) Newton = .TRUE. 551 552 NewtonAfterTol = ListGetConstReal( Solver % Values, & 553 'Nonlinear System Newton After Tolerance', stat ) 554 555 IF(Newton) THEN 556 CALL Info( 'PhaseChangeSolve','Steady state newton formulation', Level=4 ) 557 ELSE 558 CALL Info( 'PhaseChangeSolve','Steady state isoterm formulation', Level=4 ) 559 END IF 560 561! Find the melting point: 562!----------------------- 563 564 IF( TriplePointFixed ) THEN 565 MeltPoint = Temperature( TempPerm(Trip_node) ) 566 WRITE(Message,*) 'Melting point set to triple point temperature',MeltPoint,Trip_node 567 CALL Info('PhaseChangeSolve',Message) 568 END IF 569 570 571! Create the isotherm for T=T_0: 572!------------------------------ 573 574 IF(.NOT. Newton) THEN 575 576 IsoSurfAllocated = .FALSE. 577 xmin = HUGE(xmin) 578 xmax = -HUGE(xmax) 579 580100 NElems = 0 581 582 DO t=1,Solver % Mesh % NumberOfBulkElements 583 584 CurrentElement => Solver % Mesh % Elements(t) 585 586 n = CurrentElement % TYPE % NumberOfNodes 587 NodeIndexes => CurrentElement % NodeIndexes 588 ElementCode = CurrentElement % TYPE % ElementCode 589 590 k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material') 591 IF (.NOT. (ListGetLogical(Model % Materials(k) % Values, 'Solid', stat) .OR. & 592 ListGetLogical(Model % Materials(k) % Values, 'Liquid', stat) )) CYCLE 593 594 IF(ElementCode < 300) CYCLE 595 596 IF ( ANY( TempPerm( NodeIndexes(1:n) ) <= 0 ) ) CYCLE 597 TempDiff(1:n) = Temperature(TempPerm(NodeIndexes(1:n))) - MeltPoint 598 599 IF( ALL ( TempDiff(1:n) < 0.0 ) ) CYCLE 600 IF( ALL ( TempDiff(1:n) > 0.0 ) ) CYCLE 601 602 Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 603 Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 604 605 Vertex = ElementCode/100 606 n=0 607 DO nn=1,Vertex 608 next = MODULO(nn,Vertex) + 1 609 610 temp1 = TempDiff(nn) 611 temp2 = TempDiff(next) 612 613 IF ( ( (temp1 < 0.0) .AND. (0.0 <= temp2) ) .OR. & 614 ( (temp2 <= 0.0) .AND. (0.0 < temp1) ) ) THEN 615 616 n = n + 1 617 618 IF ( n <= 2 ) THEN 619 NElems = NElems + 1 620 IF(IsoSurfAllocated) THEN 621 IsoSurf(NElems,1) = Nodes % x(nn) + & 622 temp1 * ((Nodes % x(next) - Nodes % x(nn)) / (temp1-temp2)) 623 IsoSurf(NElems,2) = Nodes % y(nn) + & 624 temp1 * ((Nodes % y(next) - Nodes % y(nn)) / (temp1-temp2)) 625 626 xmin = MIN( IsoSurf(Nelems,1), xmin ) 627 xmax = MAX( IsoSurf(Nelems,1), xmax ) 628 END IF 629 ELSE 630 CALL Warn('PhaseChangeSolve','Wiggly Isotherm') 631 WRITE(Message,*) 'nodeindexes',nodeindexes(1:Vertex) 632 CALL Warn('PhaseChangeSolve',Message) 633 WRITE(Message,*) 'temperature',Temperature(TempPerm(NodeIndexes)) 634 CALL Warn('PhaseChangeSolve',Message) 635 PRINT *,'Nodes % x',Nodes % x(n) 636 PRINT *,'Nodes % y',Nodes % y(n) 637 END IF 638 END IF 639 END DO 640 641 IF ( n == 1 ) THEN 642 NElems = NElems - 1 643 CYCLE 644 END IF 645 646 IF (IsoSurfAllocated .AND. n == 2) THEN 647 IF ( IsoSurf(Nelems-1,TangentDirection) > IsoSurf(Nelems,TangentDirection) ) THEN 648 Temppi = IsoSurf(Nelems-1,1) 649 IsoSurf(Nelems-1,1) = IsoSurf(Nelems,1) 650 IsoSurf(Nelems,1) = Temppi 651 Temppi = IsoSurf(Nelems-1,2) 652 IsoSurf(Nelems-1,2) = IsoSurf(Nelems,2) 653 IsoSurf(Nelems,2) = Temppi 654 END IF 655 END IF 656 657 END DO 658 659 IF(Nelems == 0) CALL Fatal('PhaseChangeSolve','Isotherm is empty thus cannot map phase change surface') 660 661 IF(.NOT. IsoSurfAllocated) THEN 662 ALLOCATE( IsoSurf(Nelems+1,2)) 663 IsoSurfAllocated = .TRUE. 664 WRITE(Message,*) 'Isotherm created with number of segments',Nelems 665 CALL Info('PhaseChangeSolve',Message) 666 GOTO 100 667 END IF 668 END IF 669 670 area = 0.0 671 volume = 0.0 672 tave = 0.0 673 tabs = 0.0 674 volabs = 0.0 675 NodeDone = .FALSE. 676 MaxTempDiff = 0.0 677 678 679 DO t = 1, Solver % NumberOfActiveElements 680 681 CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements(t)) 682 n = CurrentElement % TYPE % NumberOfNodes 683 NodeIndexes => CurrentElement % NodeIndexes 684 685 ElementCode = CurrentElement % TYPE % ElementCode 686 IF(ElementCode < 200) CYCLE 687 IF(ElementCode > 203) THEN 688 CALL Fatal('PhaseChangeSolve','Implemented only for elements 202 and 203!') 689 CYCLE 690 END IF 691 692 Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 693 Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 694 Nodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes) 695 696 TempDiff(1:n) = Temperature( TempPerm(NodeIndexes(1:n)) ) - MeltPoint 697 MaxTempDiff = MAX(MaxTempDiff, MAXVAL(ABS(TempDiff(1:n)))) 698 699 DO nn=1,n 700 701 k = SurfPerm(NodeIndexes(nn)) 702 IF ( NodeDone(k) ) CYCLE 703 NodeDone(k) = .TRUE. 704 705 IF( TriplePointFixed .AND. NodeIndexes(nn) == trip_node) THEN 706 SurfaceMove(k) = 0.0d0 707 PRINT *,'triple point fixed by construction' 708 CYCLE 709 END IF 710 711 IF ( nn == 3 ) THEN 712 Surface(k)=0.5*SUM(Surface(SurfPerm(NodeIndexes(1:2)))) 713 Update = 0.0 714 END IF 715 716 IF(TangentDirection == 1) THEN 717 xx = Nodes % x(nn) 718 yy = Nodes % y(nn) 719 ELSE 720 xx = Nodes % y(nn) 721 yy = Nodes % x(nn) 722 END IF 723 724 725 IF ( .NOT. Newton ) THEN 726 727 ! Find the the contour element that has the x-coordinate in closest to the that of the 728 ! free surface 729 730 Eps = 1.0d-6 * ( xmax - xmin ) 731 732 dxmin = HUGE(dxmin) 733 dymin = HUGE(dymin) 734 stat = .FALSE. 735 736 DO i=1,Nelems-1,2 737 738 IF ( (xx > IsoSurf(i,TangentDirection) - Eps) .AND. (xx < IsoSurf(i+1,TangentDirection) + Eps)) THEN 739 dxmin = 0.0 740 d = MIN( ABS(yy - IsoSurf(i,NormalDirection)), ABS(yy - IsoSurf(i+1,NormalDirection)) ) 741 742 ! Punish for overlapping the boundaries 743 d = d + MAX(0.0d0, IsoSurf(i,TangentDirection) - xx) 744 d = d + MAX(0.0d0, xx - IsoSurf(i+1,TangentDirection) ) 745 746 IF(d <= dymin) THEN 747 stat = .TRUE. 748 dymin = d 749 imin = i 750 END IF 751 END IF 752 753 IF(.NOT. stat) THEN 754 d = MIN( ABS(xx - IsoSurf(i,TangentDirection)), ABS(xx - IsoSurf(i+1,TangentDirection)) ) 755 IF (d <= dxmin) THEN 756 dxmin = d 757 imin = i 758 END IF 759 END IF 760 END DO 761 762 i = imin 763 764 ! There may be a problem if the boundary cannot be mapped on an isotherm 765 IF (.NOT. stat) THEN 766 IF(dxmin > 1.0d-2* ABS(IsoSurf(i,TangentDirection)- IsoSurf(i+1,TangentDirection))) THEN 767 CALL Warn('PhaseChangeSolve','Isotherm error?') 768 WRITE(Message,*) 'Nodeindexes',NodeIndexes(nn) 769 CALL Warn('PhaseChangeSolve',Message) 770 WRITE(Message,*) 'x:',xx,' y:',yy 771 CALL Warn('PhaseChangeSolve',Message) 772 WRITE(Message,*) 'dxmin:',dxmin,' dymin:',dymin 773 CALL Warn('PhaseChangeSolve',Message) 774 END IF 775 END IF 776 777 IF ( ABS( IsoSurf(i+1,TangentDirection) - IsoSurf(i,TangentDirection) ) > AEPS ) THEN 778 Update = IsoSurf(i,NormalDirection) + ( xx - IsoSurf(i,TangentDirection) ) * & 779 ( IsoSurf(i+1,NormalDirection) - IsoSurf(i,NormalDirection) ) / & 780 ( IsoSurf(i+1,TangentDirection) - IsoSurf(i,TangentDirection) ) - yy 781 ELSE 782 Update = 0.5d0 * ( IsoSurf(i,NormalDirection) + IsoSurf(i+1,NormalDirection) ) - yy 783 END IF 784 785 END IF 786 787 TTemp = Temperature( TempPerm(NodeIndexes(nn)) ) 788 789 IF ( Newton ) THEN 790 dTdz = TTemp - PrevTemp(k) 791 IF ( ABS(dTdz) < AEPS ) THEN 792 CALL Warn( 'PhaseChangeSolve', 'Very small temperature update.' ) 793 dTdz = 1 794 END IF 795 Update = SurfaceMove(k) * ( MeltPoint - TTemp ) / dTdz 796 END IF 797 798 ! This enforcing is rather than by setting meltpoint to triple point temperature 799 ! IF ( NodeIndexes(nn) == Trip_node ) Update = 0 800 801 PrevTemp(k) = TTemp 802 SurfaceMove(k) = Update 803 END DO 804 805 IntegStuff = GaussPoints( CurrentElement ) 806 DO i=1,IntegStuff % n 807 808 u = IntegStuff % u(i) 809 v = IntegStuff % v(i) 810 w = IntegStuff % w(i) 811 812 stat = ElementInfo( CurrentElement, Nodes, u, v, w, detJ, Basis, dBasisdx ) 813 814 s = IntegStuff % s(i) * detJ 815 816 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 817 s = s * SUM(Basis(1:n) * Nodes % x(1:n)) * 2.0 * PI 818 END IF 819 820 area = area + S 821 volume = volume + S * SUM(Basis(1:n) * SurfaceMove(SurfPerm(NodeIndexes(1:n)))) 822 tave = tave + S * SUM(Basis(1:n) * TempDiff(1:n) ) 823 volabs = volabs + S * SUM(Basis(1:n) * ABS(SurfaceMove(SurfPerm(NodeIndexes(1:n))))) 824 tabs = tabs + S * SUM(Basis(1:n) * ABS(TempDiff(1:n)) ) 825 END DO 826 827 END DO 828 829 830! IF(.NOT. UseTAverage .AND. dt < SteadyDt) THEN 831! LocaRelax = Relax * GetCReal( Solver % Values,'Temperature Relaxation Factor') 832! ELSE 833 LocalRelax = Relax 834! END IF 835 836 ! There are several different acceleration methods which are mainly inactive 837 tave = tave / area 838 tabs = tabs / area 839 volume = volume / area 840 volabs = volabs / area 841 842 i = ListGetInteger(Solver % Values,'Lumped Newton After Iterations', Stat) 843 IF(Stat .AND. SubroutineVisited > i) THEN 844 845 j = ListGetInteger(Solver % Values,'Lumped Newton Mode', Stat) 846 SELECT CASE( j ) 847 CASE( 1 ) 848 cvol = 0.5*(prevtave+tave)/(prevtave-tave) 849 850 CASE( 2 ) 851 cvol = 0.5*(prevvolabs+volabs)/(prevvolabs-volabs) 852 853 CASE( 3 ) 854 cvol = 0.5*(prevtabs+tabs)/(prevtabs-tabs) 855 856 CASE DEFAULT 857 cvol = 0.5*(prevvolume+volume)/(prevvolume-volume) 858 859 END SELECT 860 861 IF(cvol < 0.0) THEN 862 cvol = 1.0 863 ccum = 1.0 864 END IF 865 866 clim = ListGetConstReal(Solver % Values,'Lumped Newton Limiter', Stat) 867 IF(.NOT. Stat) clim = 100.0 868 cvol = MIN(clim,cvol) 869 cvol = MAX(1.0/clim,cvol) 870 871 ccum = ccum * cvol 872 873 WRITE(Message,*) 'Lumped Newton relaxation: ', ccum 874 CALL Info('PhaseChangeSolve',Message) 875 876 LocalRelax = LocalRelax * ccum 877 END IF 878 879 SurfaceMove = LocalRelax * SurfaceMove 880 Surface = Surface + SurfaceMove 881 882 dpos = SurfaceMove(SurfPerm(Trip_node)) 883 884 MaxUpdate = MAXVAL(ABS(SurfaceMove)) / MAXVAL(ABS(Surface)) 885 886 IF ( ABS(MaxUpdate) < NewtonAfterTol ) Newton = .TRUE. 887 WRITE(Message,*) 'Maximum surface update: ', MaxUpdate 888 CALL Info('PhaseChangeSolve',Message) 889 890 WRITE(Message,*) 'Maximum temperature difference: ', MaxTempDiff 891 CALL Info('PhaseChangeSolve',Message) 892 893 Norm = SQRT( SUM( Surface**2 ) / SIZE( Surface ) ) 894 WRITE( Message, * ) 'Result Norm : ',Norm 895 CALL Info( 'PhaseChangeSolve', Message, Level=4 ) 896 Solver % Variable % Norm = Norm 897 898 prevvolume = volume 899 prevtave = tave 900 prevtabs = tabs 901 prevvolabs = volabs 902 903 IF(IsoSurfAllocated) THEN 904 IsoSurfAllocated = .FALSE. 905 DEALLOCATE(IsoSurf) 906 END IF 907 END IF 908 909!-------------------------------------------------------------------- 910! Steady algorithm end 911!-------------------------------------------------------------------- 912 913 914200 CALL ListAddConstReal(Model % Simulation,'res: Triple point temperature',Trip_temp) 915 CALL ListAddConstReal( Model % Simulation,'res: triple point movement',dpos) 916 CALL ListAddConstReal(Model % Simulation,'res: Pull Position',pos0) 917 918 IF(PullControl) THEN 919 IF(NormalDirection == 1) THEN 920 CALL ListAddConstReal( Model % Simulation,'res: Pull Velocity 1',UPull(NormalDirection)) 921 ELSE 922 CALL ListAddConstReal( Model % Simulation,'res: Pull Velocity 2',UPull(NormalDirection)) 923 END IF 924 END IF 925 926 FirstTime = .FALSE. 927 928!------------------------------------------------------------------------------ 929CONTAINS 930 931!------------------------------------------------------------------------------ 932 SUBROUTINE VelocityLocalMatrix( StiffMatrix, MassMatrix, ForceVector,& 933 Element, nCoord, Nodes ) 934 935 ! external variables: 936 REAL(KIND=dp) :: StiffMatrix(:,:), MassMatrix(:,:), ForceVector(:) 937 INTEGER :: nCoord 938 TYPE(Nodes_t) :: Nodes 939 TYPE(Element_t), POINTER :: Element 940 941 ! internal variables: 942 TYPE(Nodes_t) :: PNodes 943 TYPE(Element_t), POINTER :: Parent 944 REAL(KIND=dp) :: Basis(3*nCoord),dBasisdx(3*nCoord,3), & 945 X,Y,Z,U,V,W,S,detJ, TGrad(3,3),Flux,pu,pv,pw,pull, LocalHeat, & 946 NodalTemp(3*nCoord), xx(10),yy(10),zz(10), NodalSurf(nCoord), NodalNormal(2,nCoord) 947 REAL(KIND=dp) :: Velo, Ny, Normal(3), StabFactor, StabCoeff, DerSurf, xcoord 948 LOGICAL :: Stat 949 INTEGER :: i,j,k,l,t,p,q, n, pn 950 TYPE(GaussIntegrationPoints_t) :: IntegStuff 951 952!------------------------------------------------------------------------------ 953 954 ForceVector = 0.0d0 955 StiffMatrix = 0.0d0 956 MassMatrix = 0.0d0 957 n = nCoord 958 959 StabFactor = ListGetConstReal(Solver % Values,'Velocity Smoothing Factor') 960 961 Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 962 Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 963 NodalSurf(1:n) = Surface( SurfPerm(NodeIndexes) ) 964 965 IF(AverageNormal) THEN 966 NodalNormal(1,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))-1) 967 NodalNormal(2,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))) 968 END IF 969 970 IF(TransientAlgo) THEN 971 ALLOCATE( PNodes % x(10), PNodes % y(10), PNodes % z(10) ) 972 END IF 973 974! Numerical integration: 975! ---------------------- 976 977 IntegStuff = GaussPoints( Element ) 978 979 DO t = 1,IntegStuff % n 980 981 u = IntegStuff % u(t) 982 v = IntegStuff % v(t) 983 w = IntegStuff % w(t) 984 s = IntegStuff % s(t) 985 986! Basis function values & derivatives at the integration point: 987! ------------------------------------------------------------- 988 stat = ElementInfo( Element,Nodes,U,V,W,detJ,Basis,dBasisdx) 989 s = s * detJ 990 xcoord = SUM( Nodes % x(1:nCoord) * Basis(1:nCoord) ) 991 992! IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 993! s = s * xcoord 994! END IF 995 996 IF(AverageNormal) THEN 997 Normal(1) = SUM( Basis(1:n) * NodalNormal(1,1:n)) 998 Normal(2) = SUM( Basis(1:n) * NodalNormal(2,1:n)) 999 Normal(3) = 0.0d0 1000 ELSE 1001 Normal = NormalVector( Element, Nodes, u, v, .TRUE. ) 1002 END IF 1003 1004 LocalHeat = SUM(Basis(1:n) * LatentHeat(1:n)) 1005 StabCoeff = StabFactor * LocalHeat * Density 1006 1007 1008 IF(UseLoads) THEN 1009 ! do nothing, loads already computed 1010 1011 ELSE 1012 ! Compute the flux from normal derivaties 1013 TGrad = 0.0d0 1014 l = 0 1015 DO i=1,2 1016 1017 IF( i == 1) THEN 1018 Parent => Element % BoundaryInfo % Left 1019 ELSE 1020 Parent => Element % BoundaryInfo % Right 1021 END IF 1022 1023 k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material') 1024 IF (ListGetLogical(Model % Materials(k) % Values,'Solid',stat)) THEN 1025 IF(l == 2) CALL Fatal('PhaseChangeSolve','Both materials cannot be solid!') 1026 l = 2 1027 ELSE 1028 IF(l == 1) CALL Fatal('PhaseChangeSolve','Both materials cannot be liquid!') 1029 l = 1 1030 END IF 1031 1032 pn = Parent % TYPE % NumberOfNodes 1033 k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material') 1034 1035 Conductivity(1:pn) = ListGetReal( Model % Materials(k) % Values, & 1036 'Heat Conductivity', pn, Parent % NodeIndexes ) 1037 NodalTemp(1:pn) = Temperature( TempPerm(Parent % NodeIndexes) ) 1038 1039 stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis, dBasisdx) 1040 1041 ! Calculate the basis functions for the parent element: 1042 ! ----------------------------------------------------- 1043 DO j = 1,n 1044 DO k = 1,pn 1045 IF( NodeIndexes(j) == Parent % NodeIndexes(k) ) THEN 1046 xx(j) = Parent % TYPE % NodeU(k) 1047 yy(j) = Parent % TYPE % NodeV(k) 1048 zz(j) = Parent % TYPE % NodeW(k) 1049 EXIT 1050 END IF 1051 END DO 1052 END DO 1053 1054 pu = SUM( Basis(1:n) * xx(1:n) ) 1055 pv = SUM( Basis(1:n) * yy(1:n) ) 1056 pw = SUM( Basis(1:n) * zz(1:n) ) 1057 1058 PNodes % x(1:pn) = Solver % Mesh % Nodes % x(Parent % NodeIndexes) 1059 PNodes % y(1:pn) = Solver % Mesh % Nodes % y(Parent % NodeIndexes) 1060 PNodes % z(1:pn) = Solver % Mesh % Nodes % z(Parent % NodeIndexes) 1061 1062 stat = ElementInfo( Parent, PNodes, pu, pv, pw, detJ, Basis, dBasisdx ) 1063 1064 DO j=1,DIM 1065 TGrad(l,j) = SUM( Conductivity(1:pn) * Basis(1:pn) ) * & 1066 SUM( dBasisdx(1:pn,j) * NodalTemp(1:pn) ) 1067 END DO 1068 END DO 1069 1070 Flux = SUM( (TGrad(1,:) - TGrad(2,:)) * Normal) 1071 stat = ElementInfo( Element,Nodes,U, V, W, detJ, Basis, dBasisdx ) 1072 END IF 1073 1074 ! Assembly the matrix 1075 DO p=1,n 1076 DO q=1,n 1077 StiffMatrix(p,q) = StiffMatrix(p,q) + s * Basis(p) * Basis(q) * & 1078 Normal(NormalDirection) * Density * LocalHeat 1079 IF(NodeIndexes(p) /= Trip_node) THEN 1080 StiffMatrix(p,q) = StiffMatrix(p,q) + & 1081 s * StabCoeff * dBasisdx(q,TangentDirection) * dBasisdx(p,TangentDirection) 1082 END IF 1083 1084 ! BC for tipple node 1085! IF(NodeIndexes(p) == Trip_node) THEN 1086! StiffMatrix(p,q) = StiffMatrix(p,q) - & 1087! StabCoeff * Basis(p) * dBasisdx(q,TangentDirection) 1088! END IF 1089 END DO 1090 1091 ! transient part of heat flux 1092 IF(.NOT. (UseLoads)) THEN 1093 ForceVector(p) = ForceVector(p) - s * Basis(p) * Flux 1094 END IF 1095 1096 END DO 1097 1098 END DO 1099 1100 IF(TransientAlgo) DEALLOCATE( PNodes % x, PNodes % y, PNodes % z ) 1101 1102!------------------------------------------------------------------------------ 1103 END SUBROUTINE VelocityLocalMatrix 1104!------------------------------------------------------------------------------ 1105 1106 1107!------------------------------------------------------------------------------ 1108 SUBROUTINE SurfaceLocalMatrix( StiffMatrix, MassMatrix, ForceVector,& 1109 Element, nCoord, Nodes ) 1110 1111 ! external variables: 1112 REAL(KIND=dp) :: StiffMatrix(:,:), MassMatrix(:,:), ForceVector(:) 1113 INTEGER :: nCoord 1114 TYPE(Nodes_t) :: Nodes 1115 TYPE(Element_t), POINTER :: Element 1116 1117 ! internal variables: 1118 REAL(KIND=dp) :: Basis(3*nCoord),dBasisdx(3*nCoord,3), & 1119 X,Y,Z,U,V,W,S,detJ,pull, NodalVelo(nCoord), xcoord, NodalNormal(2,nCoord) 1120 REAL(KIND=dp) :: Velo, StabFactor, StabCoeff 1121 LOGICAL :: Stat 1122 INTEGER :: i,j,k,l,t,p,q, n, pn 1123 TYPE(GaussIntegrationPoints_t) :: IntegStuff 1124 1125!------------------------------------------------------------------------------ 1126 1127 ForceVector = 0.0d0 1128 StiffMatrix = 0.0d0 1129 MassMatrix = 0.0d0 1130 n = nCoord 1131 1132 StabFactor = ListGetConstReal(Solver % Values,'Surface Smoothing Factor') 1133 StabCoeff = StabFactor / dt 1134 1135 Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 1136 Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 1137 NodalVelo(1:n) = SurfaceMove( SurfPerm(NodeIndexes) ) 1138 1139 IF(AverageNormal) THEN 1140 NodalNormal(1,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))-1) 1141 NodalNormal(2,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))) 1142 END IF 1143 1144 1145! Numerical integration: 1146! ---------------------- 1147 1148 IntegStuff = GaussPoints( Element ) 1149 1150 DO t = 1,IntegStuff % n 1151 1152 u = IntegStuff % u(t) 1153 v = IntegStuff % v(t) 1154 w = IntegStuff % w(t) 1155 s = IntegStuff % s(t) 1156 1157! Basis function values & derivatives at the integration point: 1158! ------------------------------------------------------------- 1159 stat = ElementInfo( Element,Nodes,U,V,W,detJ,Basis,dBasisdx) 1160 s = s * detJ 1161 xcoord = SUM( Nodes % x(1:nCoord) * Basis(1:nCoord) ) 1162 1163 IF(AverageNormal) THEN 1164 Normal(1) = SUM( Basis(1:n) * NodalNormal(1,1:n) ) 1165 Normal(2) = SUM( Basis(1:n) * NodalNormal(2,1:n) ) 1166 Normal(3) = 0.0d0 1167 ELSE 1168 Normal = NormalVector( Element, Nodes, u, v, .TRUE. ) 1169 END IF 1170 1171 Velo = SUM( Basis(1:n) * NodalVelo(1:n)) + Upull(NormalDirection) 1172 Velo = SpeedUp * Velo 1173 1174 DO p=1,n 1175 DO q=1,n 1176 MassMatrix(p,q) = MassMatrix(p,q) + s * Basis(p) * Basis(q) 1177 StiffMatrix(p,q) = StiffMatrix(p,q) + & 1178 s * StabCoeff * dBasisdx(q,TangentDirection) * dBasisdx(p,TangentDirection) 1179 1180 ! BC for tipple node 1181 IF(NodeIndexes(p) == Trip_node) THEN 1182 StiffMatrix(p,q) = StiffMatrix(p,q) - & 1183 StabCoeff * Basis(p) * dBasisdx(q,TangentDirection) 1184 END IF 1185 END DO 1186 ForceVector(p) = ForceVector(p) + s * Basis(p) * Velo 1187 END DO 1188 1189 END DO 1190 1191!------------------------------------------------------------------------------ 1192 END SUBROUTINE SurfaceLocalMatrix 1193!------------------------------------------------------------------------------ 1194 1195 1196!------------------------------------------------------------------------------ 1197 SUBROUTINE FindPullBoundary() 1198 1199 REAL (KIND=dp) :: Ybot, Ymax, Ytop, Xtrip, x, y 1200 INTEGER :: t,k,i,n 1201 LOGICAL :: PullBoundary 1202 1203 IF(TangentDirection /= 1) CALL Fatal('PhaseChangeSolve',& 1204 'FindPullBoundary implemented only for lateral boundaries!') 1205 1206 PullBoundary = .FALSE. 1207 Ybot = Solver % Mesh % Nodes % y(Trip_node) 1208 Ymax = MAXVAL(Solver % Mesh % Nodes % y) 1209 Ytop = Ymax 1210 Xtrip = Solver % Mesh % Nodes % x(Trip_node) 1211 1212 DO t = Solver % Mesh % NumberOfBulkElements + 1, & 1213 Solver % Mesh % NumberOfBulkElements + Solver % Mesh % NumberOfBoundaryElements 1214 1215 CurrentElement => Solver % Mesh % Elements(t) 1216 Model % CurrentElement => CurrentElement 1217 n = CurrentElement % TYPE % NumberOfNodes 1218 NodeIndexes => CurrentElement % NodeIndexes 1219 1220 DO k=1, Model % NumberOfBCs 1221 IF ( Model % BCs(k) % Tag /= CurrentElement % BoundaryInfo % Constraint ) CYCLE 1222 IF( ListGetLogical(Model % BCs(k) % Values,'Pull Boundary',stat ) ) THEN 1223 PullBoundary = .TRUE. 1224 DO i = 1,n 1225 x = Solver % Mesh % Nodes % x(NodeIndexes(i)) 1226 y = Solver % Mesh % Nodes % y(NodeIndexes(i)) 1227 IF(y > Ybot .AND. ABS(x-Xtrip) > 1.0e-6 * (Ymax-Ybot)) Ytop = y 1228 END DO 1229 END IF 1230 END DO 1231 END DO 1232 1233 IF (PullBoundary) THEN 1234 CALL ListAddConstReal(Model % Simulation,'res: Triple point position',Ybot) 1235 CALL ListAddConstReal(Model % Simulation,'res: Full pull position',Ytop) 1236 END IF 1237 1238 END SUBROUTINE FindPullBoundary 1239 1240!------------------------------------------------------------------------------ 1241 END SUBROUTINE PhaseChangeSolve 1242!------------------------------------------------------------------------------ 1243 1244 1245!------------------------------------------------------------------------------- 1246 FUNCTION MeltingHeat(Model, Node, t) RESULT(Flux) 1247!------------------------------------------------------------------------------- 1248! This subroutine computes the heat flux resulting from solidification 1249! This 1250!------------------------------------------------------------------------------- 1251 USE Types 1252 USE Lists 1253 USE ElementDescription 1254 IMPLICIT NONE 1255!------------------------------------------------------------------------------- 1256 TYPE(Model_t) :: Model 1257 INTEGER:: Node 1258 REAL (KIND=dp):: Flux,t 1259!------------------------------------------------------------------------------- 1260 TYPE(Variable_t), POINTER :: NormalSol 1261 INTEGER:: k,n,i 1262 INTEGER, POINTER :: NodeIndexes(:), NormalsPerm(:) 1263 REAL (KIND=dp):: NodeLatentHeat, Density, NormalPull, u, v 1264 REAL (KIND=dp):: UPull(3) = (/ 0,0,0 /), Normal(3), ElemLatentHeat(4) 1265 REAL (KIND=dp), POINTER :: Normals(:) 1266 LOGICAL:: stat, NormalExist = .FALSE., Visited = .FALSE. 1267 TYPE(Nodes_t) :: Nodes 1268 TYPE(Element_t), POINTER :: CurrentElement, Parent 1269 1270!------------------------------------------------------------------------------ 1271 1272 SAVE NormalExist, Normals, NormalsPerm, Nodes 1273 1274 1275 IF(.NOT. Visited) THEN 1276 NormalSol => VariableGet( Model % Variables, 'Normals',ThisOnly=.TRUE. ) 1277 IF(ASSOCIATED(NormalSol)) THEN 1278 NormalsPerm => NormalSol % Perm 1279 Normals => NormalSol % Values 1280 NormalExist = .TRUE. 1281 ELSE 1282 n = Model % Mesh % MaxElementNodes 1283 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 1284 END IF 1285 Visited = .TRUE. 1286 END IF 1287 1288 CurrentElement => Model % CurrentElement 1289 NodeIndexes => CurrentElement % NodeIndexes 1290 n = CurrentElement % TYPE % NumberOfNodes 1291 1292 k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material') 1293 ElemLatentHeat(1:n) = ListGetReal( Model % Materials(k) % Values, 'Latent Heat', n, NodeIndexes ) 1294 1295 DO i=1,n 1296 IF(NodeIndexes(i) == Node) EXIT 1297 END DO 1298 IF(NodeIndexes(i) /= Node) CALL Fatal('PhaseProcs','Node not found') 1299 NodeLatentHeat = ElemLatentHeat(i) 1300 1301 UPull = 0.0 1302 UPull(1) = ListGetConstReal(Model % Simulation,'res: Pull Velocity 1',stat) 1303 IF(.NOT. stat) UPull(1) = ListGetConstReal(Model % Materials(k) % Values,'Convection Velocity 1',stat) 1304 1305 UPull(2) = ListGetConstReal(Model % Simulation,'res: Pull Velocity 2',stat) 1306 IF(.NOT. stat) UPull(2) = ListGetConstReal(Model % Materials(k) % Values,'Convection Velocity 2',stat) 1307 1308 Parent => CurrentElement % BoundaryInfo % Left 1309 k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material') 1310 IF (ListGetLogical(Model % Materials(k) % Values, 'Solid', stat)) THEN 1311 Density = ListGetConstReal( Model % Materials(k) % Values, 'Density' ) 1312 ELSE 1313 Parent => CurrentElement % BoundaryInfo % Right 1314 k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material') 1315 Density = ListGetConstReal( Model % Materials(k) % Values, 'Density' ) 1316 END IF 1317 1318 IF(.NOT. ASSOCIATED(Parent)) CALL Fatal('PhaseProcs','Parent not found') 1319!------------------------------------------------------------------------------ 1320 1321 IF( NormalExist ) THEN 1322 Normal(1) = Normals( 2*NormalsPerm( Node ) - 1 ) 1323 Normal(2) = Normals( 2*NormalsPerm( Node )) 1324 Normal(3) = 0.0d0 1325 ELSE 1326 Nodes % x(1:n) = Model % Nodes % x(NodeIndexes) 1327 Nodes % y(1:n) = Model % Nodes % y(NodeIndexes) 1328 Nodes % z(1:n) = Model % Nodes % z(NodeIndexes) 1329 1330 ! For line segments the normal in the center is usually sufficient 1331 u = 0.0d0 1332 v = 0.0d0 1333 1334 ! If inner boundary, Normal Target Body should be defined for the boundary 1335 ! (if not, material density will be used to determine then normal direction 1336 ! and should be defined for bodies on both sides): 1337 1338 Normal = NormalVector( CurrentElement, Nodes, u, v, .TRUE. ) 1339 END IF 1340 1341 NormalPull = SUM( Normal * UPull ) 1342 1343 Flux = NodeLatentHeat * Density * NormalPull 1344 1345!------------------------------------------------------------------------------ 1346END FUNCTION MeltingHeat 1347!------------------------------------------------------------------------------ 1348