/*
* 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"
/*
* Discretize the fluid and pipe layers
* Calculate interface conductivity km and kp
* Calculate cells surface
*/
void MeshPipe(PipeProblem *PP, PipeDicretize *PZ)
{
int l ;
int i ;
/*
* For each layer
*/
for( l = 0 ; l <= PP->N ; l++ )
{
/*
* Special case : the Pipe/fluid (layer 0)
*/
if (l == 0)
{
double HalfInterval = PP->r[l] / ( 2.0 * (double)PP->NbCell[l] - 1.0 ) ;
Size0PipeDiscretize(PZ) ;
PZ->L[PZ->N] = l ;
PZ->rm[PZ->N] = 0.0 ;
PZ->rc[PZ->N] = 0.0 ;
PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;
PZ->Cp[PZ->N] = PP->Cp[l] ;
PZ->d[PZ->N] = PP->d[l] ;
PZ->kc[PZ->N] = PP->k[l] ;
/*
* Dicretize the others intervals
*/
for( i = 2 ; i <= PP->NbCell[l] ; i++ )
{
IncrementSizePipeDiscretize(PZ) ;
PZ->L[PZ->N] = l ;
PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;
PZ->Cp[PZ->N] = PP->Cp[l] ;
PZ->d[PZ->N] = PP->d[l] ;
PZ->kc[PZ->N] = PP->k[l] ;
}
}
/*
* Special case : the last layer
*/
else if ( l == PP->N )
{
double HalfInterval = ( PP->r[l] - PP->r[l-1] ) / ( 2.0 * (double)PP->NbCell[l] - 1.0 ) ;
for( i = 2 ; i <= PP->NbCell[l] ; i++ )
{
IncrementSizePipeDiscretize(PZ) ;
PZ->L[PZ->N] = l ;
PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;
PZ->Cp[PZ->N] = PP->Cp[l] ;
PZ->d[PZ->N] = PP->d[l] ;
PZ->kc[PZ->N] = PP->k[l] ;
}
IncrementSizePipeDiscretize(PZ) ;
PZ->L[PZ->N] = l ;
PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
PZ->rp[PZ->N] = PZ->rc[PZ->N] ;
PZ->Cp[PZ->N] = PP->Cp[l] ;
PZ->d[PZ->N] = PP->d[l] ;
PZ->kc[PZ->N] = PP->k[l] ;
}
/*
* Normal case: the other layers
*/
else
{
double HalfInterval = ( PP->r[l] - PP->r[l-1] ) / ( 2.0 * (double)PP->NbCell[l] ) ;
for( i = 1 ; i <= PP->NbCell[l] ; i++ )
{
IncrementSizePipeDiscretize(PZ) ;
PZ->L[PZ->N] = l ;
PZ->rm[PZ->N] = PZ->rp[PZ->N - 1] ;
PZ->rc[PZ->N] = PZ->rm[PZ->N] + HalfInterval ;
PZ->rp[PZ->N] = PZ->rc[PZ->N] + HalfInterval ;
PZ->Cp[PZ->N] = PP->Cp[l] ;
PZ->d[PZ->N] = PP->d[l] ;
PZ->kc[PZ->N] = PP->k[l] ;
}
}
}
/*
* Calculate k[i-1/2] and k[i+1/2] conductivity (resp. km and kp)
*/
i = 0 ;
PZ->km[i] = PZ->kc[i] ;
PZ->kp[i] = PZ->kc[i] ;
for( i = 1 ; i <= (PZ->N-1) ; i++ )
{
PZ->km[i] = PZ->kp[i-1] ;
PZ->kp[i] = PZ->kc[i] * PZ->kc[i+1] * log(PZ->rc[i+1]/PZ->rc[i])
/ ( PZ->kc[i+1] * log(PZ->rp[i]/PZ->rc[i]) + PZ->kc[i] * log(PZ->rc[i+1]/PZ->rp[i]) ) ;
}
i = PZ->N ;
PZ->km[i] = PZ->kp[i-1] ;
PZ->kp[i] = PZ->kc[i] ;
/*
* Calculate V[i] = Cells Volume / unit length = Surface
*/
for( i = 0 ; i <= PZ->N ; i++ )
PZ->S[i] = C_Pi * ( PZ->rp[i]*PZ->rp[i] - PZ->rm[i]*PZ->rm[i] ) ;
}
/*
* Alpha, Beta & Gamma coefficients calculation
*/
void CalcModelCoefs(PipeProblem *PP, PipeDicretize *PZ)
{
/*
* At the pipe center, i=0
*/
int i = 0 ;
PZ->alpha[i] = 0.0 ;
PZ->beta[i] = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
PZ->gamma[i] = 2.0 * C_Pi * PZ->rp[i] * PZ->kp[i] / ( PZ->rc[i+1] - PZ->rc[i] ) ;
/*
* In the pipe layers, 2<=i<=N-1
*/
for( i = 1 ; i < PZ->N ; i++ )
{
PZ->alpha[i] = PZ->gamma[i-1] ;
PZ->beta[i] = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
PZ->gamma[i] = 2.0 * C_Pi * PZ->rp[i] * PZ->kp[i] / ( PZ->rc[i+1] - PZ->rc[i] ) ;
}
/*
* At the ambiant fluid / outer wall interface, i=N
*/
i = PZ->N ;
PZ->alpha[i] = PZ->gamma[i-1] ;
PZ->beta[i] = PZ->S[i] * PZ->d[i] * PZ->Cp[i] ;
PZ->gamma[i] = 2.0 * C_Pi * PZ->rc[i] * PP->Hout ;
/*
* AlphaSharp Coef for fluid/solid temperature calculation
*/
i = PP->NbCell[0] ;
PZ->alphaSharp = 2.0 * C_Pi * PZ->rc[i] * PZ->kc[i] / ( PZ->rc[i] - PZ->rm[i] ) ;
}
/*
* Build the Pipe discrete model
*/
void BuildModel(PipeProblem *PP, PipeDicretize *PZ)
{
int l ;
/*
* Apply discretization limitations
*/
/* For each layer */
for( l = 0 ; l <= PP->N ; l++ )
{
if ( PP->NbCell[l] < 1 ) PP->NbCell[l] = 1 ;
if ( PP->NbCell[l] > MaxCellNumberPerLayer ) PP->NbCell[l] = MaxCellNumberPerLayer ;
}
/* For time step */
if (PP->dt < MinTimeStep) PP->dt = MinTimeStep ;
/*
* Mesh the pipe and calculate alpha, beta & gamma coefs
*/
MeshPipe(PP, PZ) ;
CalcModelCoefs(PP, PZ) ;
}
/* End of mesh.c */