From 55577aee1ad7b242c65886d1e019d1d67ba471c8 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Thu, 27 Sep 2007 22:35:14 +0000 Subject: [PATCH] initial import --- apps/bvp/AxiStagnBVP.cpp | 144 ++++++++++++++++++++++++++++++++ apps/bvp/AxiStagnBVP.h | 131 +++++++++++++++++++++++++++++ apps/bvp/BoundaryValueProblem.h | 1 + apps/bvp/stagnation.mak | 86 +++++++++++++++++++ 4 files changed, 362 insertions(+) create mode 100644 apps/bvp/AxiStagnBVP.cpp create mode 100644 apps/bvp/AxiStagnBVP.h create mode 100644 apps/bvp/stagnation.mak diff --git a/apps/bvp/AxiStagnBVP.cpp b/apps/bvp/AxiStagnBVP.cpp new file mode 100644 index 000000000..36c8f4f16 --- /dev/null +++ b/apps/bvp/AxiStagnBVP.cpp @@ -0,0 +1,144 @@ +/// @file AxiStagnBVP.cpp + +#include +#include "AxiStagnBVP.h" + +AxiStagnBVP::AxiStagnBVP(int nsp, int np, double L) : + BVP::BoundaryValueProblem(nsp+4, + np, 0.0, L) { + + // specify the component bounds, error tolerances, and names. + BVP::Component u; + u.lower = -200.0; + u.upper = 200.0; + u.rtol = 1.0e-8; + u.atol = 1.0e-15; + u.name = "u"; + setComponent(0, u); // the axial velocity will be component 0 + + BVP::Component V; + V.lower = -1.0e8; + V.upper = 1.0e8; + V.rtol = 1.0e-8; + V.atol = 1.0e-15; + V.name = "V"; + setComponent(1, V); // the radial velocity will be component 1 + + BVP::Component T; + T.lower = 200.0; + T.upper = 1.0e9; + T.rtol = 1.0e-8; + T.atol = 1.0e-15; + T.name = "T"; + setComponent(2, T); // the temperature will be component 2 + + BVP::Component lambda; + lambda.lower = -1.0e20; + lambda.upper = 1.0e20; + lambda.rtol = 1.0e-8; + lambda.atol = 1.0e-15; + lambda.name = "Lambda"; + setComponent(3, lambda); // the pressure-gradient eigenvalue will be + //component 3 + BVP::Component Y; + Y.lower = -1.0e-5; + Y.upper = 1.0e2; + Y.rtol = 1.0e-8; + Y.atol = 1.0e-15; + for (k = 0; k < nsp; k++) { + Y.name = thermo->speciesName(k); + setComponent(k+4, Y); + } +} + + +// destructor +AxiStagnBVP::~AxiStagnBVP() {} + +// specify guesses for the initial values. These can be anything +// that leads to a converged solution. +doublereal initialValue(int n, int j) { + switch (n) { + case 0: + return m_u0; + case 1: + return m_u0/m_L; + case 2: + return m_Tinf; + case 4: + return 1.0; + default: + return 0.0; + } +} + +// Specify the residual. This is where the ODE system and boundary +// conditions are specified. The solver will attempt to find a solution +// x so that this function returns 0 for all n and j. +doublereal AxiStagnFlow::residual(doublereal* x, int n, int j) { + + // if n = 0, return the residual for the continuity equation + if (n == 0) { + if (isRight(j)) { + return -rho_u(x,j); // force u to zero at the right + } + else { + return -(rho_u(x, j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } + } + + else if (n == 1) { + + // if n = 1, then return the residual for radial momentum + if (isLeft(j)) { + return V(x,j); + } + else if (isRight(j)) { + return V(x,j); // force V to zero at the wall + } + else { + return (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + } + + } + + else if (n == 2) { + if (isLeft(j)) { + return T(x,j) - m_Tinlet; + } + else if (isRight(j)) { + return T(x,j) - m_Tsurf; + } + else { + setGas(x,j); + + // heat release term + const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); + const vector_fp& cp_R = m_thermo->cp_R_ref(); + + sum = 0.0; + sum2 = 0.0; + doublereal flxk; + for (k = 0; k < m_nsp; k++) { + flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += wdot(k,j)*h_RT[k]; + sum2 += flxk*cp_R[k]/m_wt[k]; + } + sum *= GasConstant * T(x,j); + dtdzj = dTdz(x,j); + sum2 *= GasConstant * dtdzj; + + rsd = - m_cp[j]*rho_u(x,j)*dtdzj + - divHeatFlux(x,j) - sum - sum2; + rsd /= (m_rho[j]*m_cp[j]); + + rsd -= rdt*(T(x,j) - T_prev(j)); + } + } + +} + + diff --git a/apps/bvp/AxiStagnBVP.h b/apps/bvp/AxiStagnBVP.h new file mode 100644 index 000000000..fc69f4c5a --- /dev/null +++ b/apps/bvp/AxiStagnBVP.h @@ -0,0 +1,131 @@ +/// @file AxiStagnBVP.h + +#include +#include "BoundaryValueProblem.h" + + +/** + * This class solves + */ +class AxiStagnBVP : public BVP::BoundaryValueProblem { + +public: + + AxiStagnBVP(int nsp, int np, double L) : BVP::BoundaryValueProblem(nsp+4, + np, 0.0, L) { + + // specify the component bounds, error tolerances, and names. + BVP::Component u; + u.lower = -200.0; + u.upper = 200.0; + u.rtol = 1.0e-8; + u.atol = 1.0e-15; + u.name = "u"; + setComponent(0, u); // the axial velocity will be component 0 + + BVP::Component V; + V.lower = -1.0e8; + V.upper = 1.0e8; + V.rtol = 1.0e-8; + V.atol = 1.0e-15; + V.name = "V"; + setComponent(1, V); // the radial velocity will be component 1 + + BVP::Component T; + T.lower = 200.0; + T.upper = 1.0e9; + T.rtol = 1.0e-8; + T.atol = 1.0e-15; + T.name = "T"; + setComponent(2, T); // the temperature will be component 2 + + BVP::Component lambda; + lambda.lower = -1.0e20; + lambda.upper = 1.0e20; + lambda.rtol = 1.0e-8; + lambda.atol = 1.0e-15; + lambda.name = "Lambda"; + setComponent(3, lambda); // the pressure-gradient eigenvalue will be + //component 3 + BVP::Component Y; + Y.lower = -1.0e-5; + Y.upper = 1.0e2; + Y.rtol = 1.0e-8; + Y.atol = 1.0e-15; + for (k = 0; k < nsp; k++) { + Y.name = thermo->speciesName(k); + setComponent(k+4, Y); + } + } + + + // destructor + virtual ~AxiStagnBVP() {} + + // specify guesses for the initial values. These can be anything + // that leads to a converged solution. + virtual doublereal initialValue(int n, int j) { + switch (n) { + case 0: + return 0.1*z(j); + case 1: + return 0.5*z(j); + default: + return 0.0; + } + } + + // Specify the residual. This is where the ODE system and boundary + // conditions are specified. The solver will attempt to find a solution + // x so that this function returns 0 for all n and j. + virtual doublereal residual(doublereal* x, int n, int j) { + + // if n = 0, return the residual for the first ODE + if (n == 0) { + if (isLeft(j)) // here we specify zeta(0) = 0 + return zeta(x,j); + else + // this implements d(zeta)/dz = u + return (zeta(x,j) - zeta(x,j-1))/(z(j)-z(j-1)) - u(x,j); + } + // if n = 1, then return the residual for the second ODE + else { + if (isLeft(j)) // here we specify u(0) = 0 + return u(x,j); + else if (isRight(j)) // and here we specify u(L) = 1 + return u(x,j) - 1.0; + else + // this implements the 2nd ODE + return cdif2(x,1,j) + 0.5*zeta(x,j)*centralFirstDeriv(x,1,j); + } + } + + +private: + + // for convenience only. Note that the compiler will inline these. + double zeta(double* x, int j) { return value(x,0,j); } + double u(double* x, int j) { return value(x,1,j); } + +}; + + +int main() { + try { + + // Specify a problem on (0,10), with an initial uniform grid of + // 6 points. + AxiStagnBVP eqs(6, 10.0); + // Solve the equations, refining the grid as needed, and print lots of diagnostic output (loglevel = 4) + eqs.solve(4); + // write the solution to a CSV file. + eqs.writeCSV(); + return 0; + } + catch (CanteraError) { + showErrors(cerr); + return -1; + } +} + + diff --git a/apps/bvp/BoundaryValueProblem.h b/apps/bvp/BoundaryValueProblem.h index 6083fd42f..106aabc4b 100644 --- a/apps/bvp/BoundaryValueProblem.h +++ b/apps/bvp/BoundaryValueProblem.h @@ -8,6 +8,7 @@ #include #include using namespace Cantera; +using namespace std; /// Namespace for the boundary value problem package. namespace BVP { diff --git a/apps/bvp/stagnation.mak b/apps/bvp/stagnation.mak new file mode 100644 index 000000000..a6629cd2a --- /dev/null +++ b/apps/bvp/stagnation.mak @@ -0,0 +1,86 @@ +#!/bin/sh + +# This Makefile builds a C++ application that uses Cantera. By +# default, the main program file is 'demo.cpp,' which prints out some +# properties of a reacting gas mixture. + +# To build program 'demo', simply type 'make', or 'make -f ' if this file is named something other than 'Makefile.' + +# Once you have verified that the demo runs, edit this file to replace +# object file 'demo.o' with your own object file or files. + + +#------------------------ edit this block --------------------------------- + +# the name of the executable program to be created +PROG_NAME = stagnation.x + +# the object files to be linked together. +OBJS = AxiStagnBVP.o + +# additional flags to be passed to the linker. If your program +# requires other external libraries, put them here +LINK_OPTIONS = -L/usr/local/lib -framework Accelerate + +#--------------------------------------------------------------------------- +# You probably don't need to edit anything below. + +# the C++ compiler +CXX = g++ + +# C++ compile flags +CXX_FLAGS = -O3 -Wall -fPIC + +# external libraries +EXT_LIBS = -luser -loneD -lzeroD -lequil -lkinetics -ltransport -lthermo -lctnumerics -lcvode -lctbase -lctmath -ltpx -lctf2c -lconverters -lctcxx + +# Ending C++ linking libraries +LCXX_END_LIBS = -lm + +# the directory where the Cantera libraries are located +CANTERA_LIBDIR=/Applications/Cantera/lib + +# the directory where Cantera include files may be found. +CANTERA_INCDIR=/Applications/Cantera/include + +# flags passed to the C++ compiler/linker for the linking step +LCXXFLAGS = -L$(CANTERA_LIBDIR) -O3 -Wall -fPIC + +# how to compile C++ source files to object files +.cpp.o: + $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + +PROGRAM = $(PROG_NAME)$(EXE_EXT) + +DEPENDS = $(OBJS:.o=.d) + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + $(CXX) -o $(PROGRAM) $(OBJS) $(LCXXFLAGS)\ + $(CANTERA_LIBS) $(LINK_OPTIONS) $(EXT_LIBS) \ + $(LCXX_END_LIBS) + +%.d: + g++ -MM -I$(CANTERA_INCDIR) $*.cpp > $*.d + +clean: + $(RM) $(OBJS) $(PROGRAM) + +depends: $(DEPENDS) + cat *.d > .depends + $(RM) $(DEPENDS) + +TAGS: + etags *.h *.cpp + +ifeq ($(wildcard .depends), .depends) +include .depends +endif + + + + + +