OGS
ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim > Class Template Referencefinal

Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >

Definition at line 42 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

#include <WellboreCompensateNeumannBoundaryConditionLocalAssembler.h>

Inheritance diagram for ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >:
[legend]

Public Member Functions

 WellboreCompensateNeumannBoundaryConditionLocalAssembler (MeshLib::Element const &e, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, WellboreCompensateNeumannBoundaryConditionData const &data)
void assemble (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table_boundary, double const, std::vector< GlobalVector * > const &x, int const process_id, GlobalMatrix *, GlobalVector &b, GlobalMatrix *) override
Public Member Functions inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >
 GenericNaturalBoundaryConditionLocalAssembler (MeshLib::Element const &e, bool is_axially_symmetric, NumLib::GenericIntegrationMethod const &integration_method)
Public Member Functions inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssemblerInterface
virtual ~GenericNaturalBoundaryConditionLocalAssemblerInterface ()=default

Private Types

using Base
using NodalVectorType = typename Base::NodalVectorType
using NodalMatrixType = typename Base::NodalMatrixType

Private Attributes

MeshLib::Element const & _element
WellboreCompensateNeumannBoundaryConditionData const & _data
Base::NodalVectorType _local_matrix_size

Additional Inherited Members

Protected Types inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType
using NodalVectorType = typename ShapeMatricesType::NodalVectorType
Protected Attributes inherited from ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >
NumLib::GenericIntegrationMethod const & _integration_method
std::vector< NAndWeight, Eigen::aligned_allocator< NAndWeight > > const _ns_and_weights
MeshLib::Element const & _element

Member Typedef Documentation

◆ Base

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::Base
private
Initial value:
GenericNaturalBoundaryConditionLocalAssembler<ShapeFunction, GlobalDim>
ProcessLib::GenericNaturalBoundaryConditionLocalAssembler::GenericNaturalBoundaryConditionLocalAssembler
GenericNaturalBoundaryConditionLocalAssembler(MeshLib::Element const &e, bool is_axially_symmetric, NumLib::GenericIntegrationMethod const &integration_method)
Definition GenericNaturalBoundaryConditionLocalAssembler.h:72

Definition at line 46 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

◆ NodalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::NodalMatrixType = typename Base::NodalMatrixType
private

Definition at line 49 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

◆ NodalVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::NodalVectorType = typename Base::NodalVectorType
private

Definition at line 48 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

Constructor & Destructor Documentation

◆ WellboreCompensateNeumannBoundaryConditionLocalAssembler()

template<typename ShapeFunction, int GlobalDim>
ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::WellboreCompensateNeumannBoundaryConditionLocalAssembler ( MeshLib::Element const & e,
std::size_t const local_matrix_size,
NumLib::GenericIntegrationMethod const & integration_method,
bool const is_axially_symmetric,
WellboreCompensateNeumannBoundaryConditionData const & data )
inline

The neumann_bc_term factor is directly integrated into the local element matrix.

Definition at line 54 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

60 : Base(e, is_axially_symmetric, integration_method),
61 _element(e),
62 _data(data),
63 _local_matrix_size(local_matrix_size)
64 {
65 }

References _data, _element, and _local_matrix_size.

Member Function Documentation

◆ assemble()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::assemble ( std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const & dof_table_boundary,
double const ,
std::vector< GlobalVector * > const & x,
int const process_id,
GlobalMatrix * ,
GlobalVector & b,
GlobalMatrix *  )
inlineoverridevirtual

Implements ProcessLib::GenericNaturalBoundaryConditionLocalAssemblerInterface.

Definition at line 67 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

