/****************************************************************
 * FUNCTION wrksim.c -  evaluate the pressure drop work loss of *
 *         non-ideal heat exchangers                            *
 * FILENAME:       hotsim.c                                     *
 * PROGRAMMERS:    Izzi Urieli (FORTRAN)                        *
 *                 Eric Malroy (translation to "C")             *
 * DATE:           122091/082096                                *
 * LAST MODIFICATION: 09/05/96 by I.Urieli                      *
 * INCLUDE: <math.h>, <stdio.h>, "define.h","adiabatic.h",      *
 *          "simple.h"                                          *
 * GLOBAL VARIABLES:                                            *
 *      tk - cold sink operation temperature (K)                *
 *      dk - cooler hydraulic diameter (m)                      *
 *      omega - operating frequency (rad/sec)                   *
 *      ak - cooler internal free flow area  (sq.m)             *
 *      vk - cooler void volume (cu.m)                          *
 *      lk - cooler effective length (m)                        *
 *      ar - regenerator internal free flow area (sq.m)         *
 *      tr - effective regenerator temperature (K)              *
 *      dr - regenerator hydraulic diameter (m)                 *
 *      vr - regenerator void volume (cu.m)                     *
 *      lr - regenerator effective length (m)                   *
 *      ah - heater internal free flow area (sq.m)              *
 *      th - hot sink operation temperature (K)                 *
 *      dh - heater hydraulic diameter (m)                      *
 *      vh - heater void volume (cu.m)                          *
 *      lh - heater effective length (m)                        *
 *      rmatrx - char describing regenerator type ("f" - foil   *
 *               & "m" - mesh)                                  *
 *      printfile - address to the print output file            *
 * PROTOTYPE: void wrksim(double var[][COL],double dvar[][COL], *
 *            double *pt_dwork)                                 *
 * PRE:  tk,dk,omega,ak,vk,lk,ar,tr,dr,vr,lr,ah,th,dh,vh,lh,    *
 *       printfile,rmatrx                                       *
 * POST: *pt_dwork is fully determined                          *
 ****************************************************************/
#include <math.h>
#include <stdio.h>
#include "define.h"
#include "adiabatic.h"
#include "simple.h"

void wrksim(double var[][COL],  /* matrix of output derivatives     */
            double dvar[][COL], /* matrix of output variables       */
            double *pt_dwork)   /* pressure drop work loss (Joules) */
{
   const int kol = 0;
   const int reg = 1;
   const int hot = 2;

   double delp[4][COL]; /* differential pressure for heater,        */
                        /* regen and cooler (array 3 X 37)          */
   double dp[COL];  /* Total differential pressure (Pa)             */
   double gk; /* mass flow rate in cooler (kg/s)                    */
   double gr; /* mass flow rate in regenerator (kg/s)               */
   double gh; /* mass flow rate in heater (kg/s)                    */
   double mu; /* dynamic viscosity (kg/(m*s))                       */
   double re; /* Reynolds number   (fabs(g)*d/mu)                   */
   double st; /* Stanton number                                     */
   double ht; /* heat transfer coefficient (W/(m^2*K))              */
   double fr; /* Reynolds friction factor (Ff*Re)                   */
              /* (Ff is the Fanning friction factor)                */
   double kgas; /* gas thermal conductivity (W/(m^2*K))             */
   double dwork; /* pressure drop work loss (Joules)                */
   double dtheta; /* incremental step of the crank angle            */
   int i; /* index for "for" loop                                   */
   int angle; /* crank angle (degrees)                              */

   dtheta = 2.0*PI/36.0;
   dwork = 0.0;
   for(i=0; i<36; i++) {
      gk = (var[GACK][i] + var[GAKR][i])*omega/(2.0*ak);
      reynum(tk,gk,dk,&mu,&kgas,&re);
      pipfr(dk,mu,re,&ht,&fr);
      delp[kol][i] = 2.0*fr*mu*vk*gk*lk/(var[MK][i]*dk*dk);
      gr = (var[GAKR][i] + var[GARH][i])*omega/(2.0*ar);
      reynum(tr,gr,dr,&mu,&kgas,&re);
      if(rmatrx == 'm')
         matfr(re,&st,&fr);
      else if(rmatrx == 'f')
         foilfr(dr,mu,re,&st,&ht,&fr);
      else {
         printf("matrix %c is undefined\n", rmatrx);
         fprintf(printfile,"matrix %c is undefined\n", rmatrx);
      }
      delp[reg][i] = 2.0*fr*mu*vr*gr*lr/(var[MR][i]*dr*dr);
      gh = (var[GARH][i] + var[GAHE][i])*omega/(2.0*ah);
      reynum(th,gh,dh,&mu,&kgas,&re);
      pipfr(dh,mu,re,&ht,&fr);
      delp[hot][i] = 2.0*fr*mu*vh*gh*lh/(var[MH][i]*dh*dh);
      dp[i] = delp[kol][i] + delp[reg][i] + delp[hot][i];
      dwork = dwork + dtheta*dp[i]*dvar[VE][i];
  }
  delp[kol][COL-1] = delp[kol][0];
  delp[reg][COL-1] = delp[reg][0];
  delp[hot][COL-1] = delp[hot][0];
  *pt_dwork = dwork;
  fprintf(printfile,
           "\ntheta Pc [Pa]    dpk [Pa]   dpr [Pa]  dph [Pa]  Pe [Pa]\n");
  printf("\ntheta Pc [Pa]    dpk [Pa]   dpr [Pa]  dph [Pa]  Pe [Pa]\n ");
  for(i=0;i<COL;i++) {
     angle=i*10;
     printf("%5d%10.3e%10.f%10.f%10.f%10.3e\n",angle,var[P][i],delp[kol][i],
             delp[reg][i],delp[hot][i],(var[P][i]+dp[i]));
     fprintf(printfile,"%5d%10.3e%10.f%10.f%10.f%10.3e\n",angle,var[P][i],
             delp[kol][i],delp[reg][i],delp[hot][i],(var[P][i]+dp[i]));
  }
}