Cooldt Code

/*
 * cooldt: A cool-down time calculator for multi-layers pipe
 *
 * Copyright 2010, 2013 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 and gamma for the steady state problem
 */
void SteadyStateSystem(PipeProblem *PP, PipeDicretize *PZ,
                       double *A, double *B, double *C, double *D)
{
    int i ;

    /*
     * Set unused coefs. to 0
     */
    for( i = 0 ; i < PP->NbCell[0]  ; i++ )
    {
        A[i] = 0.0 ;
        C[i] = 0.0 ;
        B[i] = 0.0 ;
        D[i] = 0.0 ;
    }

    /*
     * First coef. corresponding to Tinitial temperature at fluid/solid interface
     */
    i = PP->NbCell[0] ;
    A[i] = 0.0 ;
    B[i] = PZ->alphaSharp + PZ->gamma[i] ;
    C[i] = - PZ->gamma[i] ;
    D[i] = PZ->alphaSharp * PP->Tinitial ;

    /*
     * Solid layers coefs.
     */
    for( i = (PP->NbCell[0] + 1)  ; i < PZ->N ; i++ )
    {
        A[i] = - PZ->alpha[i] ;
        B[i] = PZ->alpha[i] + PZ->gamma[i] ;
        C[i] = - PZ->gamma[i] ;
        D[i] = 0.0 ;
    }

    /*
     * At the ambiant fluid / outer wall interface
     */
    i = PZ->N ;
    A[i] = - PZ->alpha[i] ;
    B[i] = PZ->alpha[i] + PZ->gamma[i] ;
    C[i] = 0.0 ;
    D[i] = PZ->gamma[i] * PP->Tambient ;

}


/*
 * Calculate steady state profile temperature
 * return the result to *T
 */
void CalcSteadyState(PipeProblem *PP, PipeDicretize *PZ, double *T)
{
    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) ;
    int i;

    /*
     * Impose the initial fluid temperature to fluid cells
     */
    for ( i = 0 ; i < PP->NbCell[0] ; i++ )
        T[i] = PP->Tinitial ;

    /*
     * Calculate the system coeficient for the steady case state
     */
    SteadyStateSystem(PP, PZ, A, B, C, D) ;

    /*
     * Solve the system for solid material ( PP->NbCell[0] = NumberOfFluidCells )
     */
    TriDiagonalSolve(A, B, C, D, PP->NbCell[0], PZ->N, T) ;

    /*
     * Free local pointers
     */
    free(A) ;
    free(B) ;
    free(C) ;
    free(D) ;
}

/* End of steady-state.c*/