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