Smart Heater Control

Numerical methods and mathematical models of Elmer
Post Reply
MrChips
Posts: 33
Joined: 12 Jul 2016, 00:16
Antispam: Yes
Location: Canada

Smart Heater Control

Post by MrChips »

I am using savescalars to save 'heater power density' of a smart heater setup. When I run the smart heater I put in a heat source value for the heater. If I don't do this it will say it is not defined when I run the solver. I assume the heat source value I put in is just used as the initial guess and the solver iterates until it reaches the temperature setpoint of the heater. However, the results it comes up with for the power density for a given temperature set-point seem to vary quite a bit depending on the initial guess for the heat source. I don't know if I am using this correctly in this transient case or there is something else going on.

If the initial guess of the heat source is 500 W/kg with a setpoint temperature of the smart heater of 923 K I get the following results:

Time (s)--------------------Power Density (W/kg)
3.600000000000E+003 9.551616693502E+002
7.200000000000E+003 5.279649006690E+002
1.080000000000E+004 4.409886676690E+002
1.440000000000E+004 3.996450815166E+002
1.800000000000E+004 3.776308097251E+002
2.160000000000E+004 3.641013732455E+002
2.520000000000E+004 3.547031022428E+002

However, if the initial guess of the heat source is 90 W/kg with a setpoint temperature of the smart heater of 923 K, I get these results:
Time (s)---------------------Power Density (W/kg)
3.600000000000E+003 5.492029620666E+002
7.200000000000E+003 1.953290766480E+002
1.080000000000E+004 1.631571755695E+002
1.440000000000E+004 1.478636815313E+002
1.800000000000E+004 1.397188460322E+002
2.160000000000E+004 1.347129985598E+002
2.520000000000E+004 1.312368548457E+002

Edit: I should note that the heater temperature in both cases is the same when I check it in post.

Here's my sif file:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Transient
  Steady State Max Iterations = 30
  Timestepping Method = BDF
  BDF Order = 1   
  ! Timestep Intervals(3) = 1 48 179 
  ! Timestep Sizes(3) = 0.001 1800 86400
  Timestep Intervals(1) = 180 
  Timestep Sizes(1) = 3600
  ! Output Intervals(3) = 1 24 1
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

!---------------------------------------------------------------------

Initial Condition 1
  Name = "InitialCondition 1"
  Temperature = 293.15
End

!---------------------------------------------------------------------

Body Force 1
  Name = "Heater"
  Heat Source = 93.86
  Smart Heater Control = Logical True
  Smart Heater Temperature = Real MATC "273.15 + 650"
  Smart Heater Control Point(3) = 0.0281585 0.00813963 0.0
End

!---------------------------------------------------------------------

Material 1
  Name = "Soil"
  Heat Conductivity = Variable Temperature
   Real
    Include TC.dat
   End
  Heat Capacity = 800
  Density = 2200
End

Material 2
  Name = "Steel"
  Sound speed = 5100.0
  Youngs modulus = 200.0e9
  Poisson ratio = 0.285
  Mesh Poisson ratio = 0.285
  Heat expansion Coefficient = 13.8e-6
  Heat Conductivity = 44.8
  Density = 7850.0
  Heat Capacity = 1265.0
End

Material 3
  Name = "Air"
  Heat Capacity = Variable Temperature
   Real
    Include CP-Air.dat
   End
  Heat Conductivity = Variable Temperature
   Real
    Include TC-Air.dat
   End
  Density = Variable Temperature
   Real
    Include D-Air.dat
   End
End

!---------------------------------------------------------------------

Body 1
  Name = "Core"
  Equation = 1
  Material = 2
  Initial condition = 1
  Body Force = 1
End

Body 2
  Name = "Casing"
  Equation = 1
  Material = 2
  Initial condition = 1
End

Body 3
  Name = "Soil"
  Equation = 1
  Material = 1
  Initial condition = 1
End

Body 4
  Name = "AirInner"
  Equation = 1
  Material = 3
  Initial condition = 1
End

Body 5
  Name = "AirOuter"
  Equation = 1
  Material = 3
  Initial condition = 1
End

!---------------------------------------------------------------------

Boundary Condition 1
  Name = "Heater"
  Target Boundaries(1) = 1  
  Save Scalars = True
  Radiation Boundary = 1
  Radiation = Diffuse Gray
  Emissivity = 0.80
End

