MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
hutils.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/hydr/hutils.cc,v 1.36 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 "hutils.h"
43 
44 /* Accumulator - begin */
45 
46 Accumulator::Accumulator(unsigned int uL, const DofOwner* pDO,
47  HydraulicFluid* hf,
48  const PressureNode* p1,
49  doublereal St, doublereal start, doublereal As,
50  doublereal A_pipe, doublereal ms,
51  doublereal h_in, doublereal h_out,
52  doublereal P0, doublereal Pmax,doublereal k,
53  doublereal Wg, doublereal Kspr,doublereal F0,
54  doublereal cs, doublereal cv, doublereal ca,
55  flag fOut)
56 : Elem(uL, fOut),
57 HydraulicElem(uL, pDO, hf, fOut),
58 pNode1(p1),
59 stroke(St), start(start), area(As), area_pipe(A_pipe), mass(ms),
60 press0(P0), press_max(Pmax),
61 Kappa(k), weight(Wg), spring(Kspr),force0(F0), h_in(h_in), h_out(h_out),
62 c_spost(cs), c_vel(cv), c_acc(ca)
63 {
64  ASSERT(pNode1 != NULL);
65  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
66  ASSERT(St > std::numeric_limits<doublereal>::epsilon());
67  ASSERT(As > std::numeric_limits<doublereal>::epsilon());
68  ASSERT(A_pipe > std::numeric_limits<doublereal>::epsilon());
69  ASSERT(ms > std::numeric_limits<doublereal>::epsilon());
70  ASSERT(P0 >= 0.);
71  ASSERT(Pmax >= 0.);
72  ASSERT(k >= 0.);
73  ASSERT(Wg >= 0.);
74  ASSERT(Kspr >= 0.);
75  ASSERT(F0 >= 0.);
76 
77  s_min_gas = 0.;
78 
79  if (Kappa > 0. && press_max > 0.) {
80  s_min_gas = stroke*pow(press0/press_max, 1./Kappa);
81  }
82 
83  s_max = .999*stroke-s_min_gas; /* ho messo 0.999 perchè se uso un accumulatore senza gas avrei un termine che va a +infinito */
84  ratio2 = area*area/(area_pipe*area_pipe); /* rapporto (area accumulatore/area tubo)^2 */
85 }
86 
87 
88 Accumulator::~Accumulator(void)
89 {
90  NO_OP;
91 }
92 
93 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
94 HydraulicElem::Type
95 Accumulator::GetHydraulicType(void) const
96 {
97  return HydraulicElem::ACCUMULATOR;
98 }
99 
100 /* Contributo al file di restart */
101 std::ostream&
102 Accumulator::Restart(std::ostream& out) const
103 {
104  return out << " Accumulator not implemented yet!" << std::endl;
105 }
106 
107 
108 unsigned int
109 Accumulator::iGetNumDof(void) const
110 {
111  return 2;
112 }
113 
114 DofOrder::Order
115 Accumulator::GetDofType(unsigned int i) const
116 {
117  ASSERT(i >= 0 && i <= 1);
118  return DofOrder::DIFFERENTIAL;
119 }
120 
121 void
122 Accumulator::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
123 {
124  *piNumRows = 3;
125  *piNumCols = 3;
126 }
127 
128 VariableSubMatrixHandler&
129 Accumulator::AssJac(VariableSubMatrixHandler& WorkMat,
130  doublereal dCoef,
131  const VectorHandler& XCurr,
132  const VectorHandler& XPrimeCurr)
133 {
134  DEBUGCOUT("Entering Accumulator::AssJac()" << std::endl);
135 
136  FullSubMatrixHandler& WM = WorkMat.SetFull();
137  WM.Resize(3, 3);
138 
139  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
140  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
141  integer iFirstIndex = iGetFirstIndex();
142 
143  WM.PutRowIndex(1, iNode1RowIndex);
144  WM.PutColIndex(1, iNode1ColIndex);
145  WM.PutRowIndex(2, iFirstIndex+1);
146  WM.PutColIndex(2, iFirstIndex+1);
147  WM.PutRowIndex(3, iFirstIndex+2);
148  WM.PutColIndex(3, iFirstIndex+2);
149  doublereal density = HF->dGetDensity();
150 
151  /* unused? doublereal p1 = pNode1->dGetX(); */
152  s = XCurr(iFirstIndex+1); /* spostamento */
153  v = XCurr(iFirstIndex+2); /* velocita' */
154 
155  doublereal Jac11 = 0;
156  doublereal Jac12 = 0;
157  doublereal Jac13 = 0;
158  doublereal Jac21 = 0;
159  doublereal Jac22 = 0;
160  doublereal Jac23 = 0;
161  doublereal Jac31 = 0;
162  doublereal Jac32 = 0;
163  doublereal Jac33 = 0;
164 
165  if (s > s_max) {
166 #ifdef HYDR_DEVEL
167  DEBUGCOUT("AssJac(): ho superato la s_max: s" << s <<std::endl);
168 #endif
169  s = s_max;
170  } else if(s < 0.) {
171 #ifdef HYDR_DEVEL
172  DEBUGCOUT("AssJac(): sono negativo: s" << s <<std::endl);
173 #endif
174  s = 0.;
175  }
176 
177  Jac11 = 0.;
178  Jac12 = 0.;
179  Jac13 = -density*area*dCoef;
180 
181  Jac21 = area;
182  Jac22 = -press0*area*Kappa*dCoef*pow(stroke/(stroke-s), Kappa+1.)/stroke
183  -c1*dCoef-cf1*dCoef+dCoef*spring;
184  Jac23 = -mass-weight-h*density*area*ratio2*fabs(v)*dCoef
185  -c2*dCoef-c3-cf2*dCoef-cf3;
186 
187  Jac31 = 0.;
188  Jac32 = -1.;
189  Jac33 = dCoef;
190 
191 #ifdef HYDR_DEVEL
192  DEBUGCOUT("Jac11: " << Jac11 << std::endl);
193  DEBUGCOUT("Jac12: " << Jac12 << std::endl);
194  DEBUGCOUT("Jac13: " << Jac13 << std::endl);
195  DEBUGCOUT("Jac21: " << Jac21 << std::endl);
196  DEBUGCOUT("Jac22: " << Jac22 << std::endl);
197  DEBUGCOUT("Jac23: " << Jac23 << std::endl);
198  DEBUGCOUT("Jac31: " << Jac31 << std::endl);
199  DEBUGCOUT("Jac32: " << Jac32 << std::endl);
200  DEBUGCOUT("Jac33: " << Jac33 << std::endl);
201 #endif
202 
203  WM.PutCoef(1, 1, Jac11);
204  WM.PutCoef(1, 2, Jac12);
205  WM.PutCoef(1, 3, Jac13);
206  WM.PutCoef(2, 1, Jac21);
207  WM.PutCoef(2, 2, Jac22);
208  WM.PutCoef(2, 3, Jac23);
209  WM.PutCoef(3, 1, Jac31);
210  WM.PutCoef(3, 2, Jac32);
211  WM.PutCoef(3, 3, Jac33);
212 
213  return WorkMat;
214 }
215 
216 
217 SubVectorHandler&
218 Accumulator::AssRes(SubVectorHandler& WorkVec,
219  doublereal dCoef,
220  const VectorHandler& XCurr,
221  const VectorHandler& XPrimeCurr)
222 {
223  DEBUGCOUT("Entering Accumulator::AssRes()" << std::endl);
224 
225  WorkVec.Resize(3);
226 
227  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
228  integer iFirstIndex = iGetFirstIndex();
229 
230  doublereal p1 = pNode1->dGetX();
231  s = XCurr(iFirstIndex+1); /* spostamento */
232  v = XCurr(iFirstIndex+2); /* velocita' */
233  sp = XPrimeCurr(iFirstIndex+1); /* velocita' */
234  vp = XPrimeCurr(iFirstIndex+2); /* accelerazione */
235 
236  doublereal Res_1 = 0.;
237  doublereal Res_2 = 0.;
238  doublereal Res_3 = 0.;
239  pgas = 0.;
240  doublereal x0;
241  doublereal x0spring;
242  doublereal density = HF->dGetDensity();
243 
244  x0 = (stroke*(pow(p1/press0, 1./Kappa)-1.))/pow(p1/press0, 1./Kappa);
245  x0spring = (p1*area-force0)/spring;
246 
247  if (s < 0.) {
248  c1 = c_spost;
249  if (sp < 0.) {
250  c2 = c_vel; /* ho v negativa devo smorzare */
251  } else {
252  c2 = 0.;
253  }
254 
255  if (vp < 0.) {
256  c3 = c_acc; /* ho acc negativa devo smorzare */
257  } else {
258  c3 = 0.;
259  }
260 
261  c4 = 0.;
262  } else {
263  c1 = 0.;
264  c2 = 0.;
265  c3 = 0.;
266  c4 = 0.;
267  }
268 
269  if (s > s_max) {
270  cf1 = c_spost;
271 
272  if (sp > 0.) {
273  cf2 = c_vel; /* ho v positiva devo smorzare */
274  } else {
275  cf2 = 0.;
276  }
277 
278  if (vp > 0.) {
279  cf3 = c_acc; /* ho acc positiva devo smorzare */
280  } else {
281  cf3 = 0.;
282  }
283 
284  cf4 = 0.;
285  } else {
286  cf1 = 0.;
287  cf2 = 0.;
288  cf3 = 0.;
289  cf4 = 0.;
290  }
291 
292  if (sp > 0.) {
293  h = h_in; /* perdita concentrata fluido entra nell'accumulatore */
294  } else {
295  h = h_out; /* perdita concentrata fluido esce dall'accumulatore */
296  }
297 
298  if (s > s_max) {
299  pgas = press_max;
300 #ifdef HYDR_DEVEL
301  DEBUGCOUT("AssJac(): ho superato la s_max: s " << s <<std::endl);
302 #endif
303  } else if (s < 0.) {
304  pgas = press0;
305 #ifdef HYDR_DEVEL
306  DEBUGCOUT("AssJac(): sono negativo: s " << s <<std::endl);
307 #endif
308  } else {
309  pgas = press0*pow(stroke/(stroke-s), Kappa);
310  }
311 
312  Res_1 = density*v*area;
313  Res_2 = (mass+weight)*vp-p1*area+pgas*area+force0+spring*s
314  +copysign(h*.5*density*area*ratio2*pow(v, 2), v)+c1*s
315  +c2*sp+c3*vp+c4+cf1*(s-s_max)+cf2*sp+cf3*vp+cf4;
316  Res_3 = sp-v;
317 
318  flow = -Res_1; /* portata nodo 1 (per l'output) */
319 
320 #ifdef HYDR_DEVEL
321  DEBUGCOUT("x0: " << x0 << std::endl);
322  DEBUGCOUT("x0spring: " << x0spring << std::endl);
323  DEBUGCOUT("s: " << s << std::endl);
324  DEBUGCOUT("sp: " << sp << std::endl);
325  DEBUGCOUT("v: " << v << std::endl);
326  DEBUGCOUT("vp: " << vp << std::endl);
327  DEBUGCOUT("p1: " << p1 << std::endl);
328  DEBUGCOUT("pgas: " << pgas << std::endl);
329  DEBUGCOUT("smorzatore: "
330  << copysign(h*.5*density*area*ratio2*pow(v, 2), v) << std::endl);
331  DEBUGCOUT("stroke: " << stroke << std::endl);
332  DEBUGCOUT("area: " << area << std::endl);
333  DEBUGCOUT("area_pipe: " << area_pipe << std::endl);
334  DEBUGCOUT("mass: " << mass << std::endl);
335  DEBUGCOUT("press0: " << press0 << std::endl);
336  DEBUGCOUT("press_max: " << press_max << std::endl);
337  DEBUGCOUT("Kappa: " << Kappa << std::endl);
338  DEBUGCOUT("weight: " << weight << std::endl);
339  DEBUGCOUT("spring: " << spring << std::endl);
340  DEBUGCOUT("force0: " << force0 << std::endl);
341  DEBUGCOUT("s_min_gas: " << s_min_gas << std::endl);
342  DEBUGCOUT("s_max: " << s_max << std::endl);
343  DEBUGCOUT("c1: " << c1 << std::endl);
344  DEBUGCOUT("c2: " << c2 << std::endl);
345  DEBUGCOUT("c3: " << c3 << std::endl);
346  DEBUGCOUT("c4: " << c4 << std::endl);
347  DEBUGCOUT("cf1: " << cf1 << std::endl);
348  DEBUGCOUT("cf2: " << cf2 << std::endl);
349  DEBUGCOUT("cf3: " << cf3 << std::endl);
350  DEBUGCOUT("cf4: " << cf4 << std::endl);
351  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
352  DEBUGCOUT("Res_2: " << Res_2 << std::endl);
353  DEBUGCOUT("Res_3: " << Res_3 << std::endl);
354 #endif
355 
356  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
357  WorkVec.PutItem(2, iFirstIndex+1, Res_2);
358  WorkVec.PutItem(3, iFirstIndex+2, Res_3);
359 
360  return WorkVec;
361 }
362 
363 
364 void
365 Accumulator::Output(OutputHandler& OH) const
366 {
367  if (bToBeOutput()) {
368  std::ostream& out = OH.Hydraulic();
369  out << std::setw(8) << GetLabel()
370  << " " << s << " " << v << " " << vp << " " << pgas
371  << " " << flow << " " << density << std::endl;
372  }
373 }
374 
375 void
376 Accumulator::SetValue(DataManager *pDM,
377  VectorHandler& X, VectorHandler& XP,
378  SimulationEntity::Hints *ph)
379 {
380  integer i = iGetFirstIndex();
381 
382  X.PutCoef(i+1, start);
383  X.PutCoef(i+2, 0.);
384  XP.PutCoef(i+1, 0.);
385  XP.PutCoef(i+2, 0.);
386 }
387 
388 /* Accumulator - end */
389 
390 
391 /* Tank - begin */
392 
393 Tank::Tank(unsigned int uL, const DofOwner* pDO, HydraulicFluid* hf,
394  const PressureNode* p1, const PressureNode* p2,
395  doublereal Ps,doublereal A_pipe,
396  doublereal A_serb, doublereal lev, doublereal s_mx,
397  doublereal s_mn, doublereal c_s, flag fOut)
398 : Elem(uL, fOut),
399 HydraulicElem(uL, pDO, hf, fOut),
400 pNode1(p1), pNode2(p2),
401 press(Ps), area_pipe(A_pipe), area_serb(A_serb),level(lev),
402 s_max(s_mx), s_min(s_mn), c_spost(c_s)
403 {
404  ASSERT(pNode1 != NULL);
405  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
406  ASSERT(pNode2 != NULL);
407  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
408  ASSERT(Ps > std::numeric_limits<doublereal>::epsilon());
409  ASSERT(A_pipe > std::numeric_limits<doublereal>::epsilon());
410  ASSERT(A_serb > std::numeric_limits<doublereal>::epsilon());
411  ASSERT(lev >= 0.);
412  ASSERT(s_mx >= 0.);
413  ASSERT(s_mn >= 0.);
414 
415  Kappa1 = 1.;
416  Kappa2 = .5;
417 }
418 
419 
420 Tank::~Tank(void)
421 {
422  NO_OP;
423 }
424 
425 
426 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
427 HydraulicElem::Type
428 Tank::GetHydraulicType(void) const
429 {
430  return HydraulicElem::TANK;
431 }
432 
433 
434 /* Contributo al file di restart */
435 std::ostream&
436 Tank::Restart(std::ostream& out) const
437 {
438  return out << "Tank not implemented yet!" << std::endl;
439 }
440 
441 
442 unsigned int
443 Tank::iGetNumDof(void) const
444 {
445  return 1;
446 }
447 
448 
449 DofOrder::Order
450 Tank::GetDofType(unsigned int i) const
451 {
452 
453  ASSERT(i >= 0 && i <= 1);
454  return DofOrder::DIFFERENTIAL;
455 }
456 
457 
458 void
459 Tank::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const {
460  *piNumRows = 3;
461  *piNumCols = 3;
462 }
463 
464 
465 VariableSubMatrixHandler&
466 Tank::AssJac(VariableSubMatrixHandler& WorkMat,
467  doublereal dCoef,
468  const VectorHandler& XCurr,
469  const VectorHandler& XPrimeCurr)
470 {
471  DEBUGCOUT("Entering Tank::AssJac()" << std::endl);
472 
473  FullSubMatrixHandler& WM = WorkMat.SetFull();
474  WM.Resize(3, 3);
475 
476  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
477  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
478  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
479  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
480  integer iFirstIndex = iGetFirstIndex();
481 
482  WM.PutRowIndex(1, iNode1RowIndex);
483  WM.PutRowIndex(2, iNode2RowIndex);
484  WM.PutColIndex(1, iNode1ColIndex);
485  WM.PutColIndex(2, iNode2ColIndex);
486 
487  WM.PutRowIndex(3, iFirstIndex+1);
488  WM.PutColIndex(3, iFirstIndex+1);
489 
490  doublereal p1 = pNode1->dGetX();
491  doublereal p2 = pNode2->dGetX();
492  /* p1 = XCurr(pNode1->iGetFirstRowIndex()+1); */
493 
494  s = XCurr(iFirstIndex+1); /* livello */
495 
496  /* salto di pressione nodo1-Tank */
497  doublereal jumpPres1S = fabs(p1-press);
498  /* evito di dividere per un numero troppo piccolo */
499  if (jumpPres1S < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
500  jumpPres1S = 1.e8*std::numeric_limits<doublereal>::epsilon();
501  }
502 
503  /* salto di pressione Tank-nodo2 */
504  doublereal jumpPresS2 = fabs(press-p2);
505  /* evito di dividere per un numero troppo piccolo */
506  if (jumpPresS2 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
507  jumpPresS2 = 1.e8*std::numeric_limits<doublereal>::epsilon();
508  }
509  doublereal density = HF->dGetDensity((p1+p2)/2.);
510 
511  doublereal Jac11 = -.5*density*area_pipe*sqrt(2./(Kappa1*density*jumpPres1S));
512  doublereal Jac12 = 0.;
513  doublereal Jac13 = 0.;
514 
515  doublereal Jac21 = 0.;
516  doublereal Jac22 = -.5*density*area_pipe*sqrt(2./(Kappa2*density*jumpPresS2));
517  doublereal Jac23 = 0.;
518 
519  doublereal Jac31 = .5*area_pipe/area_serb*sqrt(2./(Kappa1*density*jumpPres1S));
520  doublereal Jac32 = .5*area_pipe/area_serb*sqrt(2./(Kappa2*density*jumpPresS2));
521  doublereal Jac33 = -1.-c_spost*dCoef;
522 
523  if(s < s_min) {
524  /* Livello dell'olio sotto la soglia minima:
525  * esco dalla presa di emergenza */
526 #ifdef HYDR_DEVEL
527  DEBUGCOUT("Jac Esco dalla presa di emergenza: " << std::endl);
528 #endif
529  }
530  if (s > s_max && -flow1 > flow2) {
531  /* se e' pieno non puo' entrare di piu' di quella che esce */
532  Jac31 = Jac32 = 0.;
533 #ifdef HYDR_DEVEL
534  DEBUGCOUT("Jac Serbatoio pieno : " << std::endl);
535 #endif
536  }
537 
538  if (s < 0. && flow2 > -flow1) {
539  /* se e' vuoto non puo' uscire di piu' di quella che entra */
540  Jac31 = Jac32 = 0.;
541 #ifdef HYDR_DEVEL
542  DEBUGCOUT("Jac Serbatoio vuoto : "<< std::endl);
543 #endif
544  }
545 
546  WM.PutCoef(1, 1, Jac11);
547  WM.PutCoef(1, 2, Jac12);
548  WM.PutCoef(1, 3, Jac13);
549  WM.PutCoef(2, 1, Jac21);
550  WM.PutCoef(2, 2, Jac22);
551  WM.PutCoef(2, 3, Jac23);
552  WM.PutCoef(3, 1, Jac31);
553  WM.PutCoef(3, 2, Jac32);
554  WM.PutCoef(3, 3, Jac33);
555 
556  return WorkMat;
557 }
558 
559 
560 SubVectorHandler&
561 Tank::AssRes(SubVectorHandler& WorkVec,
562  doublereal dCoef,
563  const VectorHandler& XCurr,
564  const VectorHandler& XPrimeCurr)
565 {
566  DEBUGCOUT("Entering Tank::AssRes()" << std::endl);
567 
568  WorkVec.Resize(3);
569 
570  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
571  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
572 
573  doublereal p1 = pNode1->dGetX();
574  doublereal p2 = pNode2->dGetX();
575 
576  integer iFirstIndex = iGetFirstIndex();
577 
578  s = XCurr(iFirstIndex+1); /* livello */
579  sp = XPrimeCurr(iFirstIndex+1); /* velocita' del livello */
580 
581  if (s < 0.) {
582  silent_cerr("Tank(" << GetLabel() << ": negative fluid level "
583  << s << " impossible" << std::endl);
584  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
585  }
586 
587  doublereal jumpPres1S = fabs(p1-press);
588  doublereal jumpPresS2 = fabs(press-p2);
589  doublereal Res_1 = 0.;
590  doublereal Res_2 = 0.;
591  doublereal Res_3 = 0.;
592  /* unused? doublereal Res_4 = 0.; */
593  c_spost = 0.;
594 
595  doublereal density = HF->dGetDensity((p1+p2/2.));
596 
597  Res_1 = density*area_pipe*sqrt(2./(Kappa1*density))*copysign(sqrt(jumpPres1S), p1-press);
598  Res_2 = -density*area_pipe*sqrt(2./(Kappa2*density))*copysign(sqrt(jumpPresS2), press-p2);
599 
600  if (s > s_max && Res_1 > -Res_2) {
601  /* se e' pieno non puo' entrare di piu' di quella che esce */
602  Res_1 = -Res_2;
603  c_spost = 1.;
604 #ifdef HYDR_DEVEL
605  DEBUGCOUT("Res Serbatoio pieno: " << std::endl);
606 #endif
607  }
608 
609  if (s <s_min) {
610  /* livello dell'olio sotto al livello minimo:
611  * esco dalla presa di emergenza */
612 #ifdef HYDR_DEVEL
613  DEBUGCOUT("Esco dalla presa di emergenza: " << std::endl);
614 #endif
615  }
616 
617 
618  if (s < 0. && Res_2 < -Res_1) {
619  c_spost = 0.;
620  /* se e' vuoto non puo' uscire di piu' di quella che entra */
621  Res_2 = -Res_1;
622 #ifdef HYDR_DEVEL
623  DEBUGCOUT("Res Serbatoio vuoto: " << std::endl);
624 #endif
625  }
626 
627  Res_3 = sp+(-Res_1-Res_2)/(area_serb*density)+c_spost*(s-s_max);
628 
629  flow1 = -Res_1; /* per l'output */
630  flow2 = -Res_2; /* per l'output */
631 
632 #ifdef HYDR_DEVEL
633  DEBUGCOUT("Kappa1: " << Kappa1 << std::endl);
634  DEBUGCOUT("Kappa2: " << Kappa2 << std::endl);
635  DEBUGCOUT("jumpPres1S: " << jumpPres1S << std::endl);
636  DEBUGCOUT("jumpPresS2: " << jumpPresS2 << std::endl);
637  DEBUGCOUT("density: " << density << std::endl);
638  DEBUGCOUT("p1: " << p1 << std::endl);
639  DEBUGCOUT("p2: " << p2 << std::endl);
640  DEBUGCOUT("press: " << press << std::endl);
641  DEBUGCOUT("s: " << s << std::endl);
642  DEBUGCOUT("level: " << level << std::endl);
643  DEBUGCOUT("s_max: " << s_max << std::endl);
644  DEBUGCOUT("s_min: " << s_min << std::endl);
645  DEBUGCOUT("sp: " << sp << std::endl);
646  DEBUGCOUT("area_pipe: " << area_pipe << std::endl);
647  DEBUGCOUT("area_serb: " << area_serb << std::endl);
648  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)" << std::endl);
649  DEBUGCOUT("-Res_1: " << -Res_1 << " (portata nodo1) " << std::endl);
650  DEBUGCOUT("-Res_2: " << -Res_2 << " (portata nodo2) " << std::endl);
651  DEBUGCOUT("Res_3: " << Res_3 << std::endl);
652 #endif
653 
654  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
655  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
656  WorkVec.PutItem(3, iFirstIndex+1, Res_3);
657 
658  return WorkVec;
659 }
660 
661 void
662 Tank::Output(OutputHandler& OH) const
663 {
664  if (bToBeOutput()) {
665  std::ostream& out = OH.Hydraulic();
666  out << std::setw(8) << GetLabel();
667  out << " " << s << " " << sp << " " << flow1 << " " << flow2 << std::endl;
668  }
669 }
670 
671 void
672 Tank::SetValue(DataManager *pDM,
673  VectorHandler& X , VectorHandler& /* XP */ ,
674  SimulationEntity::Hints *ph)
675 {
676  integer i = iGetFirstIndex();
677  X.PutCoef(i+1, level);
678 }
679 
680 /* Tank - end */
Tank::flow2
doublereal flow2
Definition: hutils.h:155
Tank::area_serb
doublereal area_serb
Definition: hutils.h:145
Accumulator::h_in
doublereal h_in
Definition: hutils.h:58
Accumulator::cf1
doublereal cf1
Definition: hutils.h:74
Tank::press
doublereal press
Definition: hutils.h:143
Tank::pNode2
const PressureNode * pNode2
Definition: hutils.h:142
Accumulator::s
doublereal s
Definition: hutils.h:79
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Accumulator::weight
doublereal weight
Definition: hutils.h:55
flag
long int flag
Definition: mbdyn.h:43
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
Tank::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: hutils.cc:466
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Accumulator::press0
doublereal press0
Definition: hutils.h:52
Accumulator::press_max
doublereal press_max
Definition: hutils.h:53
Accumulator::~Accumulator
~Accumulator(void)
Definition: hutils.cc:88
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
Accumulator::c_acc
doublereal c_acc
Definition: hutils.h:69
Tank::s
doublereal s
Definition: hutils.h:152
Accumulator::spring
doublereal spring
Definition: hutils.h:56
Accumulator::c2
doublereal c2
Definition: hutils.h:71
Accumulator::stroke
doublereal stroke
Definition: hutils.h:47
FullSubMatrixHandler::Resize
void Resize(integer iNewRow, integer iNewCol)
Definition: submat.cc:138
Tank::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: hutils.cc:443
Tank::level
doublereal level
Definition: hutils.h:146
FullSubMatrixHandler::PutCoef
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
Tank::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: hutils.cc:561
Tank::Kappa2
doublereal Kappa2
Definition: hutils.h:150
NO_OP
#define NO_OP
Definition: myassert.h:74
Accumulator::s_max
doublereal s_max
Definition: hutils.h:65
Accumulator::Output
virtual void Output(OutputHandler &OH) const
Definition: hutils.cc:365
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
Accumulator::h
doublereal h
Definition: hutils.h:60
Accumulator::area
doublereal area
Definition: hutils.h:49
Tank::flow1
doublereal flow1
Definition: hutils.h:154
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
Accumulator::h_out
doublereal h_out
Definition: hutils.h:59
Accumulator::force0
doublereal force0
Definition: hutils.h:57
Tank::~Tank
~Tank(void)
Definition: hutils.cc:420
Accumulator::c3
doublereal c3
Definition: hutils.h:72
HydraulicFluid::dGetDensity
virtual doublereal dGetDensity(void) const =0
Tank::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: hutils.cc:436
ScalarAlgebraicNode::dGetX
virtual const doublereal & dGetX(void) const
Definition: node.h:492
Accumulator::v
doublereal v
Definition: hutils.h:80
Accumulator::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: hutils.cc:102
Accumulator::pgas
doublereal pgas
Definition: hutils.h:83
HydraulicElem::HF
HydraulicFluid * HF
Definition: preselem.h:79
Accumulator::c4
doublereal c4
Definition: hutils.h:73
Accumulator::vp
doublereal vp
Definition: hutils.h:82
Tank::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: hutils.cc:672
copysign
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
Tank::sp
doublereal sp
Definition: hutils.h:153
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
Accumulator::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: hutils.cc:129
Accumulator::sp
doublereal sp
Definition: hutils.h:81
Node::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: node.cc:82
Accumulator::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: hutils.cc:115
Accumulator::c_spost
doublereal c_spost
Definition: hutils.h:67
Accumulator::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: hutils.cc:95
Accumulator::s_min_gas
doublereal s_min_gas
Definition: hutils.h:64
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
Tank::c_spost
doublereal c_spost
Definition: hutils.h:151
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
Accumulator::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: hutils.cc:376
Tank::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: hutils.cc:459
Accumulator::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: hutils.cc:109
Accumulator::start
doublereal start
Definition: hutils.h:48
Tank::pNode1
const PressureNode * pNode1
Definition: hutils.h:141
Accumulator::Accumulator
Accumulator(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, doublereal St, doublereal start, doublereal As, doublereal A_pipe, doublereal ms, doublereal h_in, doublereal h_out, doublereal P0, doublereal Pmax, doublereal k, doublereal Wg, doublereal Kspr, doublereal force0, doublereal cs, doublereal cv, doublereal ca, flag fOut)
Definition: hutils.cc:46
Tank::Output
virtual void Output(OutputHandler &OH) const
Definition: hutils.cc:662
Elem
Definition: elem.h:75
Accumulator::cf2
doublereal cf2
Definition: hutils.h:75
Accumulator::mass
doublereal mass
Definition: hutils.h:51
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
Tank::s_max
doublereal s_max
Definition: hutils.h:147
Tank::Tank
Tank(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, doublereal press, doublereal A_pipe, doublereal A_serb, doublereal lev, doublereal s_mx, doublereal s_mn, doublereal c_s, flag fOut)
Definition: hutils.cc:393
Tank::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: hutils.cc:428
PressureNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: presnode.h:54
Accumulator::density
doublereal density
Definition: hutils.h:61
Accumulator::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: hutils.cc:122
Tank::s_min
doublereal s_min
Definition: hutils.h:148
OutputHandler::Hydraulic
std::ostream & Hydraulic(void) const
Definition: output.h:492
Accumulator::cf3
doublereal cf3
Definition: hutils.h:76
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
Accumulator::pNode1
const PressureNode * pNode1
Definition: hutils.h:46
Accumulator::Kappa
doublereal Kappa
Definition: hutils.h:54
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
Accumulator::flow
doublereal flow
Definition: hutils.h:84
Tank::Kappa1
doublereal Kappa1
Definition: hutils.h:149
Accumulator::ratio2
doublereal ratio2
Definition: hutils.h:66
Accumulator::c1
doublereal c1
Definition: hutils.h:70
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Accumulator::area_pipe
doublereal area_pipe
Definition: hutils.h:50
Accumulator::cf4
doublereal cf4
Definition: hutils.h:77
Accumulator::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: hutils.cc:218
Accumulator::c_vel
doublereal c_vel
Definition: hutils.h:68
VectorHandler::Resize
virtual void Resize(integer iNewSize)=0
Tank::area_pipe
doublereal area_pipe
Definition: hutils.h:144
Tank::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: hutils.cc:450
Node::iGetFirstColIndex
virtual integer iGetFirstColIndex(void) const
Definition: node.cc:88