72 {
73 NodalVectorType _local_rhs(_local_matrix_size);
74 _local_rhs.setZero();
75
76 unsigned const n_integration_points =
77 Base::_integration_method.getNumberOfPoints();
78
79 auto const indices_current_variable =
80 NumLib::getIndices(mesh_item_id, dof_table_boundary);
81 auto const indices_pressure = NumLib::getIndices(
82 mesh_item_id, *_data.dof_table_boundary_pressure);
83 auto const indices_velocity = NumLib::getIndices(
84 mesh_item_id, *_data.dof_table_boundary_velocity);
85 auto const indices_enthalpy = NumLib::getIndices(
86 mesh_item_id, *_data.dof_table_boundary_enthalpy);
87
88 std::vector<double> const local_pressure =
89 x[process_id]->get(indices_pressure);
90 std::vector<double> const local_velocity =
91 x[process_id]->get(indices_velocity);
92 std::vector<double> const local_enthalpy =
93 x[process_id]->get(indices_enthalpy);
94
95 auto const& medium = *_data.media_map.getMedium(_element.getID());
96 auto const& liquid_phase =
97 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
98 auto const& gas_phase =
99 medium.phase(MaterialPropertyLib::PhaseName::Gas);
100
101 ParameterLib::SpatialPosition pos;
102 pos.setElementID(_element.getID());
103
104 MaterialPropertyLib::VariableArray vars;
105
106 for (unsigned ip = 0; ip < n_integration_points; ip++)
107 {
108 auto const& n_and_weight = Base::_ns_and_weights[ip];
109 auto const& N = n_and_weight.N;
110 auto const& w = n_and_weight.weight;
111
112 double pressure_int_pt = 0.0;
113 double velocity_int_pt = 0.0;
114 double enthalpy_int_pt = 0.0;
115
116 NumLib::shapeFunctionInterpolate(local_pressure, N,
117 pressure_int_pt);
118 NumLib::shapeFunctionInterpolate(local_velocity, N,
119 velocity_int_pt);
120 NumLib::shapeFunctionInterpolate(local_enthalpy, N,
121 enthalpy_int_pt);
122
123 vars.liquid_phase_pressure = pressure_int_pt;
124 vars.enthalpy = enthalpy_int_pt;
125
126 double liquid_water_density =
127 liquid_phase
128 .property(
129 MaterialPropertyLib::PropertyType::saturation_density)
130 .template value<double>(vars, pos, 0, 0);
131
132 double const vapour_water_density =
133 gas_phase
134 .property(
135 MaterialPropertyLib::PropertyType::saturation_density)
136 .template value<double>(vars, pos, 0, 0);
137
138 double const h_sat_liq_w =
139 liquid_phase
140 .property(
141 MaterialPropertyLib::PropertyType::saturation_enthalpy)
142 .template value<double>(vars, pos, 0, 0);
143
144 double const h_sat_vap_w =
145 gas_phase
146 .property(
147 MaterialPropertyLib::PropertyType::saturation_enthalpy)
148 .template value<double>(vars, pos, 0, 0);
149
150 double const dryness =
151 std::max(0., (enthalpy_int_pt - h_sat_liq_w) /
152 (h_sat_vap_w - h_sat_liq_w));
153
154 double const T_int_pt =
155 (dryness == 0)
156 ? liquid_phase
157 .property(
158 MaterialPropertyLib::PropertyType::temperature)
159 .template value<double>(vars, pos, 0, 0)
160 : gas_phase
161 .property(MaterialPropertyLib::PropertyType::
162 saturation_temperature)
163 .template value<double>(vars, pos, 0, 0);
164
165 vars.temperature = T_int_pt;
166
167 // For the calculation of the void fraction of vapour,
168 // see Rohuani, Z., and E. Axelsson. "Calculation of volume void
169 // fraction in a subcooled and quality region." International
170 // Journal of Heat and Mass Transfer 17 (1970): 383-393.
171
172 // Profile parameter of drift flux
173 double const C_0 = 1 + 0.12 * (1 - dryness);
174
175 // For the surface tension calculation, see
176 // Cooper, J. R., and R. B. Dooley. "IAPWS release on surface
177 // tension of ordinary water substance." International Association
178 // for the Properties of Water and Steam (1994).
179 double const sigma_gl = 0.2358 *
180 std::pow((1 - T_int_pt / 647.096), 1.256) *
181 (1 - 0.625 * (1 - T_int_pt / 647.096));
182 // drift flux velocity
183 double const u_gu =
184 1.18 * (1 - dryness) *
185 std::pow((9.81) * sigma_gl *
186 (liquid_water_density - vapour_water_density),
187 0.25) /
188 std::pow(liquid_water_density, 0.5);
189
190 // solving void fraction of vapour using local Newton
191 // iteration.
192 double alpha = 0;
193 if (dryness != 0)
194 {
195 // Local Newton solver
196 using LocalJacobianMatrix =
197 Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
198 using LocalResidualVector = Eigen::Matrix<double, 1, 1>;
199 using LocalUnknownVector = Eigen::Matrix<double, 1, 1>;
200 LocalJacobianMatrix J_loc;
201
202 Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(1);
203
204 auto const update_residual = [&](LocalResidualVector& residual)
205 {
206 residual(0) = dryness * liquid_water_density *
207 (alpha * vapour_water_density +
208 (1 - alpha) * liquid_water_density) *
209 velocity_int_pt -
210 alpha * C_0 * dryness * liquid_water_density *
211 (alpha * vapour_water_density +
212 (1 - alpha) * liquid_water_density) *
213 velocity_int_pt -
214 alpha * C_0 * (1 - dryness) *
215 vapour_water_density *
216 (alpha * vapour_water_density +
217 (1 - alpha) * liquid_water_density) *
218 velocity_int_pt -
219 alpha * vapour_water_density *
220 liquid_water_density * u_gu;
221 };
222
223 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
224 {
225 jacobian(0) =
226 dryness * liquid_water_density * velocity_int_pt *
227 (vapour_water_density - liquid_water_density) -
228 (C_0 * dryness * liquid_water_density +
229 C_0 * (1 - dryness) * vapour_water_density) *
230 (2 * alpha * vapour_water_density +
231 (1 - 2 * alpha) * liquid_water_density) *
232 velocity_int_pt -
233 vapour_water_density * liquid_water_density * u_gu;
234 };
235
236 auto const update_solution =
237 [&](LocalUnknownVector const& increment)
238 {
239 // increment solution vectors
240 alpha += increment[0];
241 };
242
243 const int maximum_iterations(20);
244 const double residuum_tolerance(1.e-10);
245 const double increment_tolerance(0);
246
247 auto newton_solver = NumLib::NewtonRaphson(
248 linear_solver, update_jacobian, update_residual,
249 update_solution,
250 {maximum_iterations, residuum_tolerance,
251 increment_tolerance});
252
253 auto const success_iterations = newton_solver.solve(J_loc);
254
255 if (!success_iterations)
256 {
257 WARN(
258 "Attention! Steam void fraction has not been correctly "
259 "calculated!");
260 }
261 }
262
263 if (alpha == 0)
264 {
265 liquid_water_density =
266 liquid_phase
267 .property(MaterialPropertyLib::PropertyType::density)
268 .template value<double>(vars, pos, 0, 0);
269 }
270
271 double const mix_density = vapour_water_density * alpha +
272 liquid_water_density * (1 - alpha);
273
274 // slip parameter between two phases
275 double const gamma =
276 alpha * liquid_water_density * vapour_water_density *
277 mix_density / (1 - alpha) /
278 std::pow((alpha * C_0 * vapour_water_density +
279 (1 - alpha * C_0) * liquid_water_density),
280 2) *
281 std::pow((C_0 - 1) * velocity_int_pt + u_gu, 2);
282
283 double const neumann_ip_values =
284 _data.coefficients.pressure * mix_density * velocity_int_pt +
285 _data.coefficients.velocity *
286 (mix_density * velocity_int_pt * velocity_int_pt + gamma) +
287 _data.coefficients.enthalpy * mix_density * velocity_int_pt *
288 velocity_int_pt * velocity_int_pt * 0.5;
289 _local_rhs.noalias() += N.transpose() * neumann_ip_values * w;
290 }
291
292 b.add(indices_current_variable, _local_rhs);
293 }
WARN
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:34
ProcessLib::GenericNaturalBoundaryConditionLocalAssembler::_ns_and_weights
std::vector< NAndWeight, Eigen::aligned_allocator< NAndWeight > > const _ns_and_weights
Definition GenericNaturalBoundaryConditionLocalAssembler.h:85
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:20
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:105

