27 template <
typename Traits>
30 const unsigned dimension)
31 : _d{ap, num_int_pts, dimension}
35 template <
typename Traits>
37 const unsigned int_pt)
41 _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
44 const double M_pp = _d.ap.poro / _d.p * _d.rho_GR;
45 const double M_pT = -_d.ap.poro / _d.T * _d.rho_GR;
46 const double M_px = (_d.ap.M_react - _d.ap.M_inert) * _d.p / (R * _d.T) *
49 const double M_Tp = -_d.ap.poro;
51 _d.ap.poro * _d.rho_GR * _d.ap.cpG
53 (1.0 - _d.ap.poro) * _d.solid_density[int_pt] *
55 const double M_Tx = 0.0;
57 const double M_xp = 0.0;
58 const double M_xT = 0.0;
59 const double M_xx = _d.ap.poro * _d.rho_GR;
62 M << M_pp, M_pT, M_px, M_Tp, M_TT, M_Tx, M_xp, M_xT, M_xx;
67 template <
typename Traits>
68 typename Traits::LaplaceMatrix
72 const double eta_GR =
fluid_viscosity(_d.p, _d.T, _d.vapour_mass_fraction);
74 const double lambda_F =
76 const double lambda_S = _d.ap.solid_heat_cond;
78 using Mat =
typename Traits::MatrixDimDim;
80 typename Traits::LaplaceMatrix L =
85 Traits::blockDimDim(L, 0, 0, dim, dim) =
86 Traits::blockDimDim(_d.ap.solid_perm_tensor, 0, 0, dim, dim) *
91 Traits::blockDimDim(L, dim, dim, dim, dim) =
92 Mat::Identity(dim, dim) *
93 (_d.ap.poro * lambda_F + (1.0 - _d.ap.poro) * lambda_S);
96 Traits::blockDimDim(L, 2 * dim, 2 * dim, dim, dim) =
97 Mat::Identity(dim, dim) * (_d.ap.tortuosity * _d.ap.poro * _d.rho_GR *
98 _d.ap.diffusion_coefficient_component);
103 template <
typename Traits>
107 const double A_pp = 0.0;
108 const double A_pT = 0.0;
110 const double A_px = 0.0;
112 const double A_Tp = 0.0;
114 const double A_TT = _d.rho_GR * _d.ap.cpG;
115 const double A_Tx = 0.0;
117 const double A_xp = 0.0;
118 const double A_xT = 0.0;
119 const double A_xx = _d.rho_GR;
122 A << A_pp, A_pT, A_px, A_Tp, A_TT, A_Tx, A_xp, A_xT, A_xx;
127 template <
typename Traits>
131 const double C_pp = 0.0;
132 const double C_pT = 0.0;
134 const double C_px = 0.0;
136 const double C_Tp = 0.0;
138 const double C_TT = 0.0;
139 const double C_Tx = 0.0;
141 const double C_xp = 0.0;
142 const double C_xT = 0.0;
143 const double C_xx = (_d.ap.poro - 1.0) * _d.qR;
146 C << C_pp, C_pT, C_px, C_Tp, C_TT, C_Tx, C_xp, C_xT, C_xx;
151 template <
typename Traits>
153 const unsigned int_pt)
155 const double reaction_enthalpy =
156 _d.ap.react_sys->getEnthalpy(_d.p_V, _d.T, _d.ap.M_react);
159 (_d.ap.poro - 1.0) * _d.qR;
162 _d.rho_GR * _d.ap.poro * _d.ap.fluid_specific_heat_source +
163 (1.0 - _d.ap.poro) * _d.qR * reaction_enthalpy +
164 _d.solid_density[int_pt] * (1.0 - _d.ap.poro) *
165 _d.ap.solid_specific_heat_source;
169 (_d.ap.poro - 1.0) * _d.qR;
172 rhs << rhs_p, rhs_T, rhs_x;
177 template <
typename Traits>
180 auto const& rate = _d.reaction_adaptor->
initReaction(int_pt);
181 _d.qR = rate.reaction_rate;
182 _d.reaction_rate[int_pt] = rate.reaction_rate;
183 _d.solid_density[int_pt] = rate.solid_density;
186 template <
typename Traits>
188 const unsigned int_pt,
189 const std::vector<double>& localX,
190 typename Traits::ShapeMatrices
const& sm)
194 _d.p = _d.T = _d.vapour_mass_fraction = _d.p_V = _d.rho_GR = _d.qR =
195 std::numeric_limits<double>::quiet_NaN();
199 _d.vapour_mass_fraction);
203 _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
205 initReaction(int_pt);
209 assert(0.0 <= _d.vapour_mass_fraction && _d.vapour_mass_fraction <= 1.0);
211 _d.rho_GR =
fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
214 template <
typename Traits>
216 unsigned integration_point,
217 std::vector<double>
const& localX,
218 typename Traits::ShapeMatrices
const& sm,
220 Eigen::Map<typename Traits::LocalMatrix>& local_M,
221 Eigen::Map<typename Traits::LocalMatrix>& local_K,
222 Eigen::Map<typename Traits::LocalVector>& local_b)
224 preEachAssembleIntegrationPoint(integration_point, localX, sm);
227 static_cast<unsigned>(sm.dNdx.cols());
229 static_cast<unsigned>(sm.dNdx.rows());
233 auto const laplaceCoeffMat = getLaplaceCoeffMatrix(integration_point, D);
234 assert(laplaceCoeffMat.cols() == D *
NODAL_DOF);
235 auto const massCoeffMat = getMassCoeffMatrix(integration_point);
236 auto const advCoeffMat = getAdvectionCoeffMatrix(integration_point);
237 auto const contentCoeffMat = getContentCoeffMatrix(integration_point);
240 assert((
unsigned)sm.dNdx.rows() == _d.velocity.size() &&
241 (
unsigned)sm.dNdx.cols() == _d.velocity[0].size());
243 auto const velocity =
244 (Traits::blockDimDim(laplaceCoeffMat, 0, 0, D, D) *
245 (sm.dNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
250 assert(velocity.size() == D);
252 for (
unsigned d = 0; d < D; ++d)
254 _d.velocity[d][integration_point] = velocity[d];
257 auto const detJ_w_im_NT =
258 (sm.detJ * weight * sm.integralMeasure * sm.N.transpose()).eval();
259 auto const detJ_w_im_NT_N = (detJ_w_im_NT * sm.N).eval();
260 assert(detJ_w_im_NT_N.rows() == N && detJ_w_im_NT_N.cols() == N);
262 auto const detJ_w_im_NT_vT_dNdx =
263 (detJ_w_im_NT * velocity.transpose() * sm.dNdx).eval();
264 assert(detJ_w_im_NT_vT_dNdx.rows() == N && detJ_w_im_NT_vT_dNdx.cols() == N);
270 Traits::blockShpShp(local_K, N *
r, N *
c, N, N).noalias() +=
271 sm.detJ * weight * sm.integralMeasure * sm.dNdx.transpose() *
272 Traits::blockDimDim(laplaceCoeffMat, D *
r, D *
c, D, D) *
274 + detJ_w_im_NT_N * contentCoeffMat(
r,
c) +
275 detJ_w_im_NT_vT_dNdx * advCoeffMat(
r,
c);
276 Traits::blockShpShp(local_M, N *
r, N *
c, N, N).noalias() +=
277 detJ_w_im_NT_N * massCoeffMat(
r,
c);
281 auto const rhsCoeffVector = getRHSCoeffVector(integration_point);
285 Traits::blockShp(local_b, N *
r, N).noalias() +=
286 rhsCoeffVector(
r) * sm.N.transpose() * sm.detJ * weight *
291 template <
typename Traits>
294 if (_d.ap.iteration_in_current_timestep == 1)
296 if (_d.ap.number_of_try_of_iteration ==
299 _d.solid_density_prev_ts = _d.solid_density;
300 _d.reaction_rate_prev_ts = _d.reaction_rate;
302 _d.reaction_adaptor->preZerothTryAssemble();
306 _d.solid_density = _d.solid_density_prev_ts;
static double getMolarFraction(double xm, double M_this, double M_other)
static double dMolarFraction(double xm, double M_this, double M_other)
void initReaction(const unsigned int_pt)
TESLocalAssemblerInner(AssemblyParams const &ap, const unsigned num_int_pts, const unsigned dimension)
constexpr double IdealGasConstant
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
double fluid_heat_conductivity(const double p, const double T, const double x)
double fluid_viscosity(const double p, const double T, const double x)
double fluid_density(const double p, const double T, const double x)