1  /******************************************************************
     2   * FUNCTION: adiab - Stirling cycle machine simulation            *
     3   *                   Ideal Adiabatic analysis                     *
     4   * FILENAME adiab.c                                               *
     5   * PROGRAMMER I.Urieli    (FORTRAN)                               *
     6   *            E. Malroy   (C Translation)                         *
     7   * DATE:     11/30/94  (FORTRAN), 03/12/96 (C Translation)        *
     8   * LAST MODIFIED: 02/02/97 by I.Urieli                            *
     9   * INCLUDE:  <stdio.h>,<math.h>,"adiabatic.h","define.h"          *
    10   * GLOBAL VARIABLES: none                                         *
    11   * PRE:   all global variables in define.h have values            *
    12   * POST:  var, dvar matrices have been filled in                  *
    13   ******************************************************************/
    14  #include <stdio.h>
    15  #include <math.h>
    16  #include "adiabatic.h" 
    17  #include "define.h" 
    18  #include "odes.h"
    19  
    20  void adiab(double var[][COL],  /* matrix of variables            */
    21             double dvar[][COL]) /* matrix of derivatives          */
    22  {
    23  
    24    
    25     const double epsilon = 1.0; /* allowable error in temperature (K) */
    26     const int maxit = 20; /* maximum number of iterations to converge */
    27  
    28  /*    variables  */
    29  
    30     int i,j, iter, ninc, step;
    31     double y[ROWV];   /* dependent variable vector */
    32     double dy[ROWD];    /* dependent derivative vector */
    33     double tc0;       /* initial cycle value comp temp */
    34     double te0;       /* initial cycle value exp temp  */
    35     double terror;    /* cyclic accumulative temp error */
    36     double theta;     /* crankangle - independent variable   */
    37     double dtheta;    /* increment of crankangle */
    38  
    39     printf("\n...now doing an Ideal Adiabatic simulation...\n");
    40     fprintf(printfile,"\nIdeal Adiabatic simulation:\n");
    41     ninc = 360;
    42     step = ninc/36;
    43     dtheta = 2.0*PI/ninc;
    44     y[THE] = th;
    45     y[TCK] = tk;
    46     y[TE] = th;
    47     y[TC] = tk;
    48     iter = 1;
    49  /*     *** solve **************************************************/
    50     do {
    51        printf("\niteration %2d: initial Tc = %.1f(K), Te = %.1f(K)\n",
    52             iter,y[TC],y[TE]);
    53        fprintf(printfile,"\niteration %2d: initial Tc = %.1f(K), Te = %.1f(K)\n",
    54             iter,y[TC],y[TE]);
    55        tc0 = y[TC];
    56        te0 = y[TE];
    57        theta = 0.0;
    58        y[QK] = 0.0;
    59        y[QR] = 0.0;
    60        y[QH] = 0.0;
    61        y[WC] = 0.0;
    62        y[WE] = 0.0;
    63        y[W] = 0.0;
    64        for(i=0; i<ninc; i++) 
    65           rk4(dadiab,7,&theta,dtheta,y,dy);
    66        terror = fabs(tc0 - y[TC]) + fabs(te0 - y[TE]);
    67        printf("              final Tc = %.1f(K), Te = %.1f(K)\n", y[TC], y[TE]);
    68        printf("   temperature error (del(Tc) + del(Te)) = %.3f(K)\n", terror);
    69        fprintf(printfile, "              final Tc = %.1f(K), Te = %.1f(K)\n",
    70                 y[TC], y[TE]);
    71        fprintf(printfile, "   temperature error (del(Tc) + del(Te)) = %.3f(K)\n",
    72                 terror);
    73        iter++;
    74     }
    75     while((terror >= epsilon) && (iter <= maxit));
    76     if(iter > maxit) {
    77        printf("\n no convergence within %d iterations\n", maxit);
    78        fprintf(printfile,"\n no convergence within %d iterations\n", maxit);
    79     }        
    80  /*     *** setup output matrices **************************************/
    81     theta = 0.0;
    82     y[QK] = 0.0;
    83     y[QR] = 0.0;
    84     y[QH] = 0.0;
    85     y[WC] = 0.0;
    86     y[WE] = 0.0;
    87     y[W] = 0.0;
    88     
    89     filmat(0,y,dy,var,dvar);
    90     for(i=1; i<COL; i++) {   
    91        for(j=0; j<step; j++)
    92           rk4(dadiab,7,&theta,dtheta,y,dy);
    93        filmat(i,y,dy,var,dvar);
    94     }
    95  }
condor{t.gentry}4: