MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
pipe.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/hydr/pipe.cc,v 1.40 2017/01/12 14:46:32 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*
33  * Copyright 1999-2000 Lamberto Puggelli <puggelli@tiscalinet.it>
34  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
35  */
36 
37 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
38 
39 #include <cfloat>
40 #include <limits>
41 
42 #include "pipe.h"
43 
44 /* Pipe - begin */
45 
46 Pipe::Pipe(unsigned int uL, const DofOwner* pDO, HydraulicFluid* hf,
47  const PressureNode* p1, const PressureNode* p2,
48  doublereal Dh, doublereal A, doublereal L, flag transition,
49  doublereal q0, flag fOut)
50 : Elem(uL, fOut),
51 HydraulicElem(uL, pDO, hf, fOut),
52 pNode1(p1), pNode2(p2),
53 diameter(Dh), area(A),
54 length(L), turbulent(transition), q0(q0)
55 {
56  ASSERT(pNode1 != NULL);
57  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
58  ASSERT(pNode2 != NULL);
59  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
60  ASSERT(Dh > std::numeric_limits<doublereal>::epsilon());
61  ASSERT(A > std::numeric_limits<doublereal>::epsilon());
62  ASSERT(L> std::numeric_limits<doublereal>::epsilon());
63  viscosity = HF->dGetViscosity();
64 
65  doublereal density = HF->dGetDensity((pNode1->dGetX()+pNode2->dGetX())/2.);
66 
67  klam = 32.*length*viscosity/((diameter*diameter)*area*density);
68  ktra = .5*length/(diameter*density*area*area);
69  doublereal den = pow(diameter, 1.25)*pow(area, 1.75)*density;
70  doublereal num = .1582*pow(viscosity, .25)*length;
71  ktrb = num/den;
72 
73 #ifdef HYDR_DEVEL
74  DEBUGCOUT("Costruttore Laminare: klam: " << klam <<std::endl);
75  DEBUGCOUT("Costruttore Transizione: ktra: " << ktra <<std::endl);
76  DEBUGCOUT("Costruttore Turbolento: ktrb: " << ktrb <<std::endl);
77 #endif /* HYDR_DEVEL */
78 }
79 
80 Pipe::~Pipe(void)
81 {
82  NO_OP;
83 }
84 
85 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
86 HydraulicElem::Type Pipe::GetHydraulicType(void) const {
87  return HydraulicElem::PIPE;
88 }
89 
90 /* Contributo al file di restart */
91 std::ostream& Pipe::Restart(std::ostream& out) const
92 {
93  return out << "Pipe not implemented yet!" << std::endl;
94 }
95 
96 unsigned int Pipe::iGetNumDof(void) const {
97  return 1;
98 }
99 
100 DofOrder::Order Pipe::GetDofType(unsigned int i) const {
101  ASSERT(i == 0);
102  return DofOrder::ALGEBRAIC;
103 }
104 
105 
106 void Pipe::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const {
107  *piNumRows = 3;
108  *piNumCols = 3;
109 }
110 
111 VariableSubMatrixHandler&
112 Pipe::AssJac(VariableSubMatrixHandler& WorkMat,
113  doublereal dCoef,
114  const VectorHandler& XCurr,
115  const VectorHandler& XPrimeCurr)
116 {
117  DEBUGCOUT("Entering Pipe::AssJac()" << std::endl);
118 #ifdef HYDR_DEVEL
119  DEBUGCOUT("dblepsilon INIZIO: " << std::numeric_limits<doublereal>::epsilon() << std::endl);
120 #endif /* HYDR_DEVEL */
121 
122  FullSubMatrixHandler& WM = WorkMat.SetFull();
123  WM.ResizeReset(3, 3);
124 
125  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
126  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
127  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
128  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
129  integer iFirstIndex = iGetFirstIndex()+1;
130 
131  WM.PutRowIndex(1, iNode1RowIndex);
132  WM.PutRowIndex(2, iNode2RowIndex);
133  WM.PutColIndex(1, iNode1ColIndex);
134  WM.PutColIndex(2, iNode2ColIndex);
135  WM.PutRowIndex(3, iFirstIndex);
136  WM.PutColIndex(3, iFirstIndex);
137 
138  doublereal p1 = pNode1->dGetX();
139  doublereal p2 = pNode2->dGetX();
140  doublereal jumpPres = fabs(p1-p2);
141  doublereal q = XCurr(iFirstIndex); /* portata */
142  doublereal Jac13 = -1.;
143  doublereal Jac23 = 1.;
144  doublereal Jac31 = 1.;
145  doublereal Jac32 = -1.;
146  doublereal Jac33;
147 
148  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
149  /****************************************
150  * moto sicuramente laminare (jacobiano)
151  ***************************************/
152 #ifdef HYDR_DEVEL
153  DEBUGCOUT("Entering Pipe::AssJac() sono laminare" << std::endl);
154 #endif /* HYDR_DEVEL */
155  Jac33 = -klam;
156  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
157  /*******************************************
158  * moto sicuramente turbolento (jacobiano)
159  ******************************************/
160 #ifdef HYDR_DEVEL
161  DEBUGCOUT("Entering Pipe::AssJac() sono turbolento" << std::endl);
162 #endif /* HYDR_DEVEL */
163  /* evito di dividere per un numero troppo piccolo */
164  if (jumpPres < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
165  jumpPres = 1.e8*std::numeric_limits<doublereal>::epsilon();
166  }
167 #ifdef HYDR_DEVEL
168  DEBUGCOUT("AssJac() JUMPPRES dopo: " << jumpPres << std::endl);
169 #endif /* HYDR_DEVEL */
170 
171  Jac33 = -7./4.*ktrb*pow(fabs(q),3./4.);
172 
173  } else {
174  /***********************************
175  * moto di transizione (jacobiano)
176  **********************************/
177 #ifdef HYDR_DEVEL
178  DEBUGCOUT("Entering Pipe::AssJac() sono in transizione" << std::endl);
179 #endif /* HYDR_DEVEL */
180  if (turbulent == 0) {
181  /*******************************************************
182  * moto di transizione laminare-turbolento (jacobiano)
183  ******************************************************/
184 #ifdef HYDR_DEVEL
185  DEBUGCOUT("Sono in transizione lam->turb" << std::endl);
186 #endif /* HYDR_DEVEL */
187  if (Re < HF->dGetRe(HydraulicFluid::LOWER)*1.25) {
188  /* uso lo jacobiano laminare */
189  Jac33 = -klam;
190  } else {
191  doublereal dva = diameter/(viscosity*area);
192  doublereal c = -1.8e-5*dva;
193  doublereal dva2 = dva*dva;
194  doublereal b = 8.e-10*dva2;
195  doublereal dva3 = dva2*dva;
196  doublereal a = 7.e-13*dva3;
197  doublereal d = .0542;
198 
199  doublereal aq = fabs(q);
200  doublereal aq2 = aq*aq;
201  doublereal aq3 = aq2*aq;
202  doublereal aq4 = aq3*aq;
203 
204  Jac33 = -5.*ktra*a*aq4-4.*ktra*b*aq3-3.*ktra*c*aq2-2.*ktra*d*aq;
205  }
206  } else {
207  /*******************************************************
208  * moto di transizione turbolento-laminare (jacobiano)
209  ******************************************************/
210 #ifdef HYDR_DEVEL
211  DEBUGCOUT("Sono in transizione turb->lam" << std::endl);
212 #endif /* HYDR_DEVEL */
213  if (Re > HF->dGetRe(HydraulicFluid::UPPER)*.775) {
214  /* uso lo jacobiano turbolento per la parte finale */
215  /* evito di dividere per un numero troppo piccolo */
216  if (jumpPres < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
217  jumpPres = 1.e8*std::numeric_limits<doublereal>::epsilon();
218  }
219 #ifdef HYDR_DEVEL
220  DEBUGCOUT("AssJac() JUMPPRES dopo: " << jumpPres << std::endl);
221 #endif /* HYDR_DEVEL */
222 
223  Jac33 = -7./4.*ktrb*pow(fabs(q),3./4.);
224 
225  } else {
226 #ifdef HYDR_DEVEL
227  DEBUGCOUT("Jac turb->lam: TRATTO DI INTERPOLAZIONE" << std::endl);
228 #endif /* HYDR_DEVEL */
229  doublereal dva = diameter/(viscosity*area);
230  doublereal c = -1.8e-5*dva;
231  doublereal dva2 = dva*dva;
232  doublereal b = 2.e-9*dva2;
233  doublereal dva3 = dva2*dva;
234  doublereal a = 9.e-13*dva3;
235  doublereal d = .0528;
236 
237  doublereal aq = fabs(q);
238  doublereal aq2 = aq*aq;
239  doublereal aq3 = aq2*aq;
240  doublereal aq4 = aq3*aq;
241 
242  Jac33 = -5.*ktra*a*aq4-4.*ktra*b*aq3-3.*ktra*c*aq2-2.*ktra*d*aq;
243  }
244  }
245  }
246 
247 #ifdef HYDR_DEVEL
248  DEBUGCOUT("JAC Re: " << Re << std::endl);
249  DEBUGCOUT("JAC density: " << HF->dGetDensity((p1+p2)/2.) << std::endl);
250  DEBUGCOUT("JAC p1: " << p1 << std::endl);
251  DEBUGCOUT("JAC p2: " << p2 << std::endl);
252  DEBUGCOUT("JAC jumpPres: " << jumpPres << std::endl);
253  DEBUGCOUT("JAC length: " << length << std::endl);
254  DEBUGCOUT("JAC diameter: " << diameter << std::endl);
255  DEBUGCOUT("JAC viscosity: " << viscosity << std::endl);
256  DEBUGCOUT("JAC turbulent: " << turbulent << std::endl);
257  DEBUGCOUT("JAC q: " << q << std::endl);
258  DEBUGCOUT("JAC Jac13: " << Jac13 << std::endl);
259  DEBUGCOUT("JAC Jac23: " << Jac23 << std::endl);
260  DEBUGCOUT("JAC Jac31: " << Jac31 << std::endl);
261  DEBUGCOUT("JAC Jac32: " << Jac32 << std::endl);
262  DEBUGCOUT("JAC Jac33: " << Jac33 << std::endl);
263 #endif /* HYDR_DEVEL */
264 
265  WM.PutCoef(1, 3, Jac13);
266  WM.PutCoef(2, 3, Jac23);
267  WM.PutCoef(3, 1, Jac31);
268  WM.PutCoef(3, 2, Jac32);
269  WM.PutCoef(3, 3, Jac33);
270 
271  return WorkMat;
272 }
273 
274 SubVectorHandler&
275 Pipe::AssRes(SubVectorHandler& WorkVec,
276  doublereal dCoef,
277  const VectorHandler& XCurr,
278  const VectorHandler& XPrimeCurr)
279 {
280  DEBUGCOUT("Entering Pipe::AssRes()" << std::endl);
281  WorkVec.Resize(3);
282 
283  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
284  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
285  integer iFirstIndex = iGetFirstIndex()+1;
286 
287  doublereal q = XCurr(iFirstIndex); /* portata q */
288  doublereal p1 = pNode1->dGetX();
289  doublereal p2 = pNode2->dGetX();
290  doublereal density = HF->dGetDensity((p1+p2)/2.);
291 
292  doublereal Res_1 = q;
293  doublereal Res_2 = -q;
294  doublereal Res_3 = 0.;
295 
296  flow = q;
297  vel = flow/(density*area);
298  Re = density*fabs(vel)*diameter/viscosity;
299  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
300  /**************************************
301  * moto sicuramente laminare (residuo)
302  *************************************/
303 #ifdef HYDR_DEVEL
304  DEBUGCOUT("Entering Pipe::AssRes() SONO LAMINARE" << std::endl);
305 #endif /* HYDR_DEVEL */
306  Res_3 = klam*q-p1+p2;
307  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
308  /*****************************************
309  * moto sicuramente turbolento (residuo)
310  ****************************************/
311 #ifdef HYDR_DEVEL
312  DEBUGCOUT("Entering Pipe::AssRes() SONO TURBOLENTO" << std::endl);
313 #endif /* HYDR_DEVEL */
314  Res_3= ktrb*copysign(pow(fabs(q),7./4.),q)-p1+p2;
315  } else {
316  /*********************************
317  * moto di transizione (residuo)
318  ********************************/
319 #ifdef HYDR_DEVEL
320  DEBUGCOUT("SONO IN TRANSIZIONE " << std::endl);
321  DEBUGCOUT("Re: " << Re << std::endl);
322  DEBUGCOUT("q: " << q << std::endl);
323  DEBUGCOUT("vel: " << vel << std::endl);
324 #endif /* HYDR_DEVEL */
325  if (turbulent == 0) {
326  /*****************************************************
327  * moto di transizione laminare-turbolento (residuo)
328  ****************************************************/
329 #ifdef HYDR_DEVEL
330  DEBUGCOUT("SONO IN TRANSIZIONE lam->turb" << std::endl);
331 #endif /* HYDR_DEVEL */
332  if (Re < HF->dGetRe(HydraulicFluid::LOWER)*1.25) {
333  /* uso il residuo laminare */
334  Res_3 = klam*q-p1+p2;
335 #ifdef HYDR_DEVEL
336  DEBUGCOUT("RES lam->turb: TRATTO 64/RE" << std::endl);
337  DEBUGCOUT("Re: " << Re << std::endl);
338 #endif /* HYDR_DEVEL */
339  } else {
340 #ifdef HYDR_DEVEL
341  DEBUGCOUT("RES lam->turb:TRATTO DI INTERPOLAZIONE"<< std::endl);
342 #endif /* HYDR_DEVEL */
343  doublereal dva = diameter/(viscosity*area);
344  doublereal c = -1.8e-5*dva;
345  doublereal dva2 = dva*dva;
346  doublereal b = 8.e-10*dva2;
347  doublereal dva3 = dva2*dva;
348  doublereal a = 7.e-13*dva3;
349  doublereal d = .0542;
350 
351 #ifdef HYDR_DEVEL
352  doublereal fa = ((a*q+b)*q+c)*q+d;
353  DEBUGCOUT("fa lam->turb: " << fa << std::endl);
354 #endif /* HYDR_DEVEL */
355  doublereal aq = fabs(q);
356  doublereal aq2 = aq*aq;
357  doublereal aq3 = aq2*aq;
358  doublereal aq4 = aq3*aq;
359 
360  Res_3 = -p1+p2+q*(ktra*a*aq4+ktra*b*aq3+ktra*c*aq2+ktra*d*aq);
361  }
362  } else {
363  /****************************************************
364  * moto di transizione turbolento-laminare (residuo)
365  ***************************************************/
366 #ifdef HYDR_DEVEL
367  DEBUGCOUT("SONO IN TRANSIZIONE turb->lam" << std::endl);
368 #endif /* HYDR_DEVEL */
369  if (Re > HF->dGetRe(HydraulicFluid::UPPER)*.775) {
370  /* utilizzo il residuo turbolento */
371  Res_3= ktrb*copysign(pow(fabs(q),7./4.),q)-p1+p2;
372 
373 #ifdef HYDR_DEVEL
374  DEBUGCOUT("RES turb->lam:TRATTO 0.3164/Re^0.25" << std::endl);
375 #endif /* HYDR_DEVEL */
376  } else {
377 #ifdef HYDR_DEVEL
378  DEBUGCOUT("RES turb->lam:TRATTO DI INTERPOLAZIONE"<< std::endl);
379 #endif /* HYDR_DEVEL */
380  doublereal dva = diameter/(viscosity*area);
381  doublereal c = -1.8e-5*dva;
382  doublereal dva2 = dva*dva;
383  doublereal b = 2.e-9*dva2;
384  doublereal dva3 = dva2*dva;
385  doublereal a = 9.e-13*dva3;
386  doublereal d = .0528;
387 
388  doublereal aq = fabs(q);
389  doublereal aq2 = aq*aq;
390  doublereal aq3 = aq2*aq;
391  doublereal aq4 = aq3*aq;
392 
393  Res_3 = -p1+p2+q*(ktra*a*aq4+ktra*b*aq3+ktra*c*aq2+ktra*d*aq);
394  }
395  }
396  }
397 
398 #ifdef HYDR_DEVEL
399  DEBUGCOUT("RES density: " << density << std::endl);
400  DEBUGCOUT("RES p1: " << p1 << std::endl);
401  DEBUGCOUT("RES p2: " << p2 << std::endl);
402  DEBUGCOUT("RES jumpPres: " << fabs(p1-p2) << std::endl);
403  DEBUGCOUT("RES length: " << length << std::endl);
404  DEBUGCOUT("RES diameter: " << diameter << std::endl);
405  DEBUGCOUT("RES viscosity: " << viscosity << std::endl);
406  DEBUGCOUT("RES area : " << area << std::endl);
407  DEBUGCOUT("RES q: " << q << std::endl);
408  DEBUGCOUT("***********************************************" << std::endl);
409  DEBUGCOUT("RES velocita': " << vel << std::endl);
410  DEBUGCOUT(" se positiva il fluido va dal nodo 1 al nodo 2" << std::endl);
411  DEBUGCOUT("***********************************************" << std::endl);
412  DEBUGCOUT("RES Re: " << Re << std::endl);
413  DEBUGCOUT("RES Res_1: " << Res_1 << std::endl);
414  DEBUGCOUT("RES Res_2: " << Res_2 << std::endl);
415  DEBUGCOUT("RES Res_3: " << Res_3 << std::endl);
416 #endif /* HYDR_DEVEL */
417 
418  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
419  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
420  WorkVec.PutItem(3, iFirstIndex, Res_3);
421 
422  return WorkVec;
423 }
424 
425 void
426 Pipe::AfterConvergence(const VectorHandler& X, const VectorHandler& XP)
427 {
428  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
429 #ifdef HYDR_DEVEL
430  DEBUGCOUT("Pipe(" << GetLabel() << "): laminar" << std::endl);
431 #endif /* HYDR_DEVEL */
432  turbulent = 0;
433  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
434 #ifdef HYDR_DEVEL
435  DEBUGCOUT("Pipe(" << GetLabel() << "): turbulent" << std::endl);
436 #endif /* HYDR_DEVEL */
437  turbulent = 1;
438 #ifdef HYDR_DEVEL
439  } else {
440  DEBUGCOUT("Pipe(" << GetLabel() << "): transition ("
441  << turbulent << ")" << std::endl);
442 #endif /* HYDR_DEVEL */
443  }
444 }
445 
446 void Pipe::Output(OutputHandler& OH) const
447 {
448  if (bToBeOutput()) {
449  OH.Hydraulic()
450  << std::setw(8) << GetLabel()
451  << " " << vel << " " << flow << " " << Re << std::endl;
452  }
453 }
454 
455 
456 void Pipe::SetValue(DataManager *pDM,
457  VectorHandler& X , VectorHandler& /* XP */ ,
458  SimulationEntity::Hints *ph)
459 {
460  integer i = iGetFirstIndex();
461  X.PutCoef(i+1, q0); /* portata iniziale nodo 2 */
462 }
463 
464 /* Pipe - end */
465 
466 
467 /* Dynamic_pipe - begin */
468 
469 Dynamic_pipe::Dynamic_pipe(unsigned int uL, const DofOwner* pDO, HydraulicFluid* hf,
470  const PressureNode* p1, const PressureNode* p2,
471  doublereal Dh,
472  doublereal A, doublereal L, flag transition,
473  doublereal q0, flag fOut)
474 : Elem(uL, fOut),
475 HydraulicElem(uL, pDO, hf, fOut),
476 pNode1(p1), pNode2(p2),
477 diameter(Dh), area(A),
478 length(L), turbulent(transition), q0(q0)
479 {
480  ASSERT(pNode1 != NULL);
481  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
482  ASSERT(pNode2 != NULL);
483  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
484  ASSERT(Dh > std::numeric_limits<doublereal>::epsilon());
485  ASSERT(A > std::numeric_limits<doublereal>::epsilon());
486  ASSERT(L > std::numeric_limits<doublereal>::epsilon());
487 
488  doublereal viscosity = HF->dGetViscosity();
489  doublereal density = HF->dGetDensity((pNode1->dGetX()+pNode2->dGetX())/2.);
490 
491  klam = 8.*length*viscosity/(diameter*diameter);
492  ktra = .5/(density*area*diameter);
493  ktrb = .5*length*.1582*pow(viscosity, .25)/(pow(diameter, 1.25)*pow(area, .75));
494 }
495 
496 Dynamic_pipe::~Dynamic_pipe(void)
497 {
498  NO_OP;
499 }
500 
501 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
502 HydraulicElem::Type Dynamic_pipe::GetHydraulicType(void) const {
503  return HydraulicElem::DYNAMIC_PIPE;
504 }
505 
506 /* Contributo al file di restart */
507 std::ostream& Dynamic_pipe::Restart(std::ostream& out) const
508 {
509  return out << "Pipe not implemented yet!" << std::endl;
510 }
511 
512 unsigned int Dynamic_pipe::iGetNumDof(void) const
513 {
514  return 3;
515 }
516 
517 DofOrder::Order
518 Dynamic_pipe::GetDofType(unsigned int i) const
519 {
520  ASSERT(i >= 0 && i <= 2);
521  return DofOrder::DIFFERENTIAL;
522 }
523 
524 void
525 Dynamic_pipe::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
526 {
527  *piNumRows = 5;
528  *piNumCols = 5;
529 }
530 
531 VariableSubMatrixHandler&
532 Dynamic_pipe::AssJac(VariableSubMatrixHandler& WorkMat,
533  doublereal dCoef,
534  const VectorHandler& XCurr,
535  const VectorHandler& XPrimeCurr)
536 {
537  DEBUGCOUT("Entering Pipe::AssJac()" << std::endl);
538  DEBUGCOUT("Valore di dblepsilon INIZIO:" << std::numeric_limits<doublereal>::epsilon() << std::endl);
539  FullSubMatrixHandler& WM = WorkMat.SetFull();
540  WM.ResizeReset(5, 5);
541 
542  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
543  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
544  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
545  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
546  integer iFirstIndex = iGetFirstIndex();
547 
548  WM.PutRowIndex(1, iNode1RowIndex);
549  WM.PutRowIndex(2, iNode2RowIndex);
550  WM.PutColIndex(1, iNode1ColIndex);
551  WM.PutColIndex(2, iNode2ColIndex);
552  WM.PutRowIndex(3, iFirstIndex+1);
553  WM.PutColIndex(3, iFirstIndex+1);
554  WM.PutRowIndex(4, iFirstIndex+2);
555  WM.PutColIndex(4, iFirstIndex+2);
556  WM.PutRowIndex(5, iFirstIndex+3);
557  WM.PutColIndex(5, iFirstIndex+3);
558 
559  doublereal p1 = pNode1->dGetX();
560  doublereal p2 = pNode2->dGetX();
561 
562  doublereal q1 = XCurr(iFirstIndex+2); /* portata nodo 1 */
563  doublereal q2 = XCurr(iFirstIndex+3); /* portata nodo 2 */
564 
565  doublereal densityDPres = HF->dGetDensityDPres();
566  doublereal densityS = HF->dGetDensity(p1); /* densita' all'inizio del tubo */
567  doublereal densityE = HF->dGetDensity(p2); /* densita' alla fine del nodo */
568  doublereal kappa1 = length*area*densityDPres;
569 
570  doublereal Jac14 = dCoef;
571  doublereal Jac25 = dCoef;
572  doublereal Jac31 = -.5;
573  doublereal Jac32 = -.5;
574  doublereal Jac33 = dCoef;
575  doublereal Jac43 = kappa1;
576  doublereal Jac44 = dCoef;
577  doublereal Jac45 = dCoef;
578  doublereal Jac51 = -area+q1*q1/(densityS*densityS*area)*densityDPres;
579  doublereal Jac52 = area-q2*q2/(densityE*densityE*area)*densityDPres;
580  doublereal Jac54 = -length/2.-2.*dCoef*(q1/(area*densityS));
581  doublereal Jac55 = length/2.+2.*dCoef*(q2/(area*densityE));
582 
583  doublereal lack54;
584  doublereal lack55;
585 
586  /* usata nel caso turbolento & transizione */
587  doublereal kappa2 = length/(6.*diameter*area);
588 
589  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
590  /*******************************************
591  * moto sicuramente laminare (jacobiano)
592  ******************************************/
593  lack54 = dCoef*(-2.*klam/densityS);
594  lack55 = dCoef*(2.*klam/densityE);
595  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
596  /*******************************************
597  * moto sicuramente turbolento (jacobiano)
598  ******************************************/
599 #ifdef HYDR_DEVEL
600  DEBUGCOUT("AssJac() sono turbolento" << std::endl);
601 #endif /* HYDR_DEVEL */
602  fa = .3164/pow(Re, .25);
603 
604  lack54 = -dCoef*(fa*kappa2/densityS)*fabs(-q2+2.*q1);
605  lack55 = dCoef*(fa*kappa2/densityE)*fabs(2.*q2-q1);
606 #ifdef HYDR_DEVEL
607  DEBUGCOUT("fa " << fa << std::endl);
608 #endif /* HYDR_DEVEL */
609  } else {
610  /***********************************
611  * moto di transizione (jacobiano)
612  **********************************/
613 #ifdef HYDR_DEVEL
614  DEBUGCOUT("AssJac() sono in transizione" << std::endl);
615 #endif /* HYDR_DEVEL */
616 
617  if (turbulent == 0) {
618  /*******************************************************
619  * moto di transizione laminare-turbolento (jacobiano)
620  ******************************************************/
621 #ifdef HYDR_DEVEL
622  DEBUGCOUT("Sono in transizione lam->turb" << std::endl);
623 #endif /* HYDR_DEVEL */
624  if (Re < HF->dGetRe(HydraulicFluid::LOWER)*1.25) {
625  /* uso lo jacobiano laminare */
626  lack54 = dCoef*(-2.*klam/densityS);
627  lack55 = dCoef*(+2.*klam/densityE);
628  } else {
629  doublereal dva = diameter/(viscosity*area);
630  doublereal c = -1.8e-5*dva;
631  doublereal dva2 = dva*dva;
632  doublereal b = 8.e-10*dva2;
633  doublereal dva3 = dva2*dva;
634  doublereal a = 7.e-13*dva3;
635  doublereal d = .0542;
636 
637  doublereal qm1 = fabs(q1);
638  doublereal qm2 = fabs(q2);
639 
640  doublereal fa1 = ((a*qm1+b)*qm1+c)*qm1+d;
641  doublereal fa2 = ((a*qm2+b)*qm2+c)*qm2+d;
642 #ifdef HYDR_DEVEL
643  DEBUGCOUT("JAC fa1: " << fa1 << std::endl);
644  DEBUGCOUT("JAC fa2: " << fa2 << std::endl);
645 #endif /* HYDR_DEVEL */
646  lack54 = -dCoef*(fa1*length/(6.*diameter*area*densityS))*fabs(-q2+2.*q1);
647  lack55 = dCoef*(fa2*length/(6.*diameter*area*densityE))*fabs(2.*q2-q1);
648  }
649  } else {
650  /*******************************************************
651  * moto di transizione turbolento-laminare (jacobiano)
652  ******************************************************/
653 #ifdef HYDR_DEVEL
654  DEBUGCOUT("Sono in transizione turb->lam" << std::endl);
655 #endif /* HYDR_DEVEL */
656 
657  if (Re > HF->dGetRe(HydraulicFluid::UPPER)*.775) {
658  /* uso lo jacobiano turbolento per la parte finale */
659  fa = .3164/pow(Re, .25);
660 
661  lack54 = -dCoef*(fa*kappa2/densityS)*fabs(-q2+2.*q1);
662  lack55 = dCoef*(fa*kappa2/densityE)*fabs(2.*q2-q1);
663 #ifdef HYDR_DEVEL
664  DEBUGCOUT("fa " << fa << std::endl);
665 #endif /* HYDR_DEVEL */
666  } else {
667 #ifdef HYDR_DEVEL
668  DEBUGCOUT("Jac turb->lam: TRATTO DI INTERPOLAZIONE"<< std::endl);
669 #endif /* HYDR_DEVEL */
670  doublereal dva = diameter/(viscosity*area);
671  doublereal c = -1.8e-5*dva;
672  doublereal dva2 = dva*dva;
673  doublereal b = 2.e-9*dva2;
674  doublereal dva3 = dva2*dva;
675  doublereal a = 9.e-13*dva3;
676  doublereal d = .0528;
677 
678  doublereal qm1 = fabs(q1);
679  doublereal qm2 = fabs(q2);
680 
681  doublereal fa1 = ((a*qm1+b)*qm1+c)*qm1+d;
682  doublereal fa2 = ((a*qm2+b)*qm2+c)*qm2+d;
683 #ifdef HYDR_DEVEL
684  DEBUGCOUT("JAC fa1: " << fa1 << std::endl);
685  DEBUGCOUT("JAC fa2: " << fa2 << std::endl);
686 #endif /* HYDR_DEVEL */
687  lack54 = -dCoef*(fa1*length/(6.*diameter*area*densityS))*fabs(-q2+2.*q1);
688  lack55 = dCoef*(fa2*length/(6.*diameter*area*densityE))*fabs(2.*q2-q1);
689  }
690  }
691  }
692 
693  Jac54 += lack54;
694  Jac55 += lack55;
695 
696 #ifdef HYDR_DEVEL
697  DEBUGCOUT("JAC Re: " << Re << std::endl);
698  DEBUGCOUT("JAC density: " << HF->dGetDensity() << std::endl);
699  DEBUGCOUT("JAC p1: " << p1 << std::endl);
700  DEBUGCOUT("JAC p2: " << p2 << std::endl);
701  DEBUGCOUT("JAC q1: " << q1 << std::endl);
702  DEBUGCOUT("JAC q2: " << q2 << std::endl);
703  doublereal q1p = XPrimeCurr(iFirstIndex+2); /* derivata q nodo 1 */
704  doublereal q2p = XPrimeCurr(iFirstIndex+3); /* derivata q nodo 2 */
705  DEBUGCOUT("JAC q1p: " << q1p << std::endl);
706  DEBUGCOUT("JAC q2p: " << q2p << std::endl);
707  DEBUGCOUT("JAC length: " << length << std::endl);
708  DEBUGCOUT("JAC diameter: " << diameter << std::endl);
709  DEBUGCOUT("JAC viscosity: " << viscosity << std::endl);
710  DEBUGCOUT("JAC turbulent: " << turbulent << std::endl);
711  DEBUGCOUT("JAC Jac14: " << Jac14 << std::endl);
712  DEBUGCOUT("JAC Jac25: " << Jac25 << std::endl);
713  DEBUGCOUT("JAC Jac31: " << Jac31 << std::endl);
714  DEBUGCOUT("JAC Jac32: " << Jac32 << std::endl);
715  DEBUGCOUT("JAC Jac33: " << Jac33 << std::endl);
716  DEBUGCOUT("JAC Jac43: " << Jac43 << std::endl);
717  DEBUGCOUT("JAC Jac44: " << Jac44 << std::endl);
718  DEBUGCOUT("JAC Jac45: " << Jac45 << std::endl);
719  DEBUGCOUT("JAC Jac51: " << Jac51 << std::endl);
720  DEBUGCOUT("JAC Jac52: " << Jac52 << std::endl);
721  DEBUGCOUT("JAC Jac54: " << Jac54 << std::endl);
722  DEBUGCOUT("JAC Jac55: " << Jac55 << std::endl);
723 #endif /* HYDR_DEVEL */
724 
725  WM.PutCoef(1, 4, Jac14);
726  WM.PutCoef(2, 5, Jac25);
727  WM.PutCoef(3, 1, Jac31);
728  WM.PutCoef(3, 2, Jac32);
729  WM.PutCoef(3, 3, Jac33);
730  WM.PutCoef(4, 3, Jac43);
731  WM.PutCoef(4, 4, Jac44);
732  WM.PutCoef(4, 5, Jac45);
733  WM.PutCoef(5, 1, Jac51);
734  WM.PutCoef(5, 2, Jac52);
735  WM.PutCoef(5, 4, Jac54);
736  WM.PutCoef(5, 5, Jac55);
737 
738  return WorkMat;
739 }
740 
741 
742 SubVectorHandler&
743 Dynamic_pipe::AssRes(SubVectorHandler& WorkVec,
744  doublereal dCoef,
745  const VectorHandler& XCurr,
746  const VectorHandler& XPrimeCurr)
747 {
748  DEBUGCOUT("Entering Dynamic_pipe::AssRes()" << std::endl);
749  WorkVec.Resize(5);
750 
751  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
752  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
753 
754  doublereal p1 = pNode1->dGetX();
755  doublereal p2 = pNode2->dGetX();
756  integer iFirstIndex = iGetFirstIndex();
757 
758  doublereal pr = XCurr(iFirstIndex+1); /* pressione */
759  doublereal prp = XPrimeCurr(iFirstIndex+1); /* derivata pressione */
760 
761  doublereal q1 = XCurr(iFirstIndex+2); /* portata nodo 1 */
762  doublereal q2 = XCurr(iFirstIndex+3); /* portata nodo 2 */
763  doublereal q1p = XPrimeCurr(iFirstIndex+2); /* derivata portata nodo 1 */
764  doublereal q2p = XPrimeCurr(iFirstIndex+3); /* derivata portata nodo 2 */
765 
766  flow1 = q1; /* per l'output */
767  flow2 = q2; /* per l'output */
768  pp = prp; /* per l'output */
769 
770  doublereal densityS = HF->dGetDensity(p1); /* densita' all'inizio del tubo */
771  doublereal densityE = HF->dGetDensity(p2); /* densita' alla fine del nodo */
772  doublereal densityM = HF->dGetDensity(pr); /* densita' a meta' tubo */
773 
774 
775  densitySt = densityS; /* per l'output */
776  densityEn = densityE; /* per l'output */
777  densityMe = densityM; /* per l'output */
778 
779  viscosity = HF->dGetViscosity();
780  densityDPres = HF->dGetDensityDPres();
781 
782  doublereal kappa1 = area*length*densityDPres;
783 
784  const doublereal x1 = -1./(sqrt(3.));
785  const doublereal x2 = -x1;
786  doublereal Qx1 = (q1+q2)*x1-q1+q2;
787  doublereal Qx2 = (q1+q2)*x2-q1+q2;
788 
789  doublereal densityx1 = HF->dGetDensity(.5*(p2-p1)*x1+pr);
790  doublereal densityx2 = HF->dGetDensity(.5*(p2-p1)*x2+pr);
791 
792  doublereal Res_1 = -q1;
793  doublereal Res_2 = -q2;
794  doublereal Res_3 = .5*(p1+p2)-pr;
795  doublereal Res_4 = -q2-q1-kappa1*prp;
796  doublereal Res_5 = (length/2.)*(q1p-q2p)+q1*q1/(area*densityS)-q2*q2/(area*densityE)+area*(p1-p2);
797  doublereal lack;
798 
799  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
800 #ifdef HYDR_DEVEL
801  DEBUGCOUT("SONO IN LAMINARE" << std::endl);
802 #endif /* HYDR_DEVEL */
803  lack = -klam*(Qx1/densityx1+Qx2/densityx2);
804 
805  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
806 #ifdef HYDR_DEVEL
807  DEBUGCOUT("SONO IN TURBOLENTO" << std::endl);
808 #endif /* HYDR_DEVEL */
809  const doublereal espo = 7./4.;
810 
811  lack = -ktrb*(copysign(pow(.5*fabs(Qx1), espo)/densityx1, Qx1)
812  +copysign(pow(.5*fabs(Qx2), espo)/densityx2, Qx2));
813 #ifdef HYDR_DEVEL
814  DEBUGCOUT("SONO IN TURBOLENTO lack: " << lack << std::endl);
815  DEBUGCOUT("SONO IN TURBOLENTO Qx1: " << Qx1 << std::endl);
816  DEBUGCOUT("SONO IN TURBOLENTO Qx2: " << Qx2 << std::endl);
817 #endif /* HYDR_DEVEL */
818  } else {
819  /*********************************
820  * moto di transizione (residuo)
821  ********************************/
822 #ifdef HYDR_DEVEL
823  DEBUGCOUT("SONO IN TRANSIZIONE" << std::endl);
824  DEBUGCOUT("Re: " << Re << std::endl);
825 #endif /* HYDR_DEVEL */
826  if (turbulent == 0) {
827  /****************************************************
828  * moto di transizione laminare-turbolento (residuo)
829  ***************************************************/
830 #ifdef HYDR_DEVEL
831  DEBUGCOUT("SONO IN TRANSIZIONE lam->turb" << std::endl);
832 #endif /* HYDR_DEVEL */
833  if (Re < HF->dGetRe(HydraulicFluid::LOWER)*1.25) {
834  /* uso il residuo laminare */
835  lack = -klam*(Qx1/densityx1+Qx2/densityx2);
836 #ifdef HYDR_DEVEL
837  DEBUGCOUT("RES lam->turb: TRATTO 64/RE" << std::endl);
838 #endif /* HYDR_DEVEL */
839  } else {
840 #ifdef HYDR_DEVEL
841  DEBUGCOUT("RES lam->turb: TRATTO DI INTERPOLAZIONE"<< std::endl);
842 #endif /* HYDR_DEVEL */
843  doublereal dva = diameter/(viscosity*area);
844  doublereal c = -1.8e-5*dva;
845  doublereal dva2 = dva*dva;
846  doublereal b = 8.e-10*dva2;
847  doublereal dva3 = dva2*dva;
848  doublereal a = 7.e-13*dva3;
849  doublereal d = .0542;
850 
851  doublereal qx1 = fabs(-q1+((q2+q1)/length)*x1); /* usata solo per calcolare il giusto fa in x1 */
852  doublereal qx2 = fabs(-q1+((q2+q1)/length)*x2); /* usata solo per calcolare il giusto fa in x2 */
853  doublereal fax1 = ((a*qx1+b)*qx1+c)*qx1+d;
854  doublereal fax2 = ((a*qx2+b)*qx2+c)*qx2+d;
855 
856  doublereal kappa3 = length/(16.*diameter*area);
857 #ifdef HYDR_DEVEL
858  DEBUGCOUT("RES lam->turb fax1: " << fax1 << std::endl);
859  DEBUGCOUT("RES lam->turb fax2: " << fax2 << std::endl);
860 #endif /* HYDR_DEVEL */
861  lack = -kappa3*(fax1*copysign(Qx1*Qx1, Qx1)/densityx1
862  +fax2*copysign(Qx2*Qx2, Qx2)/densityx2);
863  }
864  } else {
865  /****************************************************
866  * moto di transizione turbolento-laminare (residuo)
867  ***************************************************/
868 #ifdef HYDR_DEVEL
869  DEBUGCOUT("SONO IN TRANSIZIONE turb->lam" << std::endl);
870 #endif /* HYDR_DEVEL */
871  if (Re > HF->dGetRe(HydraulicFluid::UPPER)*.775) {
872  /* utilizzo il residuo turbolento */
873  const doublereal espo = 7./4.;
874 
875  lack = -ktrb*(copysign(pow(.5*fabs(Qx1), espo)/densityx1, Qx1)
876  +copysign(pow(.5*fabs(Qx2), espo)/densityx2, Qx2));
877 
878  /* DEBUGCOUT("RES turb->lam: TRATTO 0.3164/Re^0.25" << std::endl); */
879  } else {
880 #ifdef HYDR_DEVEL
881  DEBUGCOUT("RES turb->lam:TRATTO DI INTERPOLAZIONE" << std::endl);
882 #endif /* HYDR_DEVEL */
883  doublereal dva = diameter/(viscosity*area);
884  doublereal c = -1.8e-5*dva;
885  doublereal dva2 = dva*dva;
886  doublereal b = 2.e-9*dva2;
887  doublereal dva3 = dva2*dva;
888  doublereal a = 9.e-13*dva3;
889  doublereal d = .0528;
890 
891  doublereal qx1 = fabs(-q1+((q2+q1)/length)*x1); /* usata solo per calcolare il giusto fa in x1 */
892  doublereal qx2 = fabs(-q1+((q2+q1)/length)*x2); /* usata solo per calcolare il giusto fa in x2 */
893  doublereal fax1 = ((a*qx1+b)*qx1+c)*qx1+d;
894  doublereal fax2 = ((a*qx2+b)*qx2+c)*qx2+d;
895 #ifdef HYDR_DEVEL
896  DEBUGCOUT("RES turb->lam fax1: " << fax1 << std::endl);
897  DEBUGCOUT("RES turb->lam fax2: " << fax2 << std::endl);
898 #endif /* HYDR_DEVEL */
899  doublereal kappa3 = length/(16.*diameter*area);
900 
901  lack = -kappa3*(fax1*copysign(Qx1*Qx1, Qx1)/densityx1
902  +fax2*copysign(Qx2*Qx2, Qx2)/densityx2);
903  }
904  }
905  }
906 
907  Res_5 += lack; /* aggiungo la parte dovuta all'attrito */
908  VelS = -q1/(area*densityS); /* velocita' inizio */
909  VelM = ((-q1+q2)/2.)/(area*densityM); /* velocita' media */
910  VelE = q2/(area*densityE); /* velocita' fine */
911  Re = densityM*fabs(VelM)*diameter/viscosity; /* numero di Reynolds medio */
912 
913 #ifdef HYDR_DEVEL
914  DEBUGCOUT("RES density: " << HF->dGetDensity() << std::endl);
915  DEBUGCOUT("RES densityS: " << densityS << std::endl);
916  DEBUGCOUT("RES densityM: " << densityM << std::endl);
917  DEBUGCOUT("RES densityE: " << densityE << std::endl);
918  DEBUGCOUT("RES p1 : " << p1 << std::endl);
919  DEBUGCOUT("RES p2 : " << p2 << std::endl);
920  DEBUGCOUT("RES pr : " << pr << std::endl);
921  DEBUGCOUT("RES prp: " << prp << std::endl);
922  DEBUGCOUT("RES q1: " << q1 << std::endl);
923  DEBUGCOUT("RES q2: " << q2 << std::endl);
924  DEBUGCOUT("RES q1p: " << q1p << std::endl);
925  DEBUGCOUT("RES q2p: " << q2p << std::endl);
926  DEBUGCOUT("RES length: " << length << std::endl);
927  DEBUGCOUT("RES diameter: " << diameter << std::endl);
928  DEBUGCOUT("RES viscosity: " << viscosity << std::endl);
929  DEBUGCOUT("RES area: " << area << std::endl);
930  DEBUGCOUT("RES Re: " << Re << std::endl);
931  DEBUGCOUT("***********************************************"<< std::endl);
932  DEBUGCOUT("RES velocita' Start: " << VelS << std::endl);
933  DEBUGCOUT("RES velocita' Middle: " << VelM << std::endl);
934  DEBUGCOUT("RES velocita' End: " << VelE << std::endl);
935  DEBUGCOUT(" se positiva il fluido va dal nodo 1 al nodo 2 " << std::endl);
936  DEBUGCOUT("***********************************************"<< std::endl);
937  DEBUGCOUT("RES Res_1: " << Res_1 << std::endl);
938  DEBUGCOUT("RES Res_2: " << Res_2 << std::endl);
939  DEBUGCOUT("RES Res_3: " << Res_3 << std::endl);
940  DEBUGCOUT("RES Res_4: " << Res_4 << std::endl);
941  DEBUGCOUT("RES Res_5: " << Res_5 << std::endl);
942 #endif /* HYDR_DEVEL */
943 
944  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
945  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
946  WorkVec.PutItem(3, iFirstIndex+1, Res_3);
947  WorkVec.PutItem(4, iFirstIndex+2, Res_4);
948  WorkVec.PutItem(5, iFirstIndex+3, Res_5);
949 
950  return WorkVec;
951 }
952 
953 void
954 Dynamic_pipe::AfterConvergence(const VectorHandler& X, const VectorHandler& XP)
955 {
956  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
957  turbulent = 0;
958  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
959  turbulent = 1;
960  }
961 
962 #ifdef HYDR_DEVEL
963  DEBUGCOUT("Dynamic_Pipe(" << GetLabel() << "): turbulent mode = "
964  << turbulent << std::endl);
965 #endif /* HYDR_DEVEL */
966 }
967 
968 void Dynamic_pipe::Output(OutputHandler& OH) const
969 {
970  if (bToBeOutput()) {
971  std::ostream& out = OH.Hydraulic();
972  out
973  << std::setw(8) << GetLabel()
974  << " " << Re << " " << -flow1 << " " << flow2
975  << " " << densitySt << " " << densityMe << " " << densityEn
976  << " " << VelS << " " << VelM << " " << VelE
977  << " " << pp
978  << std::endl;
979  }
980 }
981 
982 void
983 Dynamic_pipe::SetValue(DataManager *pDM,
984  VectorHandler& X , VectorHandler& XP,
985  SimulationEntity::Hints *ph)
986 {
987  integer i = iGetFirstIndex();
988 
989  doublereal p1 = pNode1->dGetX();
990  doublereal p2 = pNode2->dGetX();
991 
992  X.PutCoef(i+1, .5*(p1+p2));
993  X.PutCoef(i+2, q0);
994  X.PutCoef(i+3, -q0);
995 
996  XP.PutCoef(i+1, 0.);
997  XP.PutCoef(i+2, 0.);
998  XP.PutCoef(i+3, 0.);
999 }
1000 
1001 /* Dynamic_pipe - end */
1002 
1003 
1004 /* DynamicPipe - begin */
1005 
1006 DynamicPipe::DynamicPipe(unsigned int uL,
1007  const DofOwner* pDO,
1008  HydraulicFluid* hf,
1009  const PressureNode* p1,
1010  const PressureNode* p2,
1011  doublereal Dh,
1012  doublereal A,
1013  doublereal L,
1014  flag transition,
1015  doublereal q0,
1016  flag fOut)
1017 : Elem(uL, fOut),
1018 HydraulicElem(uL, pDO, hf, fOut),
1019 pNode1(p1),
1020 pNode2(p2),
1021 diameter(Dh),
1022 area(A),
1023 length(L),
1024 turbulent(transition),
1025 q0(q0)
1026 {
1027  ASSERT(pNode1 != NULL);
1028  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
1029  ASSERT(pNode2 != NULL);
1030  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
1031 
1032  ASSERT(Dh > std::numeric_limits<doublereal>::epsilon());
1033  ASSERT(A > std::numeric_limits<doublereal>::epsilon());
1034  ASSERT(L > std::numeric_limits<doublereal>::epsilon());
1035 
1036  dKlam = 16.*length/(diameter*diameter);
1037  dKtra = .25*length/(diameter*area);
1038  dKtrb = .5*length*.1582/(pow(diameter, 1.25)*pow(area, .75));
1039 }
1040 
1041 DynamicPipe::~DynamicPipe(void)
1042 {
1043  NO_OP;
1044 }
1045 
1046 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
1047 HydraulicElem::Type DynamicPipe::GetHydraulicType(void) const {
1048  return HydraulicElem::DYNAMIC_PIPE;
1049 }
1050 
1051 /* Contributo al file di restart */
1052 std::ostream& DynamicPipe::Restart(std::ostream& out) const
1053 {
1054  return out << "Pipe not implemented yet!" << std::endl;
1055 }
1056 
1057 unsigned int DynamicPipe::iGetNumDof(void) const
1058 {
1059  return 4;
1060 }
1061 
1062 DofOrder::Order
1063 DynamicPipe::GetDofType(unsigned int i) const
1064 {
1065  ASSERT(i >= 0 && i < 4);
1066  return DofOrder::DIFFERENTIAL;
1067 }
1068 
1069 DofOrder::Order
1070 DynamicPipe::GetEqType(unsigned int i) const
1071 {
1072  ASSERT(i >= 0 && i < 4);
1073  return DofOrder::DIFFERENTIAL;
1074 }
1075 
1076 void
1077 DynamicPipe::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
1078 {
1079  *piNumRows = 6;
1080  *piNumCols = 6;
1081 }
1082 
1083 VariableSubMatrixHandler&
1084 DynamicPipe::AssJac(VariableSubMatrixHandler& WorkMat,
1085  doublereal dCoef,
1086  const VectorHandler& XCurr,
1087  const VectorHandler& XPrimeCurr)
1088 {
1089  FullSubMatrixHandler& WM = WorkMat.SetFull();
1090  WM.ResizeReset(6, 6);
1091 
1092  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1093  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1094  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
1095  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
1096 
1097  integer iFirstIndex = iGetFirstIndex();
1098 
1099  WM.PutRowIndex(1, iNode1RowIndex);
1100  WM.PutRowIndex(2, iNode2RowIndex);
1101  WM.PutRowIndex(3, iFirstIndex+1);
1102  WM.PutRowIndex(4, iFirstIndex+2);
1103  WM.PutRowIndex(5, iFirstIndex+3);
1104  WM.PutRowIndex(6, iFirstIndex+4);
1105 
1106  WM.PutColIndex(1, iNode1ColIndex);
1107  WM.PutColIndex(2, iNode2ColIndex);
1108  WM.PutColIndex(3, iFirstIndex+1);
1109  WM.PutColIndex(4, iFirstIndex+2);
1110  WM.PutColIndex(5, iFirstIndex+3);
1111  WM.PutColIndex(6, iFirstIndex+4);
1112 
1113  doublereal dRDP1 = HF->dGetDensityDPres(p1);
1114  doublereal dRDP2 = HF->dGetDensityDPres(p2);
1115  doublereal dRDP0 = HF->dGetDensityDPres(.5*(p1+p2));
1116 
1117  doublereal qq1 = .75*q1+.25*q2;
1118  doublereal qq2 = .25*q1+.75*q2;
1119 
1120  doublereal dd1 = .5*(density1+density0);
1121  doublereal dd2 = .5*(density0+density2);
1122 
1123  doublereal dLoss11 = 0.;
1124  doublereal dLoss12 = 0.;
1125  doublereal dLoss21 = 0.;
1126  doublereal dLoss22 = 0.;
1127 
1128  if ((Re < HF->dGetRe(HydraulicFluid::LOWER)) || ((turbulent == 0)
1129  && (Re < 1.25*HF->dGetRe(HydraulicFluid::LOWER)))) {
1130  /* laminar, or ascending transition ascendente (continuation) */
1131  doublereal dk = dKlam*viscosity*dCoef;
1132 
1133  dLoss11 = .75*dk/dd1;
1134  dLoss12 = .25*dk/dd1;
1135  dLoss21 = .25*dk/dd2;
1136  dLoss22 = .75*dk/dd2;
1137  } else if ((Re > HF->dGetRe(HydraulicFluid::UPPER)) || ((turbulent == 1)
1138  && (Re > .775*HF->dGetRe(HydraulicFluid::UPPER)))) {
1139  /* turbulent, or descending transition (continuation) */
1140 
1141  doublereal dk = dKtrb*1.75*pow(viscosity, .25)*dCoef;
1142  doublereal dl1 = dk*pow(fabs(qq1), .75)/dd1;
1143  doublereal dl2 = dk*pow(fabs(qq2), .75)/dd2;
1144 
1145  dLoss11 = .75*dl1;
1146  dLoss12 = .25*dl1;
1147  dLoss21 = .25*dl2;
1148  dLoss22 = .75*dl2;
1149  } else {
1150  /* transition */
1151  doublereal dva = diameter/(viscosity*area);
1152  doublereal dva2 = dva*dva;
1153  doublereal dva3 = dva2*dva;
1154 
1155  doublereal a, b, c, d;
1156 
1157  if (turbulent == 0) {
1158  /* ascending */
1159  a = 7.e-13*dva3;
1160  b = 8.e-10*dva2;
1161  c = -1.8e-5*dva;
1162  d = .0542;
1163  } else /* if (turbulent == 1) */ {
1164  /* descending */
1165  a = 9.e-13*dva3;
1166  b = 2.e-9*dva2;
1167  c = -1.8e-5*dva;
1168  d = .0528;
1169  }
1170 
1171  doublereal fqq1 = fabs(qq1);
1172  doublereal fqq2 = fabs(qq2);
1173 
1174  doublereal fax1 = (((5.*a*fqq1+4.*b)*fqq1+3.*c)*fqq1+2.*d)*fqq1;
1175  doublereal fax2 = (((5.*a*fqq2+4.*b)*fqq2+3.*c)*fqq2+2.*d)*fqq2;
1176 
1177  doublereal dk = dKtrb*dCoef;
1178 
1179  doublereal dl1 = dk*fax1/dd1;
1180  doublereal dl2 = dk*fax2/dd2;
1181 
1182  dLoss11 = .75*dl1;
1183  dLoss12 = .25*dl1;
1184  dLoss21 = .25*dl2;
1185  dLoss22 = .75*dl2;
1186  }
1187 
1188  /* primo blocco: conservazione massa */
1189  doublereal dr = .5*dCoef;
1190  WM.PutCoef(1, 5, -dr);
1191  WM.PutCoef(1, 6, -dr);
1192  WM.PutCoef(2, 5, dr);
1193  WM.PutCoef(2, 6, dr);
1194 
1195  dr = area*length/8.;
1196  WM.PutCoef(1, 3, -dr*3.*densityDPres1);
1197  WM.PutCoef(1, 4, -dr*densityDPres1);
1198  WM.PutCoef(2, 3, -dr*densityDPres2);
1199  WM.PutCoef(2, 4, -dr*3.*densityDPres2);
1200 
1201 
1202  /* secondo blocco: equazione quantita' di moto */
1203  dr = .5*dCoef*area;
1204  doublereal q12 = .5*(q1+q2);
1205  doublereal ddq12 = .5*(q12*q12)/(density0*density0*area)*dRDP0*dCoef;
1206  doublereal ddq1 = q1*q1/(density1*density1*area)*dRDP1*dCoef;
1207  doublereal ddq2 = q2*q2/(density2*density2*area)*dRDP2*dCoef;
1208  WM.PutCoef(3, 3, -dr + ddq1 - ddq12);
1209  WM.PutCoef(3, 4, dr - ddq12);
1210  WM.PutCoef(4, 3, -dr + ddq12);
1211  WM.PutCoef(4, 4, dr + ddq12 - ddq2);
1212 
1213  /* manca la viscosita' */
1214  dr = length/8.;
1215  doublereal dq12 = .5*(q1+q2)/density0;
1216  WM.PutCoef(3, 5, 3.*dr+(dq12-2.*q1/density1)/area*dCoef + dLoss11);
1217  WM.PutCoef(3, 6, dr+dq12/area*dCoef + dLoss12);
1218  WM.PutCoef(4, 5, dr-dq12/area*dCoef + dLoss21);
1219  WM.PutCoef(4, 6, 3.*dr+(2.*q2/density2-dq12)/area*dCoef + dLoss22);
1220 
1221  /* terzo blocco: definizione delle pressioni nodali interne */
1222  WM.PutCoef(5, 1, -1.);
1223  WM.PutCoef(6, 2, -1.);
1224 
1225  WM.PutCoef(5, 3, dCoef);
1226  WM.PutCoef(6, 4, dCoef);
1227 
1228  return WorkMat;
1229 }
1230 
1231 
1232 SubVectorHandler&
1233 DynamicPipe::AssRes(SubVectorHandler& WorkVec,
1234  doublereal dCoef,
1235  const VectorHandler& XCurr,
1236  const VectorHandler& XPrimeCurr)
1237 {
1238  DEBUGCOUT("Entering DynamicPipe::AssRes()" << std::endl);
1239  WorkVec.Resize(6);
1240 
1241  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1242  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1243 
1244  doublereal pn1 = pNode1->dGetX();
1245  doublereal pn2 = pNode2->dGetX();
1246 
1247  integer iFirstIndex = iGetFirstIndex();
1248 
1249  p1 = XCurr(iFirstIndex+1); /* pressione */
1250  p2 = XCurr(iFirstIndex+2); /* pressione */
1251  p1p = XPrimeCurr(iFirstIndex+1); /* derivata pressione */
1252  p2p = XPrimeCurr(iFirstIndex+2); /* derivata pressione */
1253 
1254  q1 = XCurr(iFirstIndex+3); /* portata nodo 1 */
1255  q2 = XCurr(iFirstIndex+4); /* portata nodo 2 */
1256  q1p = XPrimeCurr(iFirstIndex+3); /* derivata portata nodo 1 */
1257  q2p = XPrimeCurr(iFirstIndex+4); /* derivata portata nodo 2 */
1258 
1259  doublereal p0 = .5*(p1+p2);
1260 
1261  density1 = HF->dGetDensity(p1); /* densita' all'inizio del tubo */
1262  density2 = HF->dGetDensity(p2); /* densita' alla fine del nodo */
1263  density0 = HF->dGetDensity(p0); /* densita' a meta' tubo */
1264 
1265  /* viscosity = HF->dGetViscosity(); */
1266  densityDPres1 = HF->dGetDensityDPres(.75*p1+.25*p2);
1267  densityDPres2 = HF->dGetDensityDPres(.25*p1+.75*p2);
1268 
1269  viscosity = HF->dGetViscosity(p0);
1270 
1271  doublereal qq1 = .75*q1+.25*q2;
1272  doublereal qq2 = .25*q1+.75*q2;
1273 
1274  doublereal dd1 = .5*(density1+density0);
1275  doublereal dd2 = .5*(density0+density2);
1276 
1277  doublereal q12 = .5*(q1+q2);
1278 
1279  /* determines whether the flow will be turbulent or laminar */
1280  Re = fabs(q12)/area*diameter/viscosity;
1281 
1282  doublereal dLoss1 = 0.;
1283  doublereal dLoss2 = 0.;
1284 
1285  if ((Re < HF->dGetRe(HydraulicFluid::LOWER)) || ((turbulent == 0)
1286  && (Re < 1.25*HF->dGetRe(HydraulicFluid::LOWER)))) {
1287  /* laminar, or ascending transition (continuation) */
1288  doublereal dk = dKlam*viscosity;
1289 
1290  dLoss1 = dk*qq1/dd1;
1291  dLoss2 = dk*qq2/dd2;
1292  } else if ((Re > HF->dGetRe(HydraulicFluid::UPPER)) || ((turbulent == 1)
1293  && (Re > .775*HF->dGetRe(HydraulicFluid::UPPER)))) {
1294  /* turbulent, or descending transition (continuation) */
1295  doublereal dk = dKtrb*pow(viscosity, .25);
1296 
1297  dLoss1 = dk*qq1*pow(fabs(qq1), .75)/dd1;
1298  dLoss2 = dk*qq2*pow(fabs(qq2), .75)/dd2;
1299  } else {
1300  /* transition */
1301  doublereal dva = diameter/(viscosity*area);
1302  doublereal dva2 = dva*dva;
1303  doublereal dva3 = dva2*dva;
1304 
1305  doublereal a, b, c, d;
1306 
1307  if (turbulent == 0) {
1308  /* ascending */
1309  a = 7.e-13*dva3;
1310  b = 8.e-10*dva2;
1311  c = -1.8e-5*dva;
1312  d = .0542;
1313  } else /* if (turbulent == 1) */ {
1314  /* descending */
1315  a = 9.e-13*dva3;
1316  b = 2.e-9*dva2;
1317  c = -1.8e-5*dva;
1318  d = .0528;
1319  }
1320 
1321  doublereal fqq1 = fabs(qq1);
1322  doublereal fqq2 = fabs(qq2);
1323 
1324  doublereal fax1 = qq1*(((a*fqq1+b)*fqq1+c)*fqq1+d)*fqq1;
1325  doublereal fax2 = qq2*(((a*fqq2+b)*fqq2+c)*fqq2+d)*fqq2;
1326 
1327  dLoss1 = dKtra*fax1/dd1;
1328  dLoss2 = dKtra*fax2/dd2;
1329  }
1330 
1331  /* mass conservation */
1332  doublereal dr = area*length/8.;
1333  WorkVec.PutItem(1, iNode1RowIndex , q12+dr*densityDPres1*(3.*p1p+p2p));
1334  WorkVec.PutItem(2, iNode2RowIndex , -q12+dr*densityDPres2*(p1p+3.*p2p));
1335 
1336  /* momentum balance */
1337  doublereal dp = .5*area*(p2-p1);
1338  doublereal dq12 = q12*q12/(density0*area);
1339  dr = length/8.;
1340  WorkVec.PutItem(3, iFirstIndex+1,
1341  -dr*(3.*q1p+q2p)-dq12+q1*q1/(density1*area)-dp-dLoss1);
1342  WorkVec.PutItem(4, iFirstIndex+2,
1343  -dr*(q1p+3.*q2p)-q2*q2/(density2*area)+dq12-dp-dLoss2);
1344 
1345  /* differential pressure definition */
1346  WorkVec.PutItem(5, iFirstIndex+3, pn1-p1);
1347  WorkVec.PutItem(6, iFirstIndex+4, pn2-p2);
1348 
1349  return WorkVec;
1350 }
1351 
1352 void
1353 DynamicPipe::AfterConvergence(const VectorHandler& X, const VectorHandler& XP)
1354 {
1355  if (Re < HF->dGetRe(HydraulicFluid::LOWER)) {
1356  turbulent = 0;
1357  } else if (Re > HF->dGetRe(HydraulicFluid::UPPER)) {
1358  turbulent = 1;
1359  }
1360 
1361 #ifdef HYDR_DEVEL
1362  DEBUGCOUT("DynamicPipe(" << GetLabel() << "): turbulent mode = "
1363  << turbulent << std::endl);
1364 #endif /* HYDR_DEVEL */
1365 }
1366 
1367 void DynamicPipe::Output(OutputHandler& OH) const
1368 {
1369  if (bToBeOutput()) {
1370  std::ostream& out = OH.Hydraulic();
1371  out
1372  << std::setw(8) << GetLabel() /* 1 */
1373  << " " << p1 /* 2 */
1374  << " " << p2 /* 3 */
1375  << " " << p1p /* 4 */
1376  << " " << p2p /* 5 */
1377  << " " << q1 /* 6 */
1378  << " " << q2 /* 7 */
1379  << " " << q1p /* 8 */
1380  << " " << q2p /* 9 */
1381  << " " << density1 /* 10 */
1382  << " " << density0 /* 11 */
1383  << " " << density2 /* 12 */
1384  << " " << densityDPres1 /* 13 */
1385  << " " << densityDPres2 /* 14 */
1386  << " " << Re /* 15 */
1387  << " " << turbulent /* 16 */
1388  << std::endl;
1389  }
1390 }
1391 
1392 void
1393 DynamicPipe::SetValue(DataManager *pDM,
1394  VectorHandler& X , VectorHandler& XP,
1395  SimulationEntity::Hints *ph)
1396 {
1397  integer i = iGetFirstIndex();
1398 
1399  doublereal p1 = pNode1->dGetX();
1400  doublereal p2 = pNode2->dGetX();
1401 
1402  X.PutCoef(i+1, p1);
1403  X.PutCoef(i+2, p2);
1404  X.PutCoef(i+3, q0);
1405  X.PutCoef(i+4, -q0);
1406 
1407  XP.PutCoef(i+1, 0.);
1408  XP.PutCoef(i+2, 0.);
1409  XP.PutCoef(i+3, 0.);
1410  XP.PutCoef(i+4, 0.);
1411 }
1412 
1413 /* DynamicPipe - end */
Dynamic_pipe::length
doublereal length
Definition: pipe.h:123
Dynamic_pipe::pNode1
const PressureNode * pNode1
Definition: pipe.h:119
Pipe::pNode1
const PressureNode * pNode1
Definition: pipe.h:47
DynamicPipe::pNode1
const PressureNode * pNode1
Definition: pipe.h:203
DynamicPipe::p2p
doublereal p2p
Definition: pipe.h:218
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
DynamicPipe::q2
doublereal q2
Definition: pipe.h:221
DynamicPipe::q2p
doublereal q2p
Definition: pipe.h:223
Dynamic_pipe::VelE
doublereal VelE
Definition: pipe.h:142
DynamicPipe::dKtrb
doublereal dKtrb
Definition: pipe.h:234
flag
long int flag
Definition: mbdyn.h:43
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
Pipe::vel
doublereal vel
Definition: pipe.h:56
Pipe::turbulent
flag turbulent
Definition: pipe.h:53
Dynamic_pipe::Dynamic_pipe
Dynamic_pipe(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, doublereal Dh, doublereal A, doublereal L, flag transition, doublereal q0, flag fOut)
Definition: pipe.cc:469
Dynamic_pipe::densityEn
doublereal densityEn
Definition: pipe.h:135
HydraulicFluid::dGetRe
virtual doublereal dGetRe(Re which)
Definition: hfluid.h:95
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
Dynamic_pipe::~Dynamic_pipe
~Dynamic_pipe(void)
Definition: pipe.cc:496
HydraulicFluid::dGetDensityDPres
virtual doublereal dGetDensityDPres(void) const =0
Dynamic_pipe::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: pipe.cc:954
Dynamic_pipe::diameter
doublereal diameter
Definition: pipe.h:121
Pipe::Output
virtual void Output(OutputHandler &OH) const
Definition: pipe.cc:446
Dynamic_pipe::fa
doublereal fa
Definition: pipe.h:143
DynamicPipe::p1p
doublereal p1p
Definition: pipe.h:217
Dynamic_pipe::Output
virtual void Output(OutputHandler &OH) const
Definition: pipe.cc:968
DynamicPipe::q0
doublereal q0
Definition: pipe.h:211
p0
static double * p0
Definition: test_strext_socket_lib.c:100
FullSubMatrixHandler::PutCoef
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
Dynamic_pipe::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: pipe.cc:512
Dynamic_pipe::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: pipe.cc:983
DynamicPipe::q1
doublereal q1
Definition: pipe.h:220
DynamicPipe::density1
doublereal density1
Definition: pipe.h:225
Pipe::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:275
NO_OP
#define NO_OP
Definition: myassert.h:74
Pipe::area
doublereal area
Definition: pipe.h:51
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
Pipe::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: pipe.cc:456
DynamicPipe::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: pipe.cc:1052
DynamicPipe::length
doublereal length
Definition: pipe.h:208
DynamicPipe::DynamicPipe
DynamicPipe(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, doublereal Dh, doublereal A, doublereal L, flag transition, doublereal q0, flag fOut)
Definition: pipe.cc:1006
DynamicPipe::diameter
doublereal diameter
Definition: pipe.h:206
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
SubVectorHandler::PutItem
virtual void PutItem(integer iSubRow, integer iRow, const doublereal &dCoef)
Definition: submat.h:1445
DynamicPipe::GetEqType
virtual DofOrder::Order GetEqType(unsigned int i) const
Definition: pipe.cc:1070
DynamicPipe::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: pipe.cc:1077
DynamicPipe::densityDPres1
doublereal densityDPres1
Definition: pipe.h:228
Pipe::ktra
doublereal ktra
Definition: pipe.h:61
Pipe::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: pipe.cc:100
HydraulicFluid::dGetDensity
virtual doublereal dGetDensity(void) const =0
ScalarAlgebraicNode::dGetX
virtual const doublereal & dGetX(void) const
Definition: node.h:492
DynamicPipe::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: pipe.cc:1393
DynamicPipe::Re
doublereal Re
Definition: pipe.h:213
Pipe::q0
doublereal q0
Definition: pipe.h:54
Pipe::Pipe
Pipe(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, doublereal Dh, doublereal A, doublereal L, flag transition, doublereal q0, flag fOut)
Definition: pipe.cc:46
Dynamic_pipe::area
doublereal area
Definition: pipe.h:122
DynamicPipe::q1p
doublereal q1p
Definition: pipe.h:222
HydraulicElem::HF
HydraulicFluid * HF
Definition: preselem.h:79
DynamicPipe::dKtra
doublereal dKtra
Definition: pipe.h:235
Pipe::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: pipe.cc:86
DynamicPipe::dKlam
doublereal dKlam
Definition: pipe.h:233
Pipe::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: pipe.cc:106
Pipe::~Pipe
~Pipe(void)
Definition: pipe.cc:80
DynamicPipe::pNode2
const PressureNode * pNode2
Definition: pipe.h:204
DynamicPipe::densityDPres2
doublereal densityDPres2
Definition: pipe.h:229
Dynamic_pipe::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: pipe.cc:525
DynamicPipe::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: pipe.cc:1057
copysign
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
Dynamic_pipe::VelM
doublereal VelM
Definition: pipe.h:141
Dynamic_pipe::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:743
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
DynamicPipe::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: pipe.cc:1047
Node::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: node.cc:82
Dynamic_pipe::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: pipe.cc:507
Pipe::length
doublereal length
Definition: pipe.h:52
Dynamic_pipe::klam
doublereal klam
Definition: pipe.h:131
Dynamic_pipe::turbulent
flag turbulent
Definition: pipe.h:124
Dynamic_pipe::Re
doublereal Re
Definition: pipe.h:129
Dynamic_pipe::densitySt
doublereal densitySt
Definition: pipe.h:133
Pipe::ktrb
doublereal ktrb
Definition: pipe.h:59
Dynamic_pipe::viscosity
doublereal viscosity
Definition: pipe.h:138
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
Pipe::pNode2
const PressureNode * pNode2
Definition: pipe.h:48
HydraulicFluid::dGetViscosity
virtual doublereal dGetViscosity(void) const =0
Dynamic_pipe::pNode2
const PressureNode * pNode2
Definition: pipe.h:120
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
Pipe::diameter
doublereal diameter
Definition: pipe.h:49
Pipe::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: pipe.cc:91
Pipe::flow
doublereal flow
Definition: pipe.h:57
DynamicPipe::~DynamicPipe
~DynamicPipe(void)
Definition: pipe.cc:1041
DynamicPipe::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:1233
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
Pipe::viscosity
doublereal viscosity
Definition: pipe.h:50
Dynamic_pipe::pp
doublereal pp
Definition: pipe.h:145
Pipe::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: pipe.cc:426
Pipe::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:112
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
Dynamic_pipe::flow1
doublereal flow1
Definition: pipe.h:127
DynamicPipe::viscosity
doublereal viscosity
Definition: pipe.h:231
DynamicPipe::turbulent
flag turbulent
Definition: pipe.h:210
Dynamic_pipe::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: pipe.cc:518
Elem
Definition: elem.h:75
Dynamic_pipe::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: pipe.cc:502
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
DynamicPipe::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: pipe.cc:1353
Dynamic_pipe::q0
doublereal q0
Definition: pipe.h:125
Dynamic_pipe::densityDPres
doublereal densityDPres
Definition: pipe.h:136
Dynamic_pipe::densityMe
doublereal densityMe
Definition: pipe.h:134
DynamicPipe::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:1084
PressureNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: presnode.h:54
DynamicPipe::density2
doublereal density2
Definition: pipe.h:227
DynamicPipe::area
doublereal area
Definition: pipe.h:207
Pipe::Re
doublereal Re
Definition: pipe.h:58
DynamicPipe::density0
doublereal density0
Definition: pipe.h:226
Dynamic_pipe::flow2
doublereal flow2
Definition: pipe.h:128
a
static const doublereal a
Definition: hfluid_.h:289
OutputHandler::Hydraulic
std::ostream & Hydraulic(void) const
Definition: output.h:492
Dynamic_pipe::ktra
doublereal ktra
Definition: pipe.h:132
DynamicPipe::p2
doublereal p2
Definition: pipe.h:216
Dynamic_pipe::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: pipe.cc:532
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
Dynamic_pipe::VelS
doublereal VelS
Definition: pipe.h:140
doublereal
double doublereal
Definition: colamd.c:52
Pipe::klam
doublereal klam
Definition: pipe.h:60
integer
long int integer
Definition: colamd.c:51
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Dynamic_pipe::ktrb
doublereal ktrb
Definition: pipe.h:130
Pipe::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: pipe.cc:96
DynamicPipe::p1
doublereal p1
Definition: pipe.h:215
VectorHandler::Resize
virtual void Resize(integer iNewSize)=0
DynamicPipe::Output
virtual void Output(OutputHandler &OH) const
Definition: pipe.cc:1367
DynamicPipe::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: pipe.cc:1063
Node::iGetFirstColIndex
virtual integer iGetFirstColIndex(void) const
Definition: node.cc:88