Boundary Condition 2
  Name = "InnerCasing"
  Target Boundaries(1) = 4
  Radiation Boundary = 1
  Radiation = Diffuse Gray
  Emissivity = 0.80
End

!---------------------------------------------------------------------

Equation 1
  Name = "Equation 1"
  Active Solvers(3) = 1 2 3
End

!---------------------------------------------------------------------

Solver 1
  Equation = Heat Equation
  Procedure = "HeatSolve" "HeatSolver"
  Variable = Temperature
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Calculate Loads = True
  Calculate Boundary Weights = True
  Calculate Flux = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-10
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = Diagonal
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Precondition Recompute = 1
  ! Update View Factors = True
  ! Update Gebhardt Factors = True
End

Solver 2
  Exec Solver = Before Saving
  Equation = SaveScalars
  Procedure = "SaveData" "SaveScalars"
  Filename = save scalars.dat
  Variable 1 = Time
  ! Variable 2 = Temperature
  ! Operator 2 = max
  ! Target Variable 2 = String maxTemp
  ! Variable 3 = Heat Source  
End

Solver 3
  Equation = Result Output
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output Format = Vtu
  vtu time collection = logical true
  Output File Name = case
  Exec Solver = Before Saving
End


kevinarden
Posts: 2237
Joined: 25 Jan 2019, 01:28
Antispam: Yes

Re: Smart Heater Control

Post by kevinarden »

I thought the smart heater control was just an on/off switch to attempt to maintain a set temperature tolerance. There is some information in the models manual on the smart heater. Not sure it is actually regulating the heat source, other than on/off.

https://www.nic.funet.fi/pub/sci/physic ... Manual.pdf
MrChips
Posts: 33
Joined: 12 Jul 2016, 00:16
Antispam: Yes
Location: Canada

Re: Smart Heater Control

Post by MrChips »

Am I supposed to be able to see it switch on and off with each timestep? I guess I figured if it was switching on and off that I would be able to see some variation in the power density as this happens. Maybe this just doesn't work so well with transient and radiation heating?
MrChips
Posts: 33
Joined: 12 Jul 2016, 00:16
Antispam: Yes
Location: Canada

Re: Smart Heater Control

Post by MrChips »

Thanks Kevin. I did see this thread, but could not find any information on the "Smart Heater Timescale" Peter refers to. I assumed that since the thread was so old that the information was antiquated.
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Smart Heater Control

Post by raback »

Hi

The smart heater is really smart in the sense that tunes the heat load exactly so that the linearized system gives the correct value for the control temperature. Now in steady-state this works perfectly. If there are nonlinearities they just add some iterations but the overall operation works fine.

Now for transient cases things are not so easy. The problem is that if you ask to exactly reproduce a given temperature that will result to an oscillating heating power. The smaller the timestep the more difficult things become.

What should be done instead is to have some more relaxed temperature control. So you follow the temperature but don't expect to meet the target exactly. Instead you have some feedback that is slow enough to not overshoot. The keyword "Smart Heater Time Scale" may be used to activate this. It really is impossible to tune a heater with timescale Dt to exactly meet the desired temperature if dt<Dt. The downside of this method is that it does not give you exactly what you ask for. Therefore this may also be lacking in documentation and testing.

If you have some equipment in your lab or factory that you need to simulate it may be best to study what kind of feedback they use in temperature control and try to mimic that with a small piece of code.

-Peter
MrChips
Posts: 33
Joined: 12 Jul 2016, 00:16
Antispam: Yes
Location: Canada

Re: Smart Heater Control

Post by MrChips »

Thanks for the explanation Peter.

In the lab setting I use an on/off heater control that uses a low and high temperature setpoint. I used a UDF in prior runs to turn the heater on and off with these temperature setpoints. It works reasonably well except that the timesteps need to be small enough to catch when it's hit the low/high temperature. If the timesteps are too large I overshoot the temperature setpoints by too much. If the timesteps are too small it just takes longer than it should to run.

I was thinking of modifying an "AdaptiveSolver" that mzenker posted back in 2012 so that I could use the max temperature as a variable to save or throw out a timestep and modify the size of a timestep. Something like this:

Code: Select all

