/*
* Cooldt cool-down time calculator for multi-layers pipe
*
* Copyright 2010, 2013, 2014, Benjamin DEGLO DE BESSES. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are
* permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this list of
* conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright notice, this list
* of conditions and the following disclaimer in the documentation and/or other materials
* provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY Benjamin DEGLO DE BESSES "AS IS" AND ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
* FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Benjamin DEGLO DE BESSES OR
* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
* ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*/
#include "cooldt.h"
/*
* Calculate linear system coefficients A, B, C & D
* based on alpha, beta, gamma and time step (dt)
*/
void TransientSystem(PipeProblem *PP, PipeDicretize *PZ, double *T, double dt,
double *A, double *B, double *C, double *D)
{
/*
* At the pipe center
*/
int i = 0 ;
A[i] = 0.0 ;
B[i] = PZ->beta[i] + dt * PZ->gamma[i] ;
C[i] = - dt * PZ->gamma[i] ;
D[i] = PZ->gamma[i] * ( T[i+1] - T[i] ) ;
/*
* In fluid and solids
*/
for( i = 1 ; i < PZ->N ; i++ )
{
A[i] = - dt * PZ->alpha[i];
B[i] = PZ->beta[i] + dt * ( PZ->alpha[i] + PZ->gamma[i] ) ;
C[i] = - dt * PZ->gamma[i] ;
D[i] = PZ->alpha[i] * ( T[i-1] - T[i] ) + PZ->gamma[i] * ( T[i+1] - T[i] ) ;
}
/*
* At the ambiant fluid / outer wall interface
*/
i = PZ->N ;
A[i] = - dt * PZ->alpha[i];
B[i] = PZ->beta[i] + dt * ( PZ->alpha[i] + PZ->gamma[i] ) ;
C[i] = 0.0 ;
D[i] = PZ->alpha[i] * ( T[i-1] - T[i] ) + PZ->gamma[i] * ( PP->Tambient - T[i] ) ;
}
/*
* Calculate temperature change derivative / time for each cell (=dT/dt)
* return the result to *DTodt
*/
void CalcDTodt(PipeProblem *PP, PipeDicretize *PZ, double *T, double dt, double *DTodt)
{
size_t DoublesSize = (PZ->N + 1) * sizeof(double) ;
double *A = (double*)malloc(DoublesSize) ;
double *B = (double*)malloc(DoublesSize) ;
double *C = (double*)malloc(DoublesSize) ;
double *D = (double*)malloc(DoublesSize) ;
/*
* Calculate coef, then the next temperature
*/
TransientSystem(PP, PZ, T, dt, A, B, C, D) ;
TriDiagonalSolve(A, B, C, D, 0, PZ->N, DTodt) ;
free(A) ;
free(B) ;
free(C) ;
free(D) ;
}
/*
* Calculate next temperature with exponential decay approximation
*/
void TimeIntegration(PipeProblem *PP, PipeDicretize *PZ, double *T)
{
int i ;
double *DTodt = (double*)malloc((PZ->N + 1) * sizeof(double)) ;
double Tau ;
/*
* Step 1: Calculate DT/dt
*/
CalcDTodt(PP, PZ, T, PP->dt, DTodt) ;
/*
* Step 2: Calculate Tau[i] then new T[i]
*/
for( i = 0 ; i <= PZ->N ; i++ )
{
/*
* Time constant for each cell temperature decay
*/
Tau = - ( T[i] + PP->dt * DTodt[i] - PP->Tambient ) / DTodt[i];
/*
* Apply the exponential decay approximation
*/
T[i] = PP->Tambient + ( T[i] - PP->Tambient ) * exp( - PP->dt / Tau ) ;
}
/*
* Free local pointers
*/
free(DTodt) ;
}
/*
* Return the temperature at the fluid/solid interface
*/
double TfluidSolidInterface(PipeDicretize *PZ, int NbFluidCell, double *T)
{
return ( ( ( PZ->alphaSharp + PZ->gamma[NbFluidCell] ) * T[NbFluidCell]
- PZ->gamma[NbFluidCell] * T[NbFluidCell+1] ) / PZ->alphaSharp ) ;
}
/*
* Calculate the transient problem solution
* return the cool down time
*/
double CalcTransient(PipeProblem *PP, PipeDicretize *PZ, double *Tcurrent, TfluidSolid *TfS)
{
double *Tprev = (double*)malloc((PZ->N + 1) * sizeof(double)) ; /* Previous temperature */
double TfScurrent ;
double time = 0.0 ;
int i ;
/*
* Save the initial temperature at fluid/solid interface
*/
Size0TfluidSolid(TfS) ;
TfS->tn[TfS->n] = time ;
TfScurrent = TfS->T[TfS->n] = TfluidSolidInterface(PZ, PP->NbCell[0], Tcurrent) ;
/*
* While the temperature have not reach the target
*/
while ( TfScurrent >= PP->Tfinal )
{
/*
* Save the previous temperature
*/
for( i = 0 ; i <= PZ->N ; i++)
Tprev[i] = Tcurrent[i] ;
/*
* Calculate next temperature
*/
TimeIntegration(PP, PZ, Tcurrent) ;
/*
* Increment time
*/
time += PP->dt ;
/*
* Calculate temperature at fluid/solid interface
*/
TfScurrent = TfluidSolidInterface(PZ, PP->NbCell[0], Tcurrent) ;
/*
* Check if the temperature change since the last save is significant
*/
if ( ((TfS->T[TfS->n] - TfScurrent ) > ((PP->Tinitial - PP->Tfinal) / NbTfluidSolidSave))
&& (TfScurrent >= PP->Tfinal) )
{
/*
* It's significant: Save the actual temperature at fluid/solid interface
*/
IncrementSizeTfluidSolid(TfS) ;
TfS->tn[TfS->n] = time ;
TfS->T[TfS->n] = TfScurrent ;
}
}
/*
* Linear interpolation to retreive the CDT at the exact target temperature
*/
double x = ( PP->Tfinal - TfScurrent )
/ ( TfluidSolidInterface(PZ, PP->NbCell[0], Tprev) - TfScurrent ) ;
time -= PP->dt * x ;
/*
* Do the same for the temperature profile
*/
for( i = 0 ; i <= PZ->N ; i++)
Tcurrent[i] += x * ( Tprev[i] - Tcurrent[i] ) ;
/*
* Save the last temperature at fluid/solid interface
*/
IncrementSizeTfluidSolid(TfS) ;
TfS->tn[TfS->n] = time ;
TfS->T[TfS->n] = TfluidSolidInterface(PZ, PP->NbCell[0], Tcurrent) ;
/*
* Free local pointers
*/
free(Tprev) ;
/*
* Return the cool down time
*/
return(time) ;
}
/* End of transient.c */