From 1d6b4f4345b182a3a44815b0a2e15074a2714198 Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Fri, 18 Aug 2023 15:53:29 -0400 Subject: [PATCH] [1D] Split StFlow evaluation into separate functions for each equation - pull radiation computation into a separate method - Various comment formatting/cleaning up - eliminate evalResidual and evalRightBoundary methods - moved wdot eval to update_thermo() method - removed setGas() call as it is redundant --- include/cantera/oneD/IonFlow.h | 14 +- include/cantera/oneD/StFlow.h | 39 ++- src/oneD/IonFlow.cpp | 52 ++- src/oneD/StFlow.cpp | 604 ++++++++++++++++++--------------- 4 files changed, 396 insertions(+), 313 deletions(-) diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index 5b4701bb3..76dc0cf31 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -68,11 +68,17 @@ public: protected: /*! - * This function overloads the original function. The residual function - * of electric field is added. + * This function overloads the original electric field function. + * The residual function of a non-zero electric field is added. */ - void evalResidual(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) override; + void evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) override; + /*! + * This function overloads the original species function. + * A Neumann boundary for the charged species at the left boundary is added. + */ + void evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) override; void updateTransport(double* x, size_t j0, size_t j1) override; void updateDiffFluxes(const double* x, size_t j0, size_t j1) override; //! Solving phase one: the fluxes of charged species are turned off diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 54cd21bed..2429ac4a1 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -26,7 +26,7 @@ enum offset , c_offset_V //! strain rate , c_offset_T //! temperature , c_offset_L //! (1/r)dP/dr - , c_offset_E //! electric poisson's equation + , c_offset_E //! electric field equation , c_offset_Y //! mass fractions }; @@ -232,6 +232,10 @@ public: return m_qdotRadiation[j]; } + //! Computes the radiative heat loss vector over points jmin to jmax and stores + //! the data in the qdotRadiation variable. + void computeRadiation(double* x, size_t jmin, size_t jmax); + //! Set the emissivities for the boundary values /*! * Reads the emissivities for the left and right boundary values in the @@ -300,13 +304,6 @@ public: */ void eval(size_t j, double* x, double* r, integer* mask, double rdt) override; - //! Evaluate all residual components at the right boundary. - virtual void evalRightBoundary(double* x, double* res, int* diag, double rdt); - - //! Evaluate the residual corresponding to the continuity equation at all - //! interior grid points. - virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt); - //! Index of the species on the left boundary with the largest mass fraction size_t leftExcessSpecies() const { return m_kExcessLeft; @@ -336,11 +333,30 @@ protected: //! to be updated are defined. virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax); - //! Evaluate the residual function. This function is called in eval - //! after updateProperties is called. - virtual void evalResidual(double* x, double* rsd, int* diag, + //! Evaluate the residual function for the continuity equation. + virtual void evalContinuity(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! Evaluate the residual function for the momentum equation. + virtual void evalMomentum(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax); + //! Evaluate the residual function for the energy equation. + virtual void evalEnergy(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! Evaluate the residual function for the species equation. + virtual void evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! Evaluate the residual function for the lambda equation. + virtual void evalLambda(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! Evaluate the residual function for the lambda equation. + virtual void evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + /** * Update the thermodynamic properties from point j0 to point j1 * (inclusive), based on solution x. @@ -352,6 +368,7 @@ protected: m_wtm[j] = m_thermo->meanMolecularWeight(); m_cp[j] = m_thermo->cp_mass(); m_thermo->getPartialMolarEnthalpies(&m_hk(0, j)); + m_kin->getNetProductionRates(&m_wdot(0, j)); } } diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index 066f2f449..060c68bd1 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -199,41 +199,61 @@ void IonFlow::setSolvingStage(const size_t stage) } } -void IonFlow::evalResidual(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) +void IonFlow::evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) { - StFlow::evalResidual(x, rsd, diag, rdt, jmin, jmax); + //----------------------------------------------- + // Electric field by Gauss's law + // + // dE/dz = e/eps_0 * sum(q_k*n_k) + // + // E = -dV/dz + //----------------------------------------------- + StFlow::evalElectricField(x, rsd, diag, rdt, jmin, jmax); if (m_stage != 2) { return; } for (size_t j = jmin; j <= jmax; j++) { if (j == 0) { - // enforcing the flux for charged species is difficult - // since charged species are also affected by electric - // force, so Neumann boundary condition is used. - for (size_t k : m_kCharge) { - rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1); - } rsd[index(c_offset_E, j)] = E(x,0); diag[index(c_offset_E, j)] = 0; } else if (j == m_points - 1) { rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; diag[index(c_offset_E, j)] = 0; } else { - //----------------------------------------------- - // Electric field by Gauss's law - // - // dE/dz = e/eps_0 * sum(q_k*n_k) - // - // E = -dV/dz - //----------------------------------------------- rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; diag[index(c_offset_E, j)] = 0; } } } +void IonFlow::evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //----------------------------------------------- + // Species equations + //----------------------------------------------- + StFlow::evalSpecies(x, rsd, diag, rdt, jmin, jmax); + if (m_stage != 2) { + return; + } + + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + // enforcing the flux for charged species is difficult + // since charged species are also affected by electric + // force, so Neumann boundary condition is used. + for (size_t k : m_kCharge) { + rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1); + } + } else if (j == m_points - 1) { // right boundary + } else { // interior points + } + } +} + + void IonFlow::solveElectricField(size_t j) { bool changed = false; diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index f1406f7d1..70f1810cf 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -317,7 +317,22 @@ void StFlow::eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) } updateProperties(jg, x, jmin, jmax); - evalResidual(x, rsd, diag, rdt, jmin, jmax); + + //---------------------------------------------------- + // Evaluate the residual equations at all required + // grid points. + // These residuals may be modified by a boundary object. + // The boundary object connected will modify + // these equations by subtracting its values for V, T, mdot, etc. + // As a result, these residual equations will force the solution + // variables to the values for the boundary object. + //---------------------------------------------------- + evalContinuity(x, rsd, diag, rdt, jmin, jmax); + evalMomentum(x, rsd, diag, rdt, jmin, jmax); + evalEnergy(x, rsd, diag, rdt, jmin, jmax); + evalSpecies(x, rsd, diag, rdt, jmin, jmax); + evalLambda(x, rsd, diag, rdt, jmin, jmax); + evalElectricField(x, rsd, diag, rdt, jmin, jmax); } void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax) @@ -344,190 +359,6 @@ void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax) updateDiffFluxes(x, j0, j1); } -void StFlow::evalResidual(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - //---------------------------------------------------- - // evaluate the residual equations at all required - // grid points - //---------------------------------------------------- - - // calculation of qdotRadiation (see docstring of enableRadiation) - if (m_do_radiation) { - // variable definitions for the Planck absorption coefficient and the - // radiation calculation: - double k_P_ref = 1.0*OneAtm; - - // polynomial coefficients: - const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, - 0.51382, -1.86840e-5}; - const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, - 56.310, -5.8169}; - - // calculation of the two boundary values - double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4); - double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4); - - // loop over all grid points - for (size_t j = jmin; j < jmax; j++) { - // helping variable for the calculation - double radiative_heat_loss = 0; - - // calculation of the mean Planck absorption coefficient - double k_P = 0; - // absorption coefficient for H2O - if (m_kRadiating[1] != npos) { - double k_P_H2O = 0; - for (size_t n = 0; n <= 5; n++) { - k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n); - } - k_P_H2O /= k_P_ref; - k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O; - } - // absorption coefficient for CO2 - if (m_kRadiating[0] != npos) { - double k_P_CO2 = 0; - for (size_t n = 0; n <= 5; n++) { - k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n); - } - k_P_CO2 /= k_P_ref; - k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2; - } - - // calculation of the radiative heat loss term - radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) - - boundary_Rad_left - boundary_Rad_right); - - // set the radiative heat loss vector - m_qdotRadiation[j] = radiative_heat_loss; - } - } - - for (size_t j = jmin; j <= jmax; j++) { - //---------------------------------------------- - // left boundary - //---------------------------------------------- - - if (j == 0) { - // these may be modified by a boundary object - - // Continuity. This propagates information right-to-left, since - // rho_u at point 0 is dependent on rho_u at point 1, but not on - // mdot from the inlet. - rsd[index(c_offset_U,0)] = - -(rho_u(x,1) - rho_u(x,0))/m_dz[0] - -(density(1)*V(x,1) + density(0)*V(x,0)); - - // the inlet (or other) object connected to this one will modify - // these equations by subtracting its values for V, T, and mdot. As - // a result, these residual equations will force the solution - // variables to the values for the boundary object - rsd[index(c_offset_V,0)] = V(x,0); - rsd[index(c_offset_T,0)] = T(x,0); - if (m_usesLambda) { - rsd[index(c_offset_L, 0)] = -rho_u(x, 0); - } else { - rsd[index(c_offset_L, 0)] = lambda(x, 0); - diag[index(c_offset_L, 0)] = 0; - } - - // The default boundary condition for species is zero flux. However, - // the boundary object may modify this. - double sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { - sum += Y(x,k,0); - rsd[index(c_offset_Y + k, 0)] = - -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); - } - rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum; - - // set residual of poisson's equ to zero - rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)]; - } else if (j == m_points - 1) { - evalRightBoundary(x, rsd, diag, rdt); - } else { // interior points - evalContinuity(j, x, rsd, diag, rdt); - // set residual of poisson's equ to zero - rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; - - //------------------------------------------------ - // Radial momentum equation - // - // \rho dV/dt + \rho u dV/dz + \rho V^2 - // = d(\mu dV/dz)/dz - lambda - //------------------------------------------------- - if (m_usesLambda) { - rsd[index(c_offset_V,j)] = - (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)); - diag[index(c_offset_V, j)] = 1; - } else { - rsd[index(c_offset_V, j)] = V(x, j); - diag[index(c_offset_V, j)] = 0; - } - - //------------------------------------------------- - // Species equations - // - // \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz - // = M_k\omega_k - //------------------------------------------------- - getWdot(x,j); - for (size_t k = 0; k < m_nsp; k++) { - double convec = rho_u(x,j)*dYdz(x,k,j); - double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) - / (z(j+1) - z(j-1)); - rsd[index(c_offset_Y + k, j)] - = (m_wt[k]*(wdot(k,j)) - - convec - diffus)/m_rho[j] - - rdt*(Y(x,k,j) - Y_prev(k,j)); - diag[index(c_offset_Y + k, j)] = 1; - } - - //----------------------------------------------- - // energy equation - // - // \rho c_p dT/dt + \rho c_p u dT/dz - // = d(k dT/dz)/dz - // - sum_k(\omega_k h_k_ref) - // - sum_k(J_k c_p_k / M_k) dT/dz - //----------------------------------------------- - if (m_do_energy[j]) { - - setGas(x,j); - double dtdzj = dTdz(x,j); - double sum = 0.0; - - grad_hk(x, j); - for (size_t k = 0; k < m_nsp; k++) { - double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); - sum += wdot(k,j)*m_hk(k,j); - sum += flxk * m_dhk_dz(k,j) / m_wt[k]; - } - - rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj - - divHeatFlux(x,j) - sum; - rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); - rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); - rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j])); - diag[index(c_offset_T, j)] = 1; - } else { - // residual equations if the energy equation is disabled - rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); - diag[index(c_offset_T, j)] = 0; - } - - if (m_usesLambda) { - rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1); - } else { - rsd[index(c_offset_L, j)] = lambda(x, j); - } - diag[index(c_offset_L, j)] = 0; - } - } -} - void StFlow::updateTransport(double* x, size_t j0, size_t j1) { if (m_do_multicomponent) { @@ -558,23 +389,6 @@ void StFlow::updateTransport(double* x, size_t j0, size_t j1) } } -void StFlow::show(const double* x) -{ - writelog(" Pressure: {:10.4g} Pa\n", m_press); - - Domain1D::show(x); - - if (m_do_radiation) { - writeline('-', 79, false, true); - writelog("\n z radiative heat loss"); - writeline('-', 79, false, true); - for (size_t j = 0; j < m_points; j++) { - writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]); - } - writelog("\n"); - } -} - void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1) { if (m_do_multicomponent) { @@ -617,6 +431,312 @@ void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1) } } +void StFlow::evalContinuity(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //---------------------------------------------- + // Continuity equation + // + // d(\rho u)/dz + 2\rho V = 0 + //---------------------------------------------- + // This propagates information right-to-left, since rho_u at point 0 is dependent + // on rho_u at point 1, but not on mdot from the inlet. + // The default value for continuity is zero u at the boundary. + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + rsd[index(c_offset_U,0)] = + -(rho_u(x,1) - rho_u(x,0))/m_dz[0] + -(density(1)*V(x,1) + density(0)*V(x,0)); + } else if (j == m_points - 1) { // right boundary + if (m_usesLambda) { + rsd[index(c_offset_U, j)] = rho_u(x, j); + } else { + rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1); + } + } else { // interior points + if (m_usesLambda) { + // Note that this propagates the mass flow rate information to the left + // (j+1 -> j) from the value specified at the right boundary. The + // lambda information propagates in the opposite direction. + rsd[index(c_offset_U,j)] = + -(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 (m_isFree) { + // terms involving V are zero as V=0 by definition + if (grid(j) > m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]; + } else if (grid(j) == m_zfixed) { + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) + - m_rho[0]*0.3); // why 0.3? + } + } else if (grid(j) < m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]; + } + } else { + // unstrained with fixed mass flow rate + rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1); + } + diag[index(c_offset_U, j)] = 0; // Algebraic constraint + } + } +} + +void StFlow::evalMomentum(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //------------------------------------------------ + // Radial momentum equation + // + // \rho dV/dt + \rho u dV/dz + \rho V^2 + // = d(\mu dV/dz)/dz - lambda + //------------------------------------------------- + // The default boundary condition is zero V. + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + rsd[index(c_offset_V,0)] = V(x,0); + } else if (j == m_points - 1) { // right boundary + rsd[index(c_offset_V,j)] = V(x,j); + diag[index(c_offset_V,j)] = 0; + } else { // interior points + if (m_usesLambda) { + rsd[index(c_offset_V,j)] = + (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)); + diag[index(c_offset_V, j)] = 1; + } else { + rsd[index(c_offset_V, j)] = V(x, j); + diag[index(c_offset_V, j)] = 0; + } + } + } +} + +void StFlow::evalEnergy(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //----------------------------------------------- + // Energy equation + // + // \rho c_p dT/dt + \rho c_p u dT/dz + // = d(k dT/dz)/dz + // - sum_k(\omega_k h_k_ref) + // - sum_k(J_k c_p_k / M_k) dT/dz + //----------------------------------------------- + // The default boundary condition for energy is a zero T. + + // Calculation of qdotRadiation (see docstring of enableRadiation) + if (m_do_radiation) { + computeRadiation(x, jmin, jmax); + } + + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + rsd[index(c_offset_T,0)] = T(x,0); + } else if (j == m_points - 1) { // right boundary + rsd[index(c_offset_T, j)] = T(x, j); + } else { // interior points + if (m_do_energy[j]) { + double dtdzj = dTdz(x,j); + double sum = 0.0; + + grad_hk(x, j); + for (size_t k = 0; k < m_nsp; k++) { + double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += wdot(k,j)*m_hk(k,j); + sum += flxk * m_dhk_dz(k,j) / m_wt[k]; + } + + rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj + - divHeatFlux(x,j) - sum; + rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j])); + diag[index(c_offset_T, j)] = 1; + } else { + // residual equations if the energy equation is disabled + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; + } + } + } +} + +void StFlow::evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //------------------------------------------------- + // Species equations + // + // \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz + // = M_k\omega_k + //------------------------------------------------- + // The default boundary condition for species is zero flux. + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); + } + rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum; + } else if (j == m_points - 1) { // right boundary + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); + } + rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum; + diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0; + } else { // interior points + for (size_t k = 0; k < m_nsp; k++) { + double convec = rho_u(x,j)*dYdz(x,k,j); + double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + / (z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] + = (m_wt[k]*(wdot(k,j)) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + } + } + } +} + +void StFlow::evalLambda(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //------------------------------------------------- + // Lambda equation + //------------------------------------------------- + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + if (m_usesLambda) { + rsd[index(c_offset_L, 0)] = -rho_u(x, 0); + } else { + rsd[index(c_offset_L, 0)] = lambda(x, 0); + diag[index(c_offset_L, 0)] = 0; + } + } else if (j == m_points - 1) { // right boundary + rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1); + diag[index(c_offset_L, j)] = 0; + } else { // interior points + if (m_usesLambda) { + rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1); + } else { + rsd[index(c_offset_L, j)] = lambda(x, j); + } + diag[index(c_offset_L, j)] = 0; + } + } +} + +void StFlow::evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + //----------------------------------------------- + // Electric field (by Gauss's law) + // + // dE/dz = 0 + // + // E = -dV/dz + //----------------------------------------------- + // The default value for E at the boundary is zero, and E is zero + // in the domain as well. + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { // left boundary + rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)]; + } else if (j == m_points - 1) { // right boundary + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; + } else { // interior points + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; + } + } +} + +void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax) +{ + // Variable definitions for the Planck absorption coefficient and the + // radiation calculation: + double k_P_ref = 1.0*OneAtm; + + // Polynomial coefficients: + const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, + 0.51382, -1.86840e-5}; + const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, + 56.310, -5.8169}; + + // Calculation of the two boundary values + double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4); + double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4); + + for (size_t j = jmin; j < jmax; j++) { + // calculation of the mean Planck absorption coefficient + double k_P = 0; + // Absorption coefficient for H2O + if (m_kRadiating[1] != npos) { + double k_P_H2O = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n); + } + k_P_H2O /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O; + } + // Absorption coefficient for CO2 + if (m_kRadiating[0] != npos) { + double k_P_CO2 = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n); + } + k_P_CO2 /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2; + } + + // Calculation of the radiative heat loss term + double radiative_heat_loss = 0; + radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) + - boundary_Rad_left - boundary_Rad_right); + + // set the radiative heat loss vector + m_qdotRadiation[j] = radiative_heat_loss; + } +} + +void StFlow::grad_hk(const double* x, size_t j) +{ + for(size_t k = 0; k < m_nsp; k++) { + if (u(x, j) > 0.0) { + m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1]; + } + else { + m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j]; + } + } +} + +void StFlow::show(const double* x) +{ + writelog(" Pressure: {:10.4g} Pa\n", m_press); + + Domain1D::show(x); + + if (m_do_radiation) { + writeline('-', 79, false, true); + writelog("\n z radiative heat loss"); + writeline('-', 79, false, true); + for (size_t j = 0; j < m_points; j++) { + writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]); + } + writelog("\n"); + } +} + string StFlow::componentName(size_t n) const { switch (n) { @@ -937,84 +1057,4 @@ void StFlow::fixTemperature(size_t j) } } -void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt) -{ - size_t j = m_points - 1; - - // the boundary object connected to the right of this one may modify or - // replace these equations. The default boundary conditions are zero u, V, - // and T, and zero diffusive flux for all species. - - rsd[index(c_offset_V,j)] = V(x,j); - diag[index(c_offset_V,j)] = 0; - double sum = 0.0; - // set residual of poisson's equ to zero - rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; - for (size_t k = 0; k < m_nsp; k++) { - sum += Y(x,k,j); - rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); - } - rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum; - diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0; - if (m_usesLambda) { - rsd[index(c_offset_U, j)] = rho_u(x, j); - } else { - rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1); - } - - rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1); - diag[index(c_offset_L, j)] = 0; - rsd[index(c_offset_T, j)] = T(x, j); -} - -void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt) -{ - //algebraic constraint - diag[index(c_offset_U, j)] = 0; - //---------------------------------------------- - // Continuity equation - // - // d(\rho u)/dz + 2\rho V = 0 - //---------------------------------------------- - if (m_usesLambda) { - // Note that this propagates the mass flow rate information to the left - // (j+1 -> j) from the value specified at the right boundary. The - // lambda information propagates in the opposite direction. - rsd[index(c_offset_U,j)] = - -(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 (m_isFree) { - // terms involving V are zero as V=0 by definition - if (grid(j) > m_zfixed) { - rsd[index(c_offset_U,j)] = - - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]; - } else if (grid(j) == m_zfixed) { - if (m_do_energy[j]) { - rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); - } else { - rsd[index(c_offset_U,j)] = (rho_u(x,j) - - m_rho[0]*0.3); // why 0.3? - } - } else if (grid(j) < m_zfixed) { - rsd[index(c_offset_U,j)] = - - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]; - } - } else { - // unstrained with fixed mass flow rate - rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1); - } -} - -void StFlow::grad_hk(const double* x, size_t j) -{ - for(size_t k = 0; k < m_nsp; k++) { - if (u(x, j) > 0.0) { - m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1]; - } - else { - m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j]; - } - } -} - } // namespace