Flow123d  master-f44eb46
heat_model.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file heat_model.cc
15  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
16  * @author Jan Stebel
17  */
18 
19 #include "input/input_type.hh"
20 #include "mesh/mesh.h"
21 #include "mesh/accessors.hh"
22 //#include "transport/transport_operator_splitting.hh"
23 #include "heat_model.hh"
24 #include "tools/unit_si.hh"
25 #include "coupling/balance.hh"
26 #include "fields/field_model.hh"
27 #include "fields/field_constant.hh"
28 
29 
30 
31 using namespace std;
32 using namespace Input::Type;
33 
34 
35 /*******************************************************************************
36  * Functors of FieldModels
37  */
38 using Sclr = double;
39 using Vect = arma::vec3;
40 using Tens = arma::mat33;
41 
42 // Functor computing velocity norm
44  inline Sclr operator() (Vect vel) {
45  return arma::norm(vel, 2);
46  }
47 };
48 
49 /**
50  * Functor computing mass matrix coefficients:
51  * cross_section * (porosity*fluid_density*fluid_heat_capacity + (1.-porosity)*solid_density*solid_heat_capacity)
52  */
54  inline Sclr operator() (Sclr csec, Sclr por, Sclr f_rho, Sclr f_c, Sclr s_rho, Sclr s_c) {
55  return csec * (por*f_rho*f_c + (1.-por)*s_rho*s_c);
56  }
57 };
58 
59 /**
60  * Functor computing sources density output:
61  * cross_section * (porosity*fluid_thermal_source + (1-porosity)*solid_thermal_source)
62  */
64  inline Sclr operator() (Sclr csec, Sclr por, Sclr f_source, Sclr s_source) {
65  return csec * (por * f_source + (1.-por) * s_source);
66  }
67 };
68 
69 /**
70  * Functor computing sources sigma output:
71  * cross_section * (porosity*fluid_density*fluid_heat_capacity*fluid_heat_exchange_rate + (1-porosity)*solid_density*solid_heat_capacity*solid_thermal_source)
72  */
74  inline Sclr operator() (Sclr csec, Sclr por, Sclr f_rho, Sclr f_cap, Sclr f_sigma, Sclr s_rho, Sclr s_cap, Sclr s_sigma) {
75  return csec * (por * f_rho * f_cap * f_sigma + (1.-por) * s_rho * s_cap * s_sigma);
76  }
77 };
78 
79 /**
80  * Functor computing sources concentration output for positive heat_sources_sigma (otherwise return 0):
81  * cross_section * (porosity*fluid_density*fluid_heat_capacity*fluid_heat_exchange_rate*fluid_ref_temperature + (1-porosity)*solid_density*solid_heat_capacity*solid_heat_exchange_rate*solid_ref_temperature)
82  */
84  inline Sclr operator() (Sclr csec, Sclr por, Sclr f_rho, Sclr f_cap, Sclr f_sigma, Sclr f_temp, Sclr s_rho, Sclr s_cap, Sclr s_sigma, Sclr s_temp, Sclr sigma) {
85  if (fabs(sigma) > numeric_limits<double>::epsilon())
86  return csec * (por * f_rho * f_cap * f_sigma * f_temp + (1.-por) * s_rho * s_cap * s_sigma * s_temp);
87  else {
88  return 0;
89  }
90  }
91 };
92 
93 /**
94  * Functor computing advection coefficient
95  * velocity * fluid_density * fluid_heat_capacity
96  */
98  inline Vect operator() (Sclr f_rho, Sclr f_cap, Vect velocity) {
99  return velocity * f_rho * f_cap;
100  }
101 };
102 
103 /**
104  * Functor computing diffusion coefficient (see notes in function)
105  */
107  inline Tens operator() (Vect velocity, Sclr v_norm, Sclr f_rho, Sclr disp_l, Sclr disp_t, Sclr f_cond, Sclr s_cond, Sclr c_sec, Sclr por) {
108  // result
109  Tens dif_coef;
110 
111  // dispersive part of thermal diffusion
112  // Note that the velocity vector is in fact the Darcian flux,
113  // so to obtain |v| we have to divide vnorm by porosity and cross_section.
114  if ( fabs(v_norm) > 0 )
115  for (int i=0; i<3; i++)
116  for (int j=0; j<3; j++)
117  dif_coef(i,j) = ( (velocity(i)/v_norm) * (velocity(j)/v_norm) * (disp_l-disp_t) + disp_t*(i==j?1:0))*v_norm*f_rho*f_cond;
118  else
119  dif_coef.zeros();
120 
121  // conductive part of thermal diffusion
122  dif_coef += c_sec * (por*f_cond + (1.-por)*s_cond) * arma::eye(3,3);
123  return dif_coef;
124  }
125 };
126 
127 
128 
129 
130 
131 
132 
133 
135 {
136  *this+=bc_type
137  .name("bc_type")
138  .description(
139  "Type of boundary condition.")
140  .units( UnitSI::dimensionless() )
141  .input_default("\"inflow\"")
142  .input_selection( get_bc_type_selection() )
144  *this+=bc_dirichlet_value
145  .name("bc_temperature")
146  .description("Boundary value of temperature.")
147  .units( UnitSI().K() )
148  .input_default("0.0")
149  .flags_add(in_rhs);
150  *this+=bc_flux
151  .disable_where(bc_type, { bc_dirichlet, bc_inflow })
152  .name("bc_flux")
153  .description("Flux in Neumann boundary condition.")
154  .units( UnitSI().kg().m().s(-1).md() )
155  .input_default("0.0")
156  .flags_add(FieldFlag::in_rhs);
157  *this+=bc_robin_sigma
158  .disable_where(bc_type, { bc_dirichlet, bc_inflow })
159  .name("bc_robin_sigma")
160  .description("Conductivity coefficient in Robin boundary condition.")
161  .units( UnitSI().m(4).s(-1).md() )
162  .input_default("0.0")
164 
165  *this+=init_condition
166  .name("init_temperature")
167  .description("Initial temperature.")
168  .units( UnitSI().K() )
169  .input_default("0.0");
170 
171  *this+=porosity
172  .name("porosity")
173  .description("Porosity.")
174  .units( UnitSI::dimensionless() )
175  .input_default("1.0")
176  .flags_add(in_main_matrix & in_time_term)
177  .set_limits(0.0);
178 
179  *this+=water_content
180  .name("water_content")
181  .units( UnitSI::dimensionless() )
182  .input_default("1.0")
183  .flags_add(input_copy & in_main_matrix & in_time_term);
184 
185  *this += flow_flux.name("flow_flux")
186  .flags( FieldFlag::input_copy )
187  .flags_add(in_time_term & in_main_matrix & in_rhs);
188 
189  *this+=fluid_density
190  .name("fluid_density")
191  .description("Density of fluid.")
192  .units( UnitSI().kg().m(-3) )
193  .input_default("1000")
194  .flags_add(in_main_matrix & in_time_term);
195 
196  *this+=fluid_heat_capacity
197  .name("fluid_heat_capacity")
198  .description("Heat capacity of fluid.")
199  .units( UnitSI::J() * UnitSI().kg(-1).K(-1) )
200  .flags_add(in_main_matrix & in_time_term);
201 
202  *this+=fluid_heat_conductivity
203  .name("fluid_heat_conductivity")
204  .description("Heat conductivity of fluid.")
205  .units( UnitSI::W() * UnitSI().m(-1).K(-1) )
206  .flags_add(in_main_matrix)
207  .set_limits(0.0);
208 
209 
210  *this+=solid_density
211  .name("solid_density")
212  .description("Density of solid (rock).")
213  .units( UnitSI().kg().m(-3) )
214  .input_default("0.0")
215  .flags_add(in_time_term);
216 
217  *this+=solid_heat_capacity
218  .name("solid_heat_capacity")
219  .description("Heat capacity of solid (rock).")
220  .units( UnitSI::J() * UnitSI().kg(-1).K(-1) )
221  .input_default("0.0")
222  .flags_add(in_time_term);
223 
224  *this+=solid_heat_conductivity
225  .name("solid_heat_conductivity")
226  .description("Heat conductivity of solid (rock).")
227  .units( UnitSI::W() * UnitSI().m(-1).K(-1) )
228  .input_default("0.0")
229  .flags_add(in_main_matrix)
230  .set_limits(0.0);
231 
232  *this+=disp_l
233  .name("disp_l")
234  .description("Longitudinal heat dispersivity in fluid.")
235  .units( UnitSI().m() )
236  .input_default("0.0")
237  .flags_add(in_main_matrix);
238 
239  *this+=disp_t
240  .name("disp_t")
241  .description("Transverse heat dispersivity in fluid.")
242  .units( UnitSI().m() )
243  .input_default("0.0")
244  .flags_add(in_main_matrix);
245 
246  *this+=fluid_thermal_source
247  .name("fluid_thermal_source")
248  .description("Density of thermal source in fluid.")
249  .units( UnitSI::W() * UnitSI().m(-3) )
250  .input_default("0.0")
251  .flags_add(in_rhs);
252 
253  *this+=solid_thermal_source
254  .name("solid_thermal_source")
255  .description("Density of thermal source in solid.")
256  .units( UnitSI::W() * UnitSI().m(-3) )
257  .input_default("0.0")
258  .flags_add(in_rhs);
259 
260  *this+=fluid_heat_exchange_rate
261  .name("fluid_heat_exchange_rate")
262  .description("Heat exchange rate of source in fluid.")
263  .units( UnitSI().s(-1) )
264  .input_default("0.0")
265  .flags_add(in_rhs);
266 
267  *this+=solid_heat_exchange_rate
268  .name("solid_heat_exchange_rate")
269  .description("Heat exchange rate of source in solid.")
270  .units( UnitSI().s(-1) )
271  .input_default("0.0")
272  .flags_add(in_rhs);
273 
274  *this+=fluid_ref_temperature
275  .name("fluid_ref_temperature")
276  .description("Reference temperature of source in fluid.")
277  .units( UnitSI().K() )
278  .input_default("0.0")
279  .flags_add(in_rhs);
280 
281  *this+=solid_ref_temperature
282  .name("solid_ref_temperature")
283  .description("Reference temperature in solid.")
284  .units( UnitSI().K() )
285  .input_default("0.0")
286  .flags_add(in_rhs);
287 
288  *this+=cross_section
289  .name("cross_section")
290  .units( UnitSI().m(3).md() )
291  .flags(input_copy & in_time_term & in_main_matrix);
292 
293  *this+=output_field
294  .name("temperature")
295  .description("Temperature solution.")
296  .units( UnitSI().K() )
297  .flags(equation_result);
298 
299 
300  // initiaization of FieldModels
301  *this += v_norm.name("v_norm")
302  .description("Velocity norm field.")
303  .input_default("0.0")
304  .units( UnitSI().m().s(-1) );
305 
306  *this += mass_matrix_coef.name("mass_matrix_coef")
307  .description("Matrix coefficients computed by model in mass assemblation.")
308  .input_default("0.0")
309  .units( UnitSI().m(3).md() );
310 
311  *this += retardation_coef.name("retardation_coef")
312  .description("Retardation coefficients computed by model in mass assemblation.")
313  .input_default("0.0")
314  .units( UnitSI().m(3).md() );
315 
316  *this += sources_conc_out.name("sources_conc_out")
317  .description("Concentration sources output.")
318  .input_default("0.0")
319  .units( UnitSI().kg().m(-3) );
320 
321  *this += sources_density_out.name("sources_density_out")
322  .description("Concentration sources output - density of substance source, only positive part is used..")
323  .input_default("0.0")
324  .units( UnitSI().kg().s(-1).md() );
325 
326  *this += sources_sigma_out.name("sources_sigma_out")
327  .description("Concentration sources - Robin type, in_flux = sources_sigma * (sources_conc - mobile_conc).")
328  .input_default("0.0")
329  .units( UnitSI().s(-1).m(3).md() );
330 
331  *this += advection_coef.name("advection_coef")
332  .description("Advection coefficients model.")
333  .input_default("0.0")
334  .units( UnitSI().m().s(-1) );
335 
336  *this += diffusion_coef.name("diffusion_coef")
337  .description("Diffusion coefficients model.")
338  .input_default("0.0")
339  .units( UnitSI().m(2).s(-1) );
340 }
341 
342 
344 {
345  // create FieldModels
346  v_norm.set(Model<3, FieldValue<3>::Scalar>::create(fn_heat_v_norm(), flow_flux), 0.0);
347  mass_matrix_coef.set(
348  Model<3, FieldValue<3>::Scalar>::create(
349  fn_heat_mass_matrix(), cross_section, porosity, fluid_density, fluid_heat_capacity, solid_density, solid_heat_capacity
350  ),
351  0.0
352  );
353  retardation_coef.set(std::make_shared< FieldConstant<3, FieldValue<3>::Scalar> >(), 0.0);
354  sources_density_out.set(
355  Model<3, FieldValue<3>::Scalar>::create_multi(fn_heat_sources_dens(), cross_section, porosity, fluid_thermal_source, solid_thermal_source),
356  0.0
357  );
358  sources_sigma_out.set(
359  Model<3, FieldValue<3>::Scalar>::create_multi(
360  fn_heat_sources_sigma(), cross_section, porosity, fluid_density, fluid_heat_capacity, fluid_heat_exchange_rate, solid_density,
361  solid_heat_capacity, solid_heat_exchange_rate
362  ),
363  0.0
364  );
365  sources_conc_out.set(
366  Model<3, FieldValue<3>::Scalar>::create_multi(
367  fn_heat_sources_conc(), cross_section, porosity, fluid_density, fluid_heat_capacity, fluid_heat_exchange_rate, fluid_ref_temperature,
368  solid_density, solid_heat_capacity, solid_heat_exchange_rate, solid_ref_temperature, sources_sigma_out
369  ),
370  0.0
371  );
372  advection_coef.set(Model<3, FieldValue<3>::VectorFixed>::create(fn_heat_ad_coef(), fluid_density, fluid_heat_capacity, flow_flux), 0.0);
373  diffusion_coef.set(
374  Model<3, FieldValue<3>::TensorFixed>::create(
375  fn_heat_diff_coef(), flow_flux, v_norm, fluid_density, disp_l, disp_t, fluid_heat_conductivity, solid_heat_conductivity, cross_section, porosity
376  ),
377  0.0
378  );
379 }
380 
381 
383  return Selection("Heat_BC_Type", "Types of boundary conditions for heat transfer model.")
384  .add_value(bc_inflow, "inflow",
385  "Default heat transfer boundary condition.\n"
386  "On water inflow (($(q_w \\le 0)$)), total energy flux is given by the reference temperature 'bc_temperature'. "
387  "On water outflow we prescribe zero diffusive flux, "
388  "i.e. the energy flows out only due to advection.")
389  .add_value(bc_dirichlet, "dirichlet",
390  "Dirichlet boundary condition (($T = T_D $)).\n"
391  "The prescribed temperature (($T_D$)) is specified by the field 'bc_temperature'.")
392  .add_value(bc_total_flux, "total_flux",
393  "Total energy flux boundary condition.\n"
394  "The prescribed incoming total flux can have the general form (($\\delta(f_N+\\sigma_R(T_R-T) )$)), "
395  "where the absolute flux (($f_N$)) is specified by the field 'bc_flux', "
396  "the transition parameter (($\\sigma_R$)) by 'bc_robin_sigma', "
397  "and the reference temperature (($T_R$)) by 'bc_temperature'.")
398  .add_value(bc_diffusive_flux, "diffusive_flux",
399  "Diffusive flux boundary condition.\n"
400  "The prescribed incoming energy flux due to diffusion can have the general form (($\\delta(f_N+\\sigma_R(T_R-T) )$)), "
401  "where the absolute flux (($f_N$)) is specified by the field 'bc_flux', "
402  "the transition parameter (($\\sigma_R$)) by 'bc_robin_sigma', "
403  "and the reference temperature (($T_R$)) by 'bc_temperature'.")
404  .close();
405 }
406 
407 
409 {
410  // Return empty selection just to provide model specific selection name and description.
411  // The fields are added by TransportDG using an auxiliary selection.
412  return IT::Selection(
413  std::string(ModelEqData::name()) + "_DG_output_fields",
414  "Selection of output fields for Heat Transfer DG model.");
415 }
416 
417 
418 
419 IT::Record HeatTransferModel::get_input_type(const string &implementation, const string &description)
420 {
421  return IT::Record(
422  std::string(ModelEqData::name()) + "_" + implementation,
423  description + " for heat transfer.")
426  .declare_key("balance", Balance::get_input_type(), Default("{}"),
427  "Settings for computing balance.")
428  .declare_key("output_stream", OutputTime::get_input_type(), Default("{}"),
429  "Parameters of output stream.");
430 }
431 
432 
434  AdvectionProcessBase(mesh, in_rec)
435 {
436  time_ = new TimeGovernor(in_rec.val<Input::Record>("time"));
437  ASSERT( time_->is_default() == false ).error("Missing key 'time' in Heat_AdvectionDiffusion_DG.");
438 
439  output_stream_ = OutputTime::create_output_stream("heat", in_rec.val<Input::Record>("output_stream"), time().get_unit_conversion());
440  //output_stream_->add_admissible_field_names(in_rec.val<Input::Array>("output_fields"));
441 }
442 
443 
445 {
446  balance_ = std::make_shared<Balance>("energy", mesh_);
447  balance_->init_from_input(in_rec.val<Input::Record>("balance"), *time_);
448  // initialization of balance object
449  eq_data().subst_idx_ = {balance_->add_quantity("energy")};
450  balance_->units(UnitSI().m(2).kg().s(-2));
451 }
452 
453 
455 {}
456 
457 
458 
459 
#define ASSERT(expr)
Definition: asserts.hh:351
static Input::Type::Abstract & get_input_type()
Common specification of the input record for secondary equations.
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:53
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:39
TimeGovernor * time_
Definition: equation.hh:241
std::shared_ptr< Balance > balance_
object for calculation and writing the mass balance to file.
Definition: equation.hh:252
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
static constexpr Mask input_copy
Definition: field_flag.hh:44
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Definition: field_flag.hh:49
static IT::Selection get_output_selection()
Definition: heat_model.cc:408
vector< unsigned int > subst_idx_
List of indices used to call balance methods for a set of quantities.
Definition: heat_model.hh:228
static const Input::Type::Selection & get_bc_type_selection()
Definition: heat_model.cc:382
void init_balance(const Input::Record &in_rec)
Definition: heat_model.cc:444
~HeatTransferModel() override
Definition: heat_model.cc:454
HeatTransferModel(Mesh &mesh, const Input::Record in_rec)
Definition: heat_model.cc:433
std::shared_ptr< OutputTime > output_stream_
Definition: heat_model.hh:265
virtual ModelEqData & eq_data()=0
Derived class should implement getter for ModelEqData instance.
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
Record type proxy class.
Definition: type_record.hh:182
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Definition: mesh.h:362
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
static std::shared_ptr< OutputTime > create_output_stream(const std::string &equation_name, const Input::Record &in_rec, const std::shared_ptr< TimeUnitConversion > &time_unit_conv)
This method delete all object instances of class OutputTime stored in output_streams vector.
Definition: output_time.cc:187
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & J()
Returns Joule.
Definition: unit_si.cc:40
static UnitSI & W()
Returns Watt.
Definition: unit_si.cc:45
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
double Sclr
Discontinuous Galerkin method for equation of transport with dispersion.