Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Q2PBorehole.h"
11 : #include "RotationMatrix.h"
12 :
13 : registerMooseObject("RichardsApp", Q2PBorehole);
14 :
15 : InputParameters
16 16 : Q2PBorehole::validParams()
17 : {
18 16 : InputParameters params = PeacemanBorehole::validParams();
19 32 : params.addRequiredParam<UserObjectName>(
20 : "fluid_density",
21 : "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
22 32 : params.addRequiredParam<UserObjectName>(
23 : "fluid_relperm",
24 : "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
25 : "fluid relative permeability as a function of the saturation Variable.");
26 32 : params.addRequiredCoupledVar("other_var",
27 : "The other variable in the 2-phase system. If "
28 : "Variable=porepressure, the other_var=saturation, and "
29 : "vice-versa.");
30 32 : params.addRequiredParam<bool>("var_is_porepressure",
31 : "This flag is needed to correctly calculate the Jacobian entries. "
32 : "If set to true, this Sink will extract fluid from the phase with "
33 : "porepressure as its Variable (usually the liquid phase). If set "
34 : "to false, this Sink will extract fluid from the phase with "
35 : "saturation as its variable (usually the gas phase)");
36 32 : params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
37 16 : params.addClassDescription("Approximates a borehole in the mesh with given bottomhole pressure, "
38 : "and radii using a number of point sinks whose positions are read "
39 : "from a file. This DiracKernel is for use by Q2P models");
40 16 : return params;
41 0 : }
42 :
43 8 : Q2PBorehole::Q2PBorehole(const InputParameters & parameters)
44 : : PeacemanBorehole(parameters),
45 8 : _density(getUserObject<RichardsDensity>("fluid_density")),
46 8 : _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
47 8 : _other_var_nodal(coupledDofValues("other_var")),
48 8 : _other_var_num(coupled("other_var")),
49 16 : _var_is_pp(getParam<bool>("var_is_porepressure")),
50 16 : _viscosity(getParam<Real>("fluid_viscosity")),
51 16 : _permeability(getMaterialProperty<RealTensorValue>("permeability")),
52 8 : _num_nodes(0),
53 8 : _pp(0),
54 8 : _sat(0),
55 8 : _mobility(0),
56 8 : _dmobility_dp(0),
57 16 : _dmobility_ds(0)
58 : {
59 8 : }
60 :
61 : void
62 746 : Q2PBorehole::prepareNodalValues()
63 : {
64 746 : _num_nodes = _other_var_nodal.size();
65 :
66 : // set _pp and _sat variables
67 746 : _pp.resize(_num_nodes);
68 746 : _sat.resize(_num_nodes);
69 746 : if (_var_is_pp)
70 : {
71 3357 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
72 : {
73 2984 : _pp[nodenum] = _var.dofValues()[nodenum];
74 2984 : _sat[nodenum] = _other_var_nodal[nodenum];
75 : }
76 : }
77 : else
78 : {
79 3357 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
80 : {
81 2984 : _pp[nodenum] = _other_var_nodal[nodenum];
82 2984 : _sat[nodenum] = _var.dofValues()[nodenum];
83 : }
84 : }
85 :
86 : Real density;
87 : Real ddensity_dp;
88 : Real relperm;
89 : Real drelperm_ds;
90 746 : _mobility.resize(_num_nodes);
91 746 : _dmobility_dp.resize(_num_nodes);
92 746 : _dmobility_ds.resize(_num_nodes);
93 6714 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
94 : {
95 5968 : density = _density.density(_pp[nodenum]);
96 5968 : ddensity_dp = _density.ddensity(_pp[nodenum]);
97 5968 : relperm = _relperm.relperm(_sat[nodenum]);
98 5968 : drelperm_ds = _relperm.drelperm(_sat[nodenum]);
99 5968 : _mobility[nodenum] = density * relperm / _viscosity;
100 5968 : _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
101 5968 : _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
102 : }
103 746 : }
104 :
105 : void
106 452 : Q2PBorehole::computeResidual()
107 : {
108 452 : prepareNodalValues();
109 452 : DiracKernel::computeResidual();
110 452 : }
111 :
112 : Real
113 7232 : Q2PBorehole::computeQpResidual()
114 : {
115 7232 : const Real character = _character.value(_t, _q_point[_qp]);
116 7232 : if (character == 0.0)
117 : return 0.0;
118 :
119 : const Real bh_pressure =
120 7232 : _p_bot +
121 : _unit_weight *
122 7232 : (_q_point[_qp] -
123 7232 : _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
124 :
125 : // Get the ID we initially assigned to this point
126 7232 : const unsigned current_dirac_ptid = currentPointCachedID();
127 :
128 : // If getting the ID failed, fall back to the old bodge!
129 : // if (current_dirac_ptid == libMesh::invalid_uint)
130 : // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
131 :
132 : Real outflow(0.0); // this is the flow rate from porespace out of the system
133 :
134 : Real wc(0.0);
135 7232 : if (current_dirac_ptid > 0)
136 : // contribution from half-segment "behind" this point (must have >1 point for
137 : // current_dirac_ptid>0)
138 : {
139 3616 : wc = wellConstant(_permeability[0],
140 : _rot_matrix[current_dirac_ptid - 1],
141 3616 : _half_seg_len[current_dirac_ptid - 1],
142 3616 : _current_elem,
143 3616 : _rs[current_dirac_ptid]);
144 3616 : if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
145 : // injection, so outflow<0 || // production, so outflow>0
146 3616 : outflow +=
147 3616 : _test[_i][_qp] * std::abs(character) * wc * _mobility[_i] * (_pp[_i] - bh_pressure);
148 : }
149 :
150 7232 : if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
151 : // contribution from half-segment "ahead of" this point, or we only have one point
152 : {
153 3616 : wc = wellConstant(_permeability[0],
154 : _rot_matrix[current_dirac_ptid],
155 : _half_seg_len[current_dirac_ptid],
156 3616 : _current_elem,
157 3616 : _rs[current_dirac_ptid]);
158 3616 : if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
159 : // injection, so outflow<0 || // production, so outflow>0
160 3616 : outflow +=
161 3616 : _test[_i][_qp] * std::abs(character) * wc * _mobility[_i] * (_pp[_i] - bh_pressure);
162 : }
163 :
164 7232 : _total_outflow_mass.add(
165 7232 : outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
166 7232 : return outflow;
167 : }
168 :
169 : void
170 294 : Q2PBorehole::computeJacobian()
171 : {
172 294 : prepareNodalValues();
173 294 : DiracKernel::computeJacobian();
174 294 : }
175 :
176 : Real
177 37632 : Q2PBorehole::computeQpJacobian()
178 : {
179 37632 : return jac(_var.number());
180 : }
181 :
182 : Real
183 37632 : Q2PBorehole::computeQpOffDiagJacobian(unsigned int jvar)
184 : {
185 37632 : if (jvar == _other_var_num || jvar == _var.number())
186 37632 : return jac(jvar);
187 : return 0.0;
188 : }
189 :
190 : Real
191 75264 : Q2PBorehole::jac(unsigned int jvar)
192 : {
193 75264 : if (_i != _j)
194 : return 0.0;
195 :
196 9408 : const Real character = _character.value(_t, _q_point[_qp]);
197 9408 : if (character == 0.0)
198 : return 0.0;
199 :
200 : const Real bh_pressure =
201 9408 : _p_bot +
202 : _unit_weight *
203 9408 : (_q_point[_qp] -
204 9408 : _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
205 :
206 : const Real phi = 1;
207 :
208 : // Get the ID we initially assigned to this point
209 9408 : const unsigned current_dirac_ptid = currentPointCachedID();
210 :
211 : // If getting the ID failed, fall back to the old bodge!
212 : // if (current_dirac_ptid == libMesh::invalid_uint)
213 : // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
214 :
215 : Real outflowp(0.0);
216 :
217 : const bool deriv_wrt_pp =
218 9408 : (_var_is_pp && (jvar == _var.number())) || (!_var_is_pp && (jvar == _other_var_num));
219 :
220 : Real wc(0.0);
221 9408 : if (current_dirac_ptid > 0)
222 : // contribution from half-segment "behind" this point
223 : {
224 4704 : wc = wellConstant(_permeability[0],
225 : _rot_matrix[current_dirac_ptid - 1],
226 4704 : _half_seg_len[current_dirac_ptid - 1],
227 4704 : _current_elem,
228 4704 : _rs[current_dirac_ptid]);
229 4704 : if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
230 : {
231 : // injection, so outflow<0 || // production, so outflow>0
232 4704 : if (deriv_wrt_pp)
233 2352 : outflowp += std::abs(character) * wc *
234 2352 : (_mobility[_i] * phi + _dmobility_dp[_i] * phi * (_pp[_i] - bh_pressure));
235 : else
236 2352 : outflowp += std::abs(character) * wc * _dmobility_ds[_i] * phi * (_pp[_i] - bh_pressure);
237 : }
238 : }
239 :
240 9408 : if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
241 : // contribution from half-segment "ahead of" this point
242 : {
243 4704 : wc = wellConstant(_permeability[0],
244 : _rot_matrix[current_dirac_ptid],
245 : _half_seg_len[current_dirac_ptid],
246 4704 : _current_elem,
247 : _rs[current_dirac_ptid]);
248 4704 : if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
249 : {
250 : // injection, so outflow<0 || // production, so outflow>0
251 4704 : if (deriv_wrt_pp)
252 2352 : outflowp += std::abs(character) * wc *
253 2352 : (_mobility[_i] * phi + _dmobility_dp[_i] * phi * (_pp[_i] - bh_pressure));
254 : else
255 2352 : outflowp += std::abs(character) * wc * _dmobility_ds[_i] * phi * (_pp[_i] - bh_pressure);
256 : }
257 : }
258 :
259 9408 : return _test[_i][_qp] * outflowp;
260 : }
|