SUBROUTINE AdaptiveSolver(Model, Solver, dt, Transient)

	USE DefUtils
	IMPLICIT NONE
	TYPE(Solver_t) :: Solver
	TYPE(Model_t) :: Model
	REAL(KIND=dp) :: dt
	LOGICAL :: Transient
	!----------------------------------------------------------------
	! Local variables
	!----------------------------------------------------------------
	TYPE(ValueList_t), POINTER :: SolverParams
    TYPE(Mesh_t), POINTER :: Mesh
    TYPE(Variable_t), POINTER :: Var
	REAL(KIND=dp), POINTER :: WrkPntr(:)
	INTEGER :: SaveResults = 1
	REAL(KIND=dp) :: MaxTimeStepSize = 0.0_dp, MinTimeStepSize, CurrentTimeStepSize=0.0_dp, MinSaveInterval
	REAL(KIND=dp) :: Time, Last_Time = 0.0_dp, Last_Save_Time = 0.0_dp, Duration
	REAL(KIND=dp) :: SaveResults_real
	LOGICAL :: Found, Visited = .FALSE.
	LOGICAL :: AdaptiveTimeStepping, Converged, StopSimulation = .FALSE.
	INTEGER :: istat
    TYPE(Solver_t), POINTER :: iSolver
    REAL(KIND=dp), ALLOCATABLE :: xx(:,:), prevxx(:,:,:)
	INTEGER :: i, j, k, n
	INTEGER :: AdaptiveKeepSmallest, KeepTimeStepSize = 0
	INTEGER :: SteadyConverged
	
	SAVE Last_Time
	SAVE CurrentTimeStepSize
	SAVE KeepTimeStepSize
	SAVE Last_Save_Time
	SAVE Visited
	SAVE xx, prevxx
	SAVE StopSimulation

	CALL Info( 'AdaptiveSolver', '*************************************************' )

	!Read in solver parameters
	!-------------------------
	SolverParams => GetSolverParams()
	IF (.NOT. ASSOCIATED(SolverParams)) CALL FATAL('AdaptiveSolver','No Solver section found')
	
	IF ( .NOT. Visited ) THEN
		! allocate space for previous result
		n = Model % NumberOfSolvers
		j = 0
		k = 0
		DO i=1,n
			iSolver => Model % Solvers(i)
			IF ( ASSOCIATED( iSolver % Variable  % Values ) ) THEN
				IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
					j = MAX( j, SIZE( iSolver % Variable % PrevValues,2 ) )
				END IF
				k = MAX( k, SIZE( iSolver % Variable % Values ) )
			END IF
		END DO
		ALLOCATE( xx(n,k), prevxx( n,k,j ) )
		! initialize variables, allocate memory
		Mesh => Model % Meshes
		DO WHILE( ASSOCIATED(Mesh) )
			IF ( Mesh % OutputActive ) EXIT 
			Mesh => Mesh % Next
		END DO
		CALL SetCurrentMesh( Model, Mesh )
		NULLIFY(WrkPntr)
		ALLOCATE(WrkPntr(1),STAT=istat)
		IF( istat /= 0 ) CALL Fatal('AdaptiveSolver','Memory allocation error') 	
		CALL VariableAdd( Model % Variables, Mesh, Solver, 'AdaptiveTimeStepSize', 1, WrkPntr )
		NULLIFY(WrkPntr)
		ALLOCATE(WrkPntr(1),STAT=istat)
		IF( istat /= 0 ) CALL Fatal('AdaptiveSolver','Memory allocation error 3') 	
		CALL VariableAdd( Model % Variables, Mesh, Solver, 'SaveResults', 1, WrkPntr )
	END IF
	
	Time = GetTime()
	
	IF ( StopSimulation ) THEN

		! Set SaveResults variable to zero
		Var => VariableGet( Model % Variables, 'SaveResults' )
		IF(.NOT. ASSOCIATED(Var)) THEN
			CALL FATAL( 'AdaptiveSolver', 'No memory allocated for SaveResults')       
		END IF
		! It seems impossible to save an Integer, so we convert it to Real
		Var % Values(1) = 0.0_dp
	
		! End of simulation was reached at previous timestep - now send stop signal
		WRITE(Message, *)  'End of Simulation! Sending stop signal to solver...'
		CALL ListAddConstReal(Model % Simulation,'Exit Condition',1.0_dp)
		
	ELSE
	
		Converged = .TRUE.
		SaveResults = 1
		DO i=1,Model % NumberOfSolvers
			iSolver => Model % Solvers(i)
			SteadyConverged = iSolver % Variable % SteadyConverged
			IF( ( Time > 0 ) .AND. ( SteadyConverged == 0 ) ) THEN
				Converged = .FALSE.
				CALL Info('AdaptiveSolver', 'No convergence for this timestep')
				EXIT
			END IF
		END DO
		
		AdaptiveTimeStepping = ListGetLogical(SolverParams,'Convergence Adaptive Timestepping',Found)
		
		IF ( AdaptiveTimeStepping ) THEN
			MinTimestepSize = ListGetConstReal( SolverParams, 'Convergence Adaptive Min Timestep', Found )
			IF ( .NOT. Found ) MinTimeStepSize = 1e-15
			AdaptiveKeepSmallest = GetInteger(SolverParams, 'Convergence Adaptive Keep Smallest', Found )
			IF ( .NOT. Found ) AdaptiveKeepSmallest = 2

			IF( Converged .OR. CurrentTimeStepSize - MinTimeStepSize < 1e-10 ) THEN
				! save current result
				DO i=1,Model % NumberOFSolvers
					iSolver => CurrentModel % Solvers(i)
					IF ( ASSOCIATED( iSolver % Variable % Values ) ) THEN
						n = SIZE( iSolver % Variable % Values )
						xx(i,1:n) = iSolver % Variable % Values
						IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
							DO j=1,SIZE( iSolver % Variable % PrevValues,2 )
								prevxx(i,1:n,j) = iSolver % Variable % PrevValues(:,j)
							END DO
						END IF
					END IF
				END DO
			ELSE
				! - reset to result for last convergence
				DO i=1,Model % NumberOFSolvers
					iSolver => CurrentModel % Solvers(i)
					IF ( ASSOCIATED( iSolver % Variable % Values ) ) THEN
						n = SIZE(iSolver % Variable % Values)
						iSolver % Variable % Values = xx(i,1:n)
						IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
							DO j=1,SIZE( iSolver % Variable % PrevValues,2 )
							iSolver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
							END DO
						END IF
					END IF
				END DO
				
				! - reset time to previous value
				Var  => VariableGet( Model % Solver % Mesh % Variables, 'Time' )
				IF ( ASSOCIATED( Var ) )  Var % Values(1)  = Time - CurrentTimeStepSize
				Time = GetTime()
				Write(Message, *) 'Adaptive Timestepping: Time reset to ', Time
				CALL Info('AdaptiveSolver', Message)
				SaveResults = 0
			END IF
		END IF ! AdaptiveTimeStepping

		! Get duration 
		Duration = ListGetConstReal( SolverParams, 'Duration',Found)
		IF ( .NOT. Found ) Duration = 0.0

		! Get time minimum saving interval 
		MinSaveInterval = ListGetConstReal(SolverParams, 'Min Data Saving Interval',Found)
		IF ( .NOT. Found ) MinSaveInterval = 0.0_dp
		
		IF( SaveResults /= 0 ) THEN
			IF( ( Time > 1e-10 ) .AND. ( Time - Last_Save_Time < MinSaveInterval ) ) THEN
				SaveResults = 0
			ELSE
				Last_Save_Time = Time
			END IF
		END IF
				
		! Check if duration is exceeded since Time was saved the last time
		! Duration is multiplied by a factor to avoid problems due to rounding errors
	    IF ( Time - Last_Time >= Duration - 1e-10 ) THEN
		! End of simulation after this timestep! Save one last time...
			SaveResults = 1
			StopSimulation = .TRUE.
		END IF
				
		! Get time step size
		MaxTimeStepSize = ListGetConstReal(SolverParams, 'Adaptive TimeStep Size',Found)
		IF ( .NOT. Found ) MaxTimeStepSize = 0.0
			
		! Determine timestep size
		IF ( .NOT. Visited ) THEN	
			CurrentTimeStepSize = MaxTimeStepSize
		ELSE IF ( AdaptiveTimeStepping ) THEN
			IF ( Converged ) THEN
				! check keepsmallest, increase timestep size if desired
				IF ( KeepTimeStepSize > 0 ) THEN
					KeepTimeStepSize = KeepTimeStepSize - 1
				ELSE
					! try to double the timestep until the timestep set by the user is reached
					CurrentTimeStepSize = CurrentTimeStepSize * 2
					IF ( CurrentTimeStepSize > MaxTimeStepSize ) THEN
						CurrentTimeStepSize = MaxTimeStepSize
					ELSE
						WRITE( Message, * ), 'Adaptive Timestepping: Time step set to ', CurrentTimeStepSize
						CALL Info('AdaptiveSolver', Message)
					END IF
				END IF
			ELSE
				! set timestep size to half if MinTimestepSize is not reached
				IF( CurrentTimeStepSize - MinTimeStepSize > 1e-10 ) THEN
					CurrentTimeStepSize = CurrentTimeStepSize / 2
					IF ( CurrentTimeStepSize < MinTimeStepSize ) CurrentTimeStepSize = MinTimeStepSize
					WRITE( Message, * ), 'Adaptive Timestepping: Time step set to ', CurrentTimeStepSize
					CALL Info('AdaptiveSolver', Message)
				ELSE
					CALL Warn('AdaptiveSolver','Cannot set smaller time step!')
					CALL Warn('AdaptiveSolver','Minimum time step might be too big to reach convergence.')
					CALL Warn('AdaptiveSolver','Continuing anyway...')
				END IF
				! reset KeepTimeStepSize
				KeepTimeStepSize = AdaptiveKeepSmallest
			END IF
		END IF
			
		! save TimeStepSize
		Var => VariableGet( Model % Variables, 'AdaptiveTimeStepSize' )
		IF(.NOT. ASSOCIATED(Var)) CALL FATAL( 'AdaptiveSolver', 'No memory allocated for AdaptiveTimeStepSize')       
		Var % Values(1) = CurrentTimeStepSize
					
		! Save SaveResults variable
		Var => VariableGet( Model % Variables, 'SaveResults' )
		IF(.NOT. ASSOCIATED(Var)) CALL FATAL( 'AdaptiveSolver', 'No memory allocated for SaveResults')       
		! It seems impossible to save an Integer, so we convert it to Real
		SaveResults_real = REAL(SaveResults)
		Var % Values(1) = SaveResults_real

		IF ( .NOT. StopSimulation ) THEN
			WRITE(Message, *)  ' SaveResults', SaveResults
		ELSE
			WRITE(Message, *)  'This was the last timestep! Preparing end of Simulation...'
		END IF
	END IF


	CALL Info( 'AdaptiveSolver', Message, level=4)
	CALL Info( 'AdaptiveSolver', '*************************************************' )
	
	IF ( .NOT. Visited ) Visited = .TRUE.
	
