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
119 auto const&
N = ip_data.N;
120 auto const& dNdx = ip_data.dNdx;
121 auto const& w = ip_data.integration_weight;
122 auto& mix_density = ip_data.mix_density;
124 auto& steam_mass_frac = ip_data.dryness;
125 auto& vapor_volume_frac = ip_data.vapor_volume_fraction;
126 auto& vapor_mass_flowrate = ip_data.vapor_mass_flow_rate;
127 auto& liquid_mass_flowrate = ip_data.liquid_mass_flow_rate;
128
129 double p_int_pt = 0.0;
130 double v_int_pt = 0.0;
131 double h_int_pt = 0.0;
132
134 h_int_pt);
135
136 double p_prev_int_pt = 0.0;
137 double v_prev_int_pt = 0.0;
138 double h_prev_int_pt = 0.0;
139
141 v_prev_int_pt, h_prev_int_pt);
142
143 double vdot_int_pt = (v_int_pt - v_prev_int_pt) / dt;
144
145
146
147 const double pi = std::numbers::pi;
148
151
152 double liquid_water_density =
153 liquid_phase
155 .template value<double>(vars, pos, t, dt);
156 double const vapour_water_density =
157 gas_phase
159 .template value<double>(vars, pos, t, dt);
160
161 double const h_sat_liq_w =
162 liquid_phase
163 .property(
165 .template value<double>(vars, pos, t, dt);
166 double const h_sat_vap_w =
167 gas_phase
168 .property(
170 .template value<double>(vars, pos, t, dt);
171
172
173
174 double const dryness = std::max(
175 0., (h_int_pt - h_sat_liq_w) / (h_sat_vap_w - h_sat_liq_w));
176 steam_mass_frac = dryness;
177
178 double const T_int_pt =
179 (dryness == 0)
180 ? liquid_phase
182 .template value<double>(vars, pos, t, dt)
183 : gas_phase
186 .template value<double>(vars, pos, t, dt);
189
190
191
192
193
194
195
196 double C_0 = 1 + 0.12 * (1 - dryness);
197
198
199
200
201
202 double const sigma_gl = 0.2358 *
203 std::pow((1 - T_int_pt / 647.096), 1.256) *
204 (1 - 0.625 * (1 - T_int_pt / 647.096));
205
206 double const u_gu =
207 1.18 * (1 - dryness) *
208 std::pow((9.81) * sigma_gl *
209 (liquid_water_density - vapour_water_density),
210 0.25) /
211 std::pow(liquid_water_density, 0.5);
212
213
215 if (dryness != 0)
216 {
217
218 using LocalJacobianMatrix =
219 Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
220 using LocalResidualVector = Eigen::Matrix<double, 1, 1>;
221 using LocalUnknownVector = Eigen::Matrix<double, 1, 1>;
222 LocalJacobianMatrix J_loc;
223
224 Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(1);
225
226 auto const update_residual = [&](LocalResidualVector& residual)
227 {
229 liquid_water_density, v_int_pt, dryness, C_0,
230 u_gu, residual);
231 };
232
233 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
234 {
236 alpha, vapour_water_density, liquid_water_density, v_int_pt,
237 dryness, C_0, u_gu,
238 jacobian);
239 };
240
241 auto const update_solution =
242 [&](LocalUnknownVector const& increment)
243 {
244
245 alpha += increment[0];
246 };
247
248 const int maximum_iterations(20);
249 const double residuum_tolerance(1.e-10);
250 const double increment_tolerance(0);
251
253 linear_solver, update_jacobian, update_residual,
254 update_solution,
255 {maximum_iterations, residuum_tolerance, increment_tolerance});
256
257 auto const success_iterations = newton_solver.
solve(J_loc);
258
259 if (!success_iterations)
260 {
262 "Attention! Steam void fraction has not been correctly "
263 "calculated!");
264 }
265 }
266
267 vapor_volume_frac =
alpha;
268
269 if (alpha == 0)
270 {
271 liquid_water_density =
272 liquid_phase
274 .template value<double>(vars, pos, t, dt);
275 }
276
277 mix_density =
278 vapour_water_density *
alpha + liquid_water_density * (1 -
alpha);
279
280 auto& mix_density_prev = ip_data.mix_density_prev;
282
283 auto const rho_dot = (mix_density - mix_density_prev) / dt;
284
285 double const liquid_water_velocity_act =
286 (
alpha == 0) ? v_int_pt
287 : (1 - dryness) * mix_density * v_int_pt /
288 (1 -
alpha) / liquid_water_density;
289 double const vapor_water_velocity_act =
291 : dryness * mix_density * v_int_pt /
292 (
alpha * vapour_water_density);
293
294 vapor_mass_flowrate = vapor_water_velocity_act * vapour_water_density *
295 pi * r_i * r_i *
alpha;
296
297 liquid_mass_flowrate = liquid_water_velocity_act *
298 liquid_water_density * pi * r_i * r_i *
300
301
302
303
304
305 double const gamma =
306 alpha * liquid_water_density * vapour_water_density * mix_density /
308 std::pow((alpha * C_0 * vapour_water_density +
309 (1 - alpha * C_0) * liquid_water_density),
310 2) *
311 std::pow((C_0 - 1) * v_int_pt + u_gu, 2);
312
313 double const miu =
315 .template value<double>(vars, pos, t, dt);
316 double const Re = mix_density * v_int_pt * 2 * r_i / miu;
317
318
319
320
321
322
323 double f = 0.0;
324 if (Re > 10 && Re <= 2400)
325 f = 16 / Re;
326 else if (Re > 2400)
327 f = std::pow(std::log(xi / 3.7 / r_i) -
328 5.02 / Re * std::log(xi / 3.7 / r_i + 13 / Re),
329 -2) /
330 16;
331
332 double Q_hx = 0;
333 double const T_r_int_pt =
N.dot(T_r);
334
336 {
337
338
339
340 const double alpha_r = k_r / rho_r / c_r;
341 const double t_d = alpha_r * t / (r_i * r_i);
342
344 if (t_d < 2.8)
345 beta = std::pow((pi * t_d), -0.5) + 0.5 -
346 0.25 * std::pow((t_d / pi), 0.5) + 0.125 * t_d;
347 else
348 beta = 2 * (1 / (std::log(4 * t_d) - 2 * 0.57722) -
349 0.57722 /
350 std::pow((std::log(4 * t_d) - 2 * 0.57722), 2));
351
352 const double P_c = 2 * pi * r_i;
353 Q_hx = P_c * k_r * (T_r_int_pt - T_int_pt) / r_i * beta;
354 }
355
356
357 double const p_r_int_pt =
N.dot(p_r);
358 double const PI_int_pt =
N.dot(
PI);
359 double Q_mx = PI_int_pt * (p_int_pt - p_r_int_pt);
360
361
362
363 double Q_mom = 0;
364 double Q_ene = 0;
365 if (Q_mx != 0)
366 {
367 Q_mom = Q_mx * v_int_pt;
368
369
372 double const h_fres =
373 liquid_phase
375 .template value<double>(vars, pos, t, dt);
376 Q_ene = Q_mx * h_fres;
377 }
378
379
380 Mvv.noalias() += w *
N.transpose() * mix_density *
N;
381
382 Mhp.noalias() += -w *
N.transpose() *
N;
383 Mhh.noalias() += w *
N.transpose() * mix_density *
N;
384
385
386 Kpv.noalias() += w * dNdx.transpose() *
N * mix_density;
387
388 Kvp.noalias() += w *
N.transpose() * dNdx;
389 Kvv.noalias() += w *
N.transpose() * rho_dot *
N;
390
391 Khh.noalias() += w *
N.transpose() * mix_density * v_int_pt * dNdx;
392
393
394 Bp.noalias() += w *
N.transpose() * rho_dot + w *
N.transpose() * Q_mx;
395
396 Bv.noalias() +=
397 w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt +
398 w * dNdx.transpose() * gamma -
399 w *
N.transpose() * f * mix_density * std::abs(v_int_pt) *
400 v_int_pt / (4 * r_i) -
401 w *
N.transpose() * Q_mom;
402
403 Bh.noalias() +=
404 -1 / 2 * w *
N.transpose() * rho_dot * v_int_pt * v_int_pt -
405 w *
N.transpose() * mix_density * v_int_pt * vdot_int_pt +
406 1 / 2 * w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt *
407 v_int_pt +
408 w *
N.transpose() * (Q_hx / pi / r_i / r_i) -
409 w *
N.transpose() * Q_ene;
410
412 {
415
416 Bv.noalias() += gravity_operator * mix_density;
417 Bh.noalias() += gravity_operator * mix_density * v_int_pt;
418 }
419 }
420
421
422
423
424
425
426
427}
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
std::optional< int > solve(JacobianMatrix &jacobian) const
void setIntegrationPoint(unsigned integration_point)
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 &)
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