References _data, _element, ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::_integration_method, _local_matrix_size, ProcessLib::GenericNaturalBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::_ns_and_weights, MathLib::EigenVector::add(), MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::enthalpy, MaterialPropertyLib::Gas, NumLib::getIndices(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::saturation_density, MaterialPropertyLib::saturation_enthalpy, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), NumLib::NewtonRaphson< LinearSolver, JacobianMatrixUpdate, ResidualUpdate, SolutionUpdate >::solve(), MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, and WARN().

Member Data Documentation

◆ _data

template<typename ShapeFunction, int GlobalDim>
WellboreCompensateNeumannBoundaryConditionData const& ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::_data
private

Definition at line 297 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

Referenced by WellboreCompensateNeumannBoundaryConditionLocalAssembler(), and assemble().

◆ _element

template<typename ShapeFunction, int GlobalDim>
MeshLib::Element const& ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::_element
private

Definition at line 296 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

Referenced by WellboreCompensateNeumannBoundaryConditionLocalAssembler(), and assemble().

◆ _local_matrix_size

template<typename ShapeFunction, int GlobalDim>
Base::NodalVectorType ProcessLib::WellboreCompensateNeumannBoundaryConditionLocalAssembler< ShapeFunction, GlobalDim >::_local_matrix_size
private

Definition at line 298 of file WellboreCompensateNeumannBoundaryConditionLocalAssembler.h.

Referenced by WellboreCompensateNeumannBoundaryConditionLocalAssembler(), and assemble().

The documentation for this class was generated from the following file: