32{
33 auto const local_matrix_size = local_x.size();
34
35 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
36
38 local_M_data, local_matrix_size, local_matrix_size);
40 local_K_data, local_matrix_size, local_matrix_size);
42 local_b_data, local_matrix_size);
43
44
45 auto Mvv = local_M.template block<velocity_size, velocity_size>(
47
48 auto Mhp = local_M.template block<enthalpy_size, pressure_size>(
50 auto Mhh = local_M.template block<enthalpy_size, enthalpy_size>(
52
53 auto Kpv = local_K.template block<pressure_size, velocity_size>(
55
56 auto Kvp = local_K.template block<velocity_size, pressure_size>(
58 auto Kvv = local_K.template block<velocity_size, velocity_size>(
60
61 auto Khh = local_K.template block<enthalpy_size, enthalpy_size>(
63
67
68 unsigned const n_integration_points =
70
73
75
77
78
80 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
81 auto const& gas_phase = medium.phase("Gas");
82
83
84
86
88
89
91
92
94
95 auto const r_o = r_w - t_ca;
96
97 auto const r_i = r_o - t_p;
98
99
108 auto const k_r =
111 auto const c_r =
113
114 for (unsigned ip(0); ip < n_integration_points; ip++)
115 {
117 auto const&
N = ip_data.N;
118 auto const& dNdx = ip_data.dNdx;
119 auto const& w = ip_data.integration_weight;
120 auto& mix_density = ip_data.mix_density;
122 auto& steam_mass_frac = ip_data.dryness;
123 auto& vapor_volume_frac = ip_data.vapor_volume_fraction;
124 auto& vapor_mass_flowrate = ip_data.vapor_mass_flow_rate;
125 auto& liquid_mass_flowrate = ip_data.liquid_mass_flow_rate;
126
127 pos = {
132
133 double p_int_pt = 0.0;
134 double v_int_pt = 0.0;
135 double h_int_pt = 0.0;
136
138 h_int_pt);
139
140 double p_prev_int_pt = 0.0;
141 double v_prev_int_pt = 0.0;
142 double h_prev_int_pt = 0.0;
143
145 v_prev_int_pt, h_prev_int_pt);
146
147 double vdot_int_pt = (v_int_pt - v_prev_int_pt) / dt;
148
149
150
151 const double pi = std::numbers::pi;
152
155
156 double liquid_water_density =
157 liquid_phase
159 .template value<double>(vars, pos, t, dt);
160 double const vapour_water_density =
161 gas_phase
163 .template value<double>(vars, pos, t, dt);
164
165 double const h_sat_liq_w =
166 liquid_phase
167 .property(
169 .template value<double>(vars, pos, t, dt);
170 double const h_sat_vap_w =
171 gas_phase
172 .property(
174 .template value<double>(vars, pos, t, dt);
175
176
177
178 double const dryness = std::max(
179 0., (h_int_pt - h_sat_liq_w) / (h_sat_vap_w - h_sat_liq_w));
180 steam_mass_frac = dryness;
181
182 double const T_int_pt =
183 (dryness == 0)
184 ? liquid_phase
186 .template value<double>(vars, pos, t, dt)
187 : gas_phase
190 .template value<double>(vars, pos, t, dt);
193
194
195
196
197
198
199
200 double C_0 = 1 + 0.12 * (1 - dryness);
201
202
203
204
205
206 double const sigma_gl = 0.2358 *
207 std::pow((1 - T_int_pt / 647.096), 1.256) *
208 (1 - 0.625 * (1 - T_int_pt / 647.096));
209
210 double const u_gu =
211 1.18 * (1 - dryness) *
212 std::pow((9.81) * sigma_gl *
213 (liquid_water_density - vapour_water_density),
214 0.25) /
215 std::pow(liquid_water_density, 0.5);
216
217
219 if (dryness != 0)
220 {
221
222 using LocalJacobianMatrix =
223 Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
224 using LocalResidualVector = Eigen::Matrix<double, 1, 1>;
225 using LocalUnknownVector = Eigen::Matrix<double, 1, 1>;
226 LocalJacobianMatrix J_loc;
227
228 Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(1);
229
230 auto const update_residual = [&](LocalResidualVector& residual)
231 {
233 liquid_water_density, v_int_pt, dryness, C_0,
234 u_gu, residual);
235 };
236
237 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
238 {
240 alpha, vapour_water_density, liquid_water_density, v_int_pt,
241 dryness, C_0, u_gu,
242 jacobian);
243 };
244
245 auto const update_solution =
246 [&](LocalUnknownVector const& increment)
247 {
248
249 alpha += increment[0];
250 };
251
252 const int maximum_iterations(20);
253 const double residuum_tolerance(1.e-10);
254 const double increment_tolerance(0);
255
257 linear_solver, update_jacobian, update_residual,
258 update_solution,
259 {maximum_iterations, residuum_tolerance, increment_tolerance});
260
261 auto const success_iterations = newton_solver.
solve(J_loc);
262
263 if (!success_iterations)
264 {
266 "Attention! Steam void fraction has not been correctly "
267 "calculated!");
268 }
269 }
270
271 vapor_volume_frac =
alpha;
272
273 if (alpha == 0)
274 {
275 liquid_water_density =
276 liquid_phase
278 .template value<double>(vars, pos, t, dt);
279 }
280
281 mix_density =
282 vapour_water_density *
alpha + liquid_water_density * (1 -
alpha);
283
284 auto& mix_density_prev = ip_data.mix_density_prev;
286
287 auto const rho_dot = (mix_density - mix_density_prev) / dt;
288
289 double const liquid_water_velocity_act =
290 (
alpha == 0) ? v_int_pt
291 : (1 - dryness) * mix_density * v_int_pt /
292 (1 -
alpha) / liquid_water_density;
293 double const vapor_water_velocity_act =
295 : dryness * mix_density * v_int_pt /
296 (
alpha * vapour_water_density);
297
298 vapor_mass_flowrate = vapor_water_velocity_act * vapour_water_density *
299 pi * r_i * r_i *
alpha;
300
301 liquid_mass_flowrate = liquid_water_velocity_act *
302 liquid_water_density * pi * r_i * r_i *
304
305
306
307
308
309 double const gamma =
310 alpha * liquid_water_density * vapour_water_density * mix_density /
312 std::pow((alpha * C_0 * vapour_water_density +
313 (1 - alpha * C_0) * liquid_water_density),
314 2) *
315 std::pow((C_0 - 1) * v_int_pt + u_gu, 2);
316
317 double const miu =
319 .template value<double>(vars, pos, t, dt);
320 double const Re = mix_density * v_int_pt * 2 * r_i / miu;
321
322
323
324
325
326
327 double f = 0.0;
328 if (Re > 10 && Re <= 2400)
329 f = 16 / Re;
330 else if (Re > 2400)
331 f = std::pow(std::log(xi / 3.7 / r_i) -
332 5.02 / Re * std::log(xi / 3.7 / r_i + 13 / Re),
333 -2) /
334 16;
335
336 double Q_hx = 0;
337 double const T_r_int_pt =
N.dot(T_r);
338
340 {
341
342
343
344 const double alpha_r = k_r / rho_r / c_r;
345 const double t_d = alpha_r * t / (r_i * r_i);
346
348 if (t_d < 2.8)
349 beta = std::pow((pi * t_d), -0.5) + 0.5 -
350 0.25 * std::pow((t_d / pi), 0.5) + 0.125 * t_d;
351 else
352 beta = 2 * (1 / (std::log(4 * t_d) - 2 * 0.57722) -
353 0.57722 /
354 std::pow((std::log(4 * t_d) - 2 * 0.57722), 2));
355
356 const double P_c = 2 * pi * r_i;
357 Q_hx = P_c * k_r * (T_r_int_pt - T_int_pt) / r_i * beta;
358 }
359
360
361 double const p_r_int_pt =
N.dot(p_r);
362 double const PI_int_pt =
N.dot(
PI);
363 double Q_mx = PI_int_pt * (p_int_pt - p_r_int_pt);
364
365
366
367 double Q_mom = 0;
368 double Q_ene = 0;
369 if (Q_mx != 0)
370 {
371 Q_mom = Q_mx * v_int_pt;
372
373
376 double const h_fres =
377 liquid_phase
379 .template value<double>(vars, pos, t, dt);
380 Q_ene = Q_mx * h_fres;
381 }
382
383
384 Mvv.noalias() += w *
N.transpose() * mix_density *
N;
385
386 Mhp.noalias() += -w *
N.transpose() *
N;
387 Mhh.noalias() += w *
N.transpose() * mix_density *
N;
388
389
390 Kpv.noalias() += w * dNdx.transpose() *
N * mix_density;
391
392 Kvp.noalias() += w *
N.transpose() * dNdx;
393 Kvv.noalias() += w *
N.transpose() * rho_dot *
N;
394
395 Khh.noalias() += w *
N.transpose() * mix_density * v_int_pt * dNdx;
396
397
398 Bp.noalias() += w *
N.transpose() * rho_dot + w *
N.transpose() * Q_mx;
399
400 Bv.noalias() +=
401 w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt +
402 w * dNdx.transpose() * gamma -
403 w *
N.transpose() * f * mix_density * std::abs(v_int_pt) *
404 v_int_pt / (4 * r_i) -
405 w *
N.transpose() * Q_mom;
406
407 Bh.noalias() +=
408 -1 / 2 * w *
N.transpose() * rho_dot * v_int_pt * v_int_pt -
409 w *
N.transpose() * mix_density * v_int_pt * vdot_int_pt +
410 1 / 2 * w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt *
411 v_int_pt +
412 w *
N.transpose() * (Q_hx / pi / r_i / r_i) -
413 w *
N.transpose() * Q_ene;
414
416 {
419
420 Bv.noalias() += gravity_operator * mix_density;
421 Bh.noalias() += gravity_operator * mix_density * v_int_pt;
422 }
423 }
424
425
426
427
428
429
430
431}
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
std::optional< int > solve(JacobianMatrix &jacobian) const
static const int velocity_index
static const int pressure_index
void calculateJacobian(double const alpha, double const vapor_water_density, double const liquid_water_density, double const v_mix, double const dryness, double const C_0, double const u_gu, JacobianMatrix &Jac)
static const int enthalpy_index
void calculateResidual(double const alpha, double const vapor_water_density, double const liquid_water_density, double const v_mix, double const dryness, double const C_0, double const u_gu, ResidualVector &res)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
ParameterLib::Parameter< double > const & pressure
ParameterLib::Parameter< double > const & thermal_conductivity
ParameterLib::Parameter< double > const & temperature
ParameterLib::Parameter< double > const & density
ParameterLib::Parameter< double > const & specific_heat_capacity
ParameterLib::Parameter< double > const & casing_thickness
ParameterLib::Parameter< double > const & roughness
ParameterLib::Parameter< double > const & diameter
ParameterLib::Parameter< double > const & pipe_thickness
Eigen::VectorXd const specific_body_force
ParameterLib::Parameter< double > const & productivity_index
WellboreGeometry wellbore
ReservoirProperties reservoir_properties
bool const has_heat_exchange_with_formation
Public Member Functions inherited from ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface
1.12.0