END SUBROUTINE AdaptiveSolver
From: viewtopic.php?f=3&t=2674&hilit=adaptive+timestep

I have to learn some Fortran...

Is there a current way in Elmer to do this with adaptive time stepping?
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Smart Heater Control

Post by raback »

Hi

The timestep may be a function. How about checking min/max temperature and the related heating/cooling speed. Then you could come up with a simple heuristics to shorten the timestep if you're about to pass. I wrote a possible UDF for timestep size. Didn't test it, not even compiled. Just to give an idea what one could do when you know your system.

Code: Select all

FUNCTION Mydt(Model, n, t) RESULT(dt)
  USE DefUtils
  IMPLICIT None
  TYPE(Model_t) :: Model
  INTEGER :: n
  REAL(KIND=dp) :: t, dt

  TYPE(Variable_t), POINTER :: Var
  REAL(KIND=dp) :: dtprev,t0,tmax,tmaxprev,dtmax,dtmin,dtest
  INTEGER :: i,imax
  
  t0 = ListGetCReal( Model % Simulation,'Upper temperature limit')
  dtmin = ListGetCReal( Model % Simulation,'Minimum dt')
  dtmax = ListGetCReal( Model % Simulation,'Maximum dt')
  dtprev = GetTimestepSize()
  
  Var => VariableGet( Model % Mesh % Variables,'temperature')
  tmax = -HUGE(tmax)
  DO i=1,SIZE(Var % Values)
    IF(Var % Values(i) > tmax ) THEN
      imax = i
      tmax = Var % Values(i)
    END IF
  END DO
  ! previous temperature at the same location
  tmaxprev = Var % PrevValues(imax,1)
  
  ! estimated timestep for hitting the target temperature
  ! Note that this might be more intelligent. Just something for example.
  dtest = dtprev*(t0-tmax)/(tmax-tmaxprev)  
  dt = MAX(dtmin,MIN(dtmax,dtest))  
END IF    
-Peter
MrChips
Posts: 33
Joined: 12 Jul 2016, 00:16
Antispam: Yes
Location: Canada

Re: Smart Heater Control

Post by MrChips »

Thanks Peter. This seems promising. I am working through it (slowly). I can get this mostly to run except for when I uncomment the GetTimestepSize() line it gives a segmentation fault. Is there a function header line I need to put in for this to work? As far as I can tell it should just return an integer. Same thing if I try to use GetTime().
Post Reply