MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
valve.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/hydr/valve.cc,v 1.45 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 "valve.h"
43 
44 /* Control_valve - begin */
45 
46 Control_valve::Control_valve(unsigned int uL, const DofOwner* pDO,
47  HydraulicFluid* hf,
48  const PressureNode* p1, const PressureNode* p2,
49  const PressureNode* p3, const PressureNode* p4,
50  doublereal A_max, doublereal Loss_A, const DriveCaller* pDC,
51  flag fOut)
52 : Elem(uL, fOut),
53 HydraulicElem(uL, pDO, hf, fOut),
54 DriveOwner(pDC),
55 pNode1(p1), pNode2(p2), pNode3(p3), pNode4(p4),
56 area_max(A_max), loss_area(Loss_A),
57 A1min(2.*area_max*loss_area),
58 A2min(area_max*loss_area),
59 A3min(2.*area_max*loss_area),
60 A4min(area_max*loss_area)
61 {
62  ASSERT(pNode1 != NULL);
63  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
64  ASSERT(pNode2 != NULL);
65  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
66  ASSERT(pNode3 != NULL);
67  ASSERT(pNode3->GetNodeType() == Node::HYDRAULIC);
68  ASSERT(pNode4 != NULL);
69  ASSERT(pNode4->GetNodeType() == Node::HYDRAULIC);
70  ASSERT(A_max > std::numeric_limits<doublereal>::epsilon());
71  ASSERT(loss_area >= 0.);
72 
73  /*
74  * Cd = pi / ( pi + 2 ) ~= .611
75  *
76  * cfr. Merritt, Hydraulic Control Systems, Par. 3.4, pp. 39-45
77  * John Wiley & Sons, New York, 1967
78  */
79  Cd = M_PI / ( M_PI + 2. );
80 }
81 
82 Control_valve::~Control_valve(void)
83 {
84  NO_OP;
85 }
86 
87 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
88 HydraulicElem::Type Control_valve::GetHydraulicType(void) const
89 {
90  return HydraulicElem::CONTROL_VALVE;
91 }
92 
93 /* Contributo al file di restart */
94 std::ostream& Control_valve::Restart(std::ostream& out) const
95 {
96  return out << "Control_valve not implemented yet!" << std::endl;
97 }
98 
99 unsigned int Control_valve::iGetNumDof(void) const
100 {
101  return 0;
102 }
103 
104 DofOrder::Order Control_valve::GetDofType(unsigned int i) const
105 {
106  silent_cerr("ControlValve(" << GetLabel() << ") has no dofs!" << std::endl);
107  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
108 }
109 
110 
111 void
112 Control_valve::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
113 {
114  *piNumRows = 4;
115  *piNumCols = 4;
116 }
117 
118 
119 VariableSubMatrixHandler&
120 Control_valve::AssJac(VariableSubMatrixHandler& WorkMat,
121  doublereal dCoef,
122  const VectorHandler& XCurr,
123  const VectorHandler& XPrimeCurr)
124 {
125  DEBUGCOUT("Entering Control_valve::AssJac()" << std::endl);
126 
127  FullSubMatrixHandler& WM = WorkMat.SetFull();
128  WM.Resize(4, 4);
129 
130  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
131  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
132  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
133  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
134  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
135  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
136  integer iNode3ColIndex = pNode3->iGetFirstColIndex()+1;
137  integer iNode4ColIndex = pNode4->iGetFirstColIndex()+1;
138 
139  WM.PutRowIndex(1, iNode1RowIndex);
140  WM.PutRowIndex(2, iNode2RowIndex);
141  WM.PutRowIndex(3, iNode3RowIndex);
142  WM.PutRowIndex(4, iNode4RowIndex);
143  WM.PutColIndex(1, iNode1ColIndex);
144  WM.PutColIndex(2, iNode2ColIndex);
145  WM.PutColIndex(3, iNode3ColIndex);
146  WM.PutColIndex(4, iNode4ColIndex);
147 
148  doublereal p1 = pNode1->dGetX();
149  doublereal p2 = pNode2->dGetX();
150  doublereal p3 = pNode3->dGetX();
151  doublereal p4 = pNode4->dGetX();
152  doublereal density = HF->dGetDensity();
153 
154  doublereal jumpPres12 = fabs(p1-p2); /* salto di pressione nodo1 & nodo3 */
155  doublereal jumpPres13 = fabs(p1-p3); /* salto di pressione nodo1 & nodo4 */
156  doublereal jumpPres24 = fabs(p2-p4); /* salto di pressione nodo2 & nodo3 */
157  doublereal jumpPres34 = fabs(p3-p4); /* salto di pressione nodo2 & nodo4 */
158 
159  if (jumpPres12 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
160  jumpPres12 = 1.e8*std::numeric_limits<doublereal>::epsilon();
161  }
162  if (jumpPres13 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
163  jumpPres13 = 1.e8*std::numeric_limits<doublereal>::epsilon();
164  }
165  if (jumpPres24 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
166  jumpPres24 = 1.e8*std::numeric_limits<doublereal>::epsilon();
167  }
168  if (jumpPres34 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
169  jumpPres34 = 1.e8*std::numeric_limits<doublereal>::epsilon();
170  }
171 
172  doublereal primo = Cd*A1/sqrt(2*jumpPres12/density);
173  doublereal secondo = Cd*A2/sqrt(2*jumpPres13/density);
174  doublereal terzo = Cd*A4/sqrt(2*jumpPres24/density);
175  doublereal quarto = Cd*A3/sqrt(2*jumpPres34/density);
176 
177  doublereal Jac11 = -primo-secondo;
178  doublereal Jac12 = primo;
179  doublereal Jac13 = secondo;
180  // doublereal Jac14 = 0.;
181  doublereal Jac21 = primo;
182  doublereal Jac22 = -primo-terzo;
183  // doublereal Jac23 = 0.;
184  doublereal Jac24 = terzo;
185  doublereal Jac31 = secondo;
186  // doublereal Jac32 = 0.;
187  doublereal Jac33 = -secondo-quarto;
188  doublereal Jac34 = quarto;
189  // doublereal Jac41 = 0.;
190  doublereal Jac42 = terzo;
191  doublereal Jac43 = quarto;
192  doublereal Jac44 = -terzo-quarto;
193 
194  WM.PutCoef(1, 1, Jac11);
195  WM.PutCoef(1, 2, Jac12);
196  WM.PutCoef(1, 3, Jac13);
197  // WM.PutCoef(1, 4, Jac14);
198  WM.PutCoef(2, 1, Jac21);
199  WM.PutCoef(2, 2, Jac22);
200  // WM.PutCoef(2, 3, Jac23);
201  WM.PutCoef(2, 4, Jac24);
202  WM.PutCoef(3, 1, Jac31);
203  // WM.PutCoef(3, 2, Jac32);
204  WM.PutCoef(3, 3, Jac33);
205  WM.PutCoef(3, 4, Jac24);
206  // WM.PutCoef(4, 1, Jac41);
207  WM.PutCoef(4, 2, Jac42);
208  WM.PutCoef(4, 3, Jac43);
209  WM.PutCoef(4, 4, Jac44);
210 
211  return WorkMat;
212 }
213 
214 
215 SubVectorHandler&
216 Control_valve::AssRes(SubVectorHandler& WorkVec,
217  doublereal dCoef,
218  const VectorHandler& XCurr,
219  const VectorHandler& XPrimeCurr)
220 {
221  DEBUGCOUT("Entering Control_valve::AssRes()" << std::endl);
222 
223  WorkVec.Resize(4);
224 
225  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
226  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
227  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
228  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
229 
230  doublereal p1 = pNode1->dGetX();
231  doublereal p2 = pNode2->dGetX();
232  doublereal p3 = pNode3->dGetX();
233  doublereal p4 = pNode4->dGetX();
234  doublereal density = HF->dGetDensity();
235 
236  Stato = pGetDriveCaller()->dGet();
237 
238  doublereal dp12 = p1-p2;
239  doublereal dp13 = p1-p3;
240  doublereal dp24 = p2-p4;
241  doublereal dp34 = p3-p4;
242 
243  doublereal jumpPres12 = fabs(dp12); /* salto di pressione nodo1 & nodo2 */
244  doublereal jumpPres13 = fabs(dp13); /* salto di pressione nodo1 & nodo3 */
245  doublereal jumpPres24 = fabs(dp24); /* salto di pressione nodo2 & nodo4 */
246  doublereal jumpPres34 = fabs(dp34); /* salto di pressione nodo3 & nodo4 */
247 
248  if (Stato > 1.) {
249  Stato = 1.;
250  } else if (Stato < -1.) {
251  Stato = -1.;
252  }
253 
254  if (Stato > 0.) {
255  A1 = Stato*area_max+A1min;
256  A2 = A2min;
257  A3 = Stato*area_max+A3min;
258  A4 = A4min;
259  } else {
260  A1 = A1min;
261  A2 = -Stato*area_max+A2min;
262  A3 = A3min;
263  A4 = -Stato*area_max+A4min;
264  }
265 
266  doublereal Q12 = Cd*A1*copysign(sqrt(2.*jumpPres12*density), dp12);
267  doublereal Q13 = Cd*A2*copysign(sqrt(2.*jumpPres13*density), dp13);
268  doublereal Q24 = Cd*A4*copysign(sqrt(2.*jumpPres24*density), dp24);
269  doublereal Q34 = Cd*A3*copysign(sqrt(2.*jumpPres34*density), dp34);
270 
271  doublereal Res_1 = Q12+Q13;
272  doublereal Res_2 = -Q12+Q24;
273  doublereal Res_3 = -Q13+Q34;
274  doublereal Res_4 = -Q34-Q24;
275 
276  flow1 = -Res_1;
277  flow2 = -Res_2;
278  flow3 = -Res_3;
279  flow4 = -Res_4;
280 
281 #ifdef HYDR_DEVEL
282  DEBUGCOUT("Stato: " << Stato << std::endl);
283  DEBUGCOUT("A1: " << A1 << std::endl);
284  DEBUGCOUT("A2: " << A2 << std::endl);
285  DEBUGCOUT("A3: " << A3 << std::endl);
286  DEBUGCOUT("A4: " << A4 << std::endl);
287  DEBUGCOUT("p1: " << p1 << std::endl);
288  DEBUGCOUT("p2: " << p2 << std::endl);
289  DEBUGCOUT("p3: " << p3 << std::endl);
290  DEBUGCOUT("p4: " << p4 << std::endl);
291  DEBUGCOUT("Cd: " << Cd << std::endl);
292  DEBUGCOUT("density: " << density << std::endl);
293  DEBUGCOUT("Area_max: " << area_max << std::endl);
294  DEBUGCOUT("Loss_area in %: " << loss_area << std::endl);
295  DEBUGCOUT("Q12: " << Q12 << std::endl);
296  DEBUGCOUT("Q13: " << Q13 << std::endl);
297  DEBUGCOUT("Q24: " << Q24 << std::endl);
298  DEBUGCOUT("Q34: " << Q34 << std::endl);
299  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)"<< std::endl);
300  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
301  DEBUGCOUT("-Res_2 (portata nodo2): " << -Res_2 << std::endl);
302  DEBUGCOUT("-Res_3 (portata nodo3): " << -Res_3 << std::endl);
303  DEBUGCOUT("-Res_4 (portata nodo4): " << -Res_4 << std::endl);
304 #endif
305 
306  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
307  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
308  WorkVec.PutItem(3, iNode3RowIndex, Res_3);
309  WorkVec.PutItem(4, iNode4RowIndex, Res_4);
310 
311  return WorkVec;
312 }
313 
314 void Control_valve::Output(OutputHandler& OH) const
315 {
316  if (bToBeOutput()) {
317  std::ostream& out = OH.Hydraulic();
318  out << std::setw(8) << GetLabel()
319  << " " << Stato
320  << " " << flow1 << " " << flow2
321  << " " << flow3 << " " << flow4 << std::endl;
322  }
323 }
324 
325 /* Control_valve - end */
326 
327 
328 /* Control_valve2 - begin */
329 
330 Control_valve2::Control_valve2(unsigned int uL, const DofOwner* pDO,
331  HydraulicFluid* hf,
332  const PressureNode* p1, const PressureNode* p2,
333  const PressureNode* p3, const PressureNode* p4,
334  doublereal A_max, doublereal Loss_A,
335  const DriveCaller* pDC,
336  flag fOut)
337 : Elem(uL, fOut),
338 HydraulicElem(uL, pDO, hf, fOut),
339 DriveOwner(pDC),
340 area_max(A_max), loss_area(Loss_A), area_min(area_max*loss_area)
341 {
342  pNode[N1] = p1;
343  pNode[N2] = p2;
344  pNode[N3] = p3;
345  pNode[N4] = p4;
346 
347  ASSERT(pNode[N1] != NULL);
348  ASSERT(pNode[N1]->GetNodeType() == Node::HYDRAULIC);
349  ASSERT(pNode[N2] != NULL);
350  ASSERT(pNode[N2]->GetNodeType() == Node::HYDRAULIC);
351  ASSERT(pNode[N3] != NULL);
352  ASSERT(pNode[N3]->GetNodeType() == Node::HYDRAULIC);
353  ASSERT(pNode[N4] != NULL);
354  ASSERT(pNode[N4]->GetNodeType() == Node::HYDRAULIC);
355 
356  ASSERT(area_max > std::numeric_limits<doublereal>::epsilon());
357  ASSERT(loss_area > 1.e-9); /*
358  * se = 0. occorre fare un elemento
359  * apposta con solo 2 dof
360  */
361 
362  Cd = .611; /* coefficiente di perdita */
363 }
364 
365 Control_valve2::~Control_valve2(void)
366 {
367  NO_OP;
368 }
369 
370 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
371 HydraulicElem::Type
372 Control_valve2::GetHydraulicType(void) const
373 {
374  return HydraulicElem::CONTROL_VALVE;
375 }
376 
377 /* Contributo al file di restart */
378 std::ostream&
379 Control_valve2::Restart(std::ostream& out) const
380 {
381  return out << "Control_valve2 not implemented yet!" << std::endl;
382 }
383 
384 unsigned int
385 Control_valve2::iGetNumDof(void) const
386 {
387  return LAST_Q;
388 }
389 
390 DofOrder::Order Control_valve2::GetDofType(unsigned int i) const
391 {
392  ASSERT(i >= 0 && i < iGetNumDof());
393  return DofOrder::ALGEBRAIC;
394 }
395 
396 void
397 Control_valve2::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
398 {
399  *piNumRows = 4+LAST_Q;
400  *piNumCols = 4+LAST_Q;
401 }
402 
403 VariableSubMatrixHandler&
404 Control_valve2::AssJac(VariableSubMatrixHandler& WorkMat,
405  doublereal dCoef,
406  const VectorHandler& XCurr,
407  const VectorHandler& XPrimeCurr)
408 {
409  DEBUGCOUT("Entering Control_valve2::AssJac()" << std::endl);
410 
411  FullSubMatrixHandler& WM = WorkMat.SetFull();
412  integer iNumRows, iNumCols;
413  WorkSpaceDim(&iNumRows, &iNumCols);
414  WM.ResizeReset(iNumRows, iNumCols);
415 
416  integer iFirstIndex = iGetFirstIndex();
417 
418  for (int i = 0; i < LAST_N; i++) {
419  WM.PutRowIndex(1+i, pNode[i]->iGetFirstRowIndex()+1);
420  WM.PutColIndex(1+i, pNode[i]->iGetFirstColIndex()+1);
421  }
422 
423  for (int i = 1; i <= LAST_Q; i++) {
424  WM.PutRowIndex(4+i, iFirstIndex+i);
425  WM.PutColIndex(4+i, iFirstIndex+i);
426  }
427 
428  doublereal dKappa = 2.*Cd*Cd*density;
429 
430  /* Q12 */
431  WM.PutCoef(1, 5, -1.);
432  WM.PutCoef(2, 5, 1.);
433 
434  doublereal a = A[0]*A[0];
435  WM.PutCoef(5, 1, -a);
436  WM.PutCoef(5, 2, a);
437 
438  /* Q34 */
439  WM.PutCoef(3, 6, -1.);
440  WM.PutCoef(4, 6, 1.);
441 
442  a = A[1]*A[1];
443  WM.PutCoef(6, 3, -a);
444  WM.PutCoef(6, 4, a);
445 
446  /* Q13 */
447  WM.PutCoef(1, 7, -1.);
448  WM.PutCoef(3, 7, 1.);
449 
450  a = A[2]*A[2];
451  WM.PutCoef(7, 1, -a);
452  WM.PutCoef(7, 3, a);
453 
454  /* Q24 */
455  WM.PutCoef(2, 8, -1.);
456  WM.PutCoef(4, 8, 1.);
457 
458  a = A[3]*A[3];
459  WM.PutCoef(8, 2, -a);
460  WM.PutCoef(8, 4, a);
461 
462 #ifdef VALVE_6
463  /* Q14 */
464  WM.PutCoef(1, 9, -1.);
465  WM.PutCoef(4, 9, 1.);
466 
467  a = A[4]*A[4];
468  WM.PutCoef(9, 1, -a);
469  WM.PutCoef(9, 4, a);
470 
471  /* Q23 */
472  WM.PutCoef(2, 10, -1.);
473  WM.PutCoef(3, 10, 1.);
474 
475  a = A[5]*A[5];
476  WM.PutCoef(10, 2, -a);
477  WM.PutCoef(10, 3, a);
478 #endif /* VALVE_6 */
479 
480  for (int i = 0; i < LAST_Q; i++) {
481  WM.PutCoef(5+i, 5+i, 2.*fabs(q[i])/dKappa);
482  }
483 
484  return WorkMat;
485 }
486 
487 void
488 Control_valve2::Prepare(void)
489 {
490  doublereal p[LAST_N];
491  doublereal pm = 0.;
492 
493  for (int i = 0; i < LAST_N; i++) {
494  p[i] = pNode[i]->dGetX();
495  pm += p[i];
496  }
497  pm /= 4.;
498 
499  density = HF->dGetDensity(pm);
500 
501  dp[Q12] = p[N1]-p[N2];
502  dp[Q34] = p[N3]-p[N4];
503  dp[Q13] = p[N1]-p[N3];
504  dp[Q24] = p[N2]-p[N4];
505 #ifdef VALVE_6
506  dp[Q14] = p[N1]-p[N4];
507  dp[Q23] = p[N2]-p[N3];
508 #endif /* VALVE_6 */
509 
510  Stato = pGetDriveCaller()->dGet();
511  if (Stato > 0.) {
512  if (Stato > 1.) {
513  Stato = 1.;
514  }
515 
516  A[Q12] = Stato*area_max+2.*area_min;
517  A[Q34] = Stato*area_max+2.*area_min;
518  A[Q13] = area_min;
519  A[Q24] = area_min;
520  } else {
521  if (Stato < -1.) {
522  Stato = -1.;
523  }
524 
525  A[Q12] = 2.*area_min;
526  A[Q34] = 2.*area_min;
527  A[Q13] = -Stato*area_max+area_min;
528  A[Q24] = -Stato*area_max+area_min;
529  }
530 
531  /*
532  * vale sempre area_min (lo genero ogni volta perche' magari
533  * la si puo' far dipendere da qualche parametro
534  */
535 #ifdef VALVE_6
536  A[Q14] = area_min;
537  A[Q23] = area_min;
538 #endif /* VALVE_6 */
539 }
540 
541 SubVectorHandler&
542 Control_valve2::AssRes(SubVectorHandler& WorkVec,
543  doublereal dCoef,
544  const VectorHandler& XCurr,
545  const VectorHandler& XPrimeCurr)
546 {
547  DEBUGCOUT("Entering Control_valve2::AssRes()" << std::endl);
548 
549  integer iNumRows, iNumCols;
550  WorkSpaceDim(&iNumRows, &iNumCols);
551  WorkVec.Resize(iNumRows);
552 
553  integer iFirstIndex = iGetFirstIndex()+1;
554 
555  integer iNodeRowIndex[LAST_N];
556  for (int i = 0; i < LAST_N; i++) {
557  iNodeRowIndex[i] = pNode[i]->iGetFirstRowIndex()+1;
558  }
559 
560  Prepare();
561  doublereal dKappa = 2.*Cd*Cd*density;
562 
563  for (int i = 0; i < LAST_Q; i++) {
564  q[i] = XCurr(iFirstIndex+i);
565  }
566 
567 #ifdef VALVE_6
568  f[N1] = q[Q12]+q[Q13]+q[Q14];
569  f[N2] = -q[Q12]+q[Q23]+q[Q24];
570  f[N3] = -q[Q13]-q[Q23]+q[Q34];
571  f[N4] = -q[Q14]-q[Q24]-q[Q34];
572 #else /* !VALVE_6 */
573  f[N1] = q[Q12]+q[Q13];
574  f[N2] = -q[Q12]+q[Q24];
575  f[N3] = -q[Q13]+q[Q34];
576  f[N4] = -q[Q24]-q[Q34];
577 #endif /* !VALVE_6 */
578 
579  for (int i = 0; i < LAST_N; i++) {
580  WorkVec.PutItem(1+i, iNodeRowIndex[i], f[i]);
581  }
582 
583  for (int i = 0; i < LAST_Q; i++) {
584  WorkVec.PutItem(5+i, iFirstIndex+i,
585  A[i]*A[i]*dp[i]-q[i]*fabs(q[i])/dKappa);
586  }
587 
588  return WorkVec;
589 }
590 
591 void
592 Control_valve2::Output(OutputHandler& OH) const
593 {
594  if (bToBeOutput()) {
595  std::ostream& out = OH.Hydraulic();
596  out << std::setw(8) << GetLabel()
597  << " " << Stato
598  << " " << -f[N1] << " " << -f[N2]
599  << " " << -f[N3] << " " << -f[N4] << std::endl;
600  }
601 }
602 
603 void
604 Control_valve2::SetValue(DataManager *pDM,
605  VectorHandler& X, VectorHandler& /* XP */ ,
606  SimulationEntity::Hints *ph)
607 {
608  integer iFirstIndex = iGetFirstIndex()+1;
609 
610  const_cast<Control_valve2 *>(this)->Prepare();
611 
612  for (int i = 0; i < LAST_Q; i++) {
613 
614  /*
615  * q = sign(Dp)*Cd*A*sqrt(2.*rho*abs(Dp))
616  */
617 
618  X.PutCoef(iFirstIndex+i,
619  Cd*A[i]*copysign(sqrt(2.*density*fabs(dp[i])), dp[i]));
620  }
621 }
622 
623 /* Control_valve2 - end */
624 
625 
626 /* Dynamic_control_valve - begin */
627 
628 Dynamic_control_valve::Dynamic_control_valve(unsigned int uL,
629  const DofOwner* pDO,
630  HydraulicFluid* hf,
631  const PressureNode* p1,
632  const PressureNode* p2,
633  const PressureNode* p3,
634  const PressureNode* p4,
635  const DriveCaller* pDC,
636  doublereal s0, doublereal s_mx,
637  doublereal W, doublereal Loss_A,
638  doublereal Valve_d,
639  doublereal Valve_rho,
640  doublereal cs, doublereal cv,
641  doublereal ca, flag fOut)
642 : Elem(uL, fOut),
643 HydraulicElem(uL, pDO, hf, fOut), DriveOwner(pDC),
644 pNode1(p1), pNode2(p2), pNode3(p3), pNode4(p4),
645 start(s0), c_spost(cs), c_vel(cv), c_acc(ca),
646 width(W), loss_area(Loss_A),
647 valve_diameter(Valve_d), valve_density(Valve_rho), s_max(s_mx)
648 {
649  ASSERT(pNode1 != NULL);
650  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
651  ASSERT(pNode2 != NULL);
652  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
653  ASSERT(pNode3 != NULL);
654  ASSERT(pNode3->GetNodeType() == Node::HYDRAULIC);
655  ASSERT(pNode4 != NULL);
656  ASSERT(pNode4->GetNodeType() == Node::HYDRAULIC);
657  ASSERT(s0 >= 0.);
658  ASSERT(Valve_rho > std::numeric_limits<doublereal>::epsilon());
659  ASSERT(Valve_d > std::numeric_limits<doublereal>::epsilon());
660  ASSERT(W > std::numeric_limits<doublereal>::epsilon());
661  ASSERT(Loss_A >= 0.);
662  ASSERT(s_mx >= 0.);
663 
664  Cd = .611; /* coefficiente di perdita */
665  Mass = 1.175*valve_diameter*valve_diameter*valve_diameter*valve_density*M_PI; // massa della valvola
666 
667  // W = .005; /* larghezza del condotto (m): A=x*W 0.005; */
668  // diameter = .01; /* diametro della valvola (m) */
669  // densityvalve = 7900.; /* densita' del corpo della valvola (acciaio) (kg/m^3) */
670  // Mass = 1.175*diameter*diameter*diameter*densityvalve*M_PI; /* massa della valvola (Kg) */
671  // s_max = 2.; /* corsa massima della valvola (m) */
672 }
673 
674 
675 Dynamic_control_valve::~Dynamic_control_valve(void)
676 {
677  NO_OP;
678 }
679 
680 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
681 HydraulicElem::Type Dynamic_control_valve::GetHydraulicType(void) const
682 {
683  return HydraulicElem::DYNAMIC_CONTROL_VALVE;
684 }
685 
686 /* Contributo al file di restart */
687 std::ostream& Dynamic_control_valve::Restart(std::ostream& out) const
688 {
689  return out << "Dynamic_control_valve not implemented yet!" << std::endl;
690 }
691 
692 unsigned int Dynamic_control_valve::iGetNumDof(void) const
693 {
694  return 2;
695 }
696 
697 DofOrder::Order Dynamic_control_valve::GetDofType(unsigned int i) const
698 {
699  ASSERT(i >= 0 && i <= 1);
700  return DofOrder::DIFFERENTIAL;
701 }
702 
703 
704 void
705 Dynamic_control_valve::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
706 {
707  *piNumRows = 6;
708  *piNumCols = 6;
709 }
710 
711 VariableSubMatrixHandler&
712 Dynamic_control_valve::AssJac(VariableSubMatrixHandler& WorkMat,
713  doublereal dCoef,
714  const VectorHandler& XCurr,
715  const VectorHandler& XPrimeCurr)
716 {
717  DEBUGCOUT("Entering Control_valve::AssJac()" << std::endl);
718 
719  FullSubMatrixHandler& WM = WorkMat.SetFull();
720  WM.ResizeReset(6, 6);
721 
722  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
723  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
724  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
725  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
726  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
727  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
728  integer iNode3ColIndex = pNode3->iGetFirstColIndex()+1;
729  integer iNode4ColIndex = pNode4->iGetFirstColIndex()+1;
730  integer iFirstIndex = iGetFirstIndex();
731 
732  WM.PutRowIndex(1, iNode1RowIndex);
733  WM.PutRowIndex(2, iNode2RowIndex);
734  WM.PutRowIndex(3, iNode3RowIndex);
735  WM.PutRowIndex(4, iNode4RowIndex);
736  WM.PutColIndex(1, iNode1ColIndex);
737  WM.PutColIndex(2, iNode2ColIndex);
738  WM.PutColIndex(3, iNode3ColIndex);
739  WM.PutColIndex(4, iNode4ColIndex);
740 
741  WM.PutRowIndex(5, iFirstIndex+1);
742  WM.PutColIndex(5, iFirstIndex+1);
743  WM.PutRowIndex(6, iFirstIndex+2);
744  WM.PutColIndex(6, iFirstIndex+2);
745 
746  doublereal p1 = pNode1->dGetX();
747  doublereal p2 = pNode2->dGetX();
748  doublereal p3 = pNode3->dGetX();
749  doublereal p4 = pNode4->dGetX();
750  doublereal density = HF->dGetDensity();
751 
752  s = XCurr(iFirstIndex+1); /* spostamento */
753  v = XCurr(iFirstIndex+2); /* velocita' */
754 
755  doublereal jumpPres12 = fabs(p1-p2); /* salto di pressione nodo1 & nodo2 */
756  doublereal jumpPres13 = fabs(p1-p3); /* salto di pressione nodo1 & nodo3 */
757  doublereal jumpPres24 = fabs(p2-p4); /* salto di pressione nodo2 & nodo4 */
758  doublereal jumpPres34 = fabs(p3-p4); /* salto di pressione nodo3 & nodo4 */
759 
760  if (jumpPres12 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
761  jumpPres12 = 1.e8*std::numeric_limits<doublereal>::epsilon();
762  }
763  if (jumpPres13 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
764  jumpPres13 = 1.e8*std::numeric_limits<doublereal>::epsilon();
765  }
766  if (jumpPres24 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
767  jumpPres24 = 1.e8*std::numeric_limits<doublereal>::epsilon();
768  }
769  if (jumpPres34 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
770  jumpPres34 = 1.e8*std::numeric_limits<doublereal>::epsilon();
771  }
772  doublereal Jac15;
773  doublereal Jac25;
774  doublereal Jac35;
775  doublereal Jac45;
776 
777  doublereal primo = Cd*A1/sqrt(2.*jumpPres12/density);
778  doublereal secondo = Cd*A2/sqrt(2.*jumpPres13/density);
779  doublereal terzo = Cd*A4/sqrt(2.*jumpPres24/density);
780  doublereal quarto = Cd*A3/sqrt(2.*jumpPres34/density);
781 
782  doublereal Jac11 = -primo-secondo;
783  doublereal Jac12 = primo;
784  doublereal Jac13 = secondo;
785  doublereal Jac14 = 0.;
786  doublereal Jac21 = primo;
787  doublereal Jac22 = -primo-terzo;
788  doublereal Jac23 = 0.;
789  doublereal Jac24 = terzo;
790  doublereal Jac31 = secondo;
791  doublereal Jac32 = 0.;
792  doublereal Jac33 = -secondo-quarto;
793  doublereal Jac34 = quarto;
794  doublereal Jac41 = 0.;
795  doublereal Jac42 = terzo;
796  doublereal Jac43 = quarto;
797  doublereal Jac44 = -terzo-quarto;
798 
799  doublereal costante = density*Cd*width*valve_diameter;
800 
801 #if 0
802  doublereal Jac51 = -.2*costante*sp/sqrt(density*deltaP)-.43*width*s;
803  doublereal Jac52 = -Jac51;
804  doublereal Jac53 = Jac51;
805  doublereal Jac54 = -Jac51;
806  doublereal Jac55 = -.4*valve_diameter*Cd*width*sqrt(density*deltaP)-dCoef*.43*width*deltaP-dCoef*c1-c2-dCoef*cf1-cf2;
807  doublereal Jac56 = -Mass-c3-cf3;
808 #endif
809 
810  doublereal Jac51 = -.2*costante*sp/sqrt(density*deltaP)-.43*width*s;
811  doublereal Jac52 = 0.;
812  doublereal Jac53 = 0.;
813  doublereal Jac54 = 0.;
814  doublereal Jac55 = -.4*valve_diameter*Cd*width*sqrt(density*deltaP)-dCoef*.43*width*deltaP-dCoef*c1-c2-dCoef*cf1-cf2;
815  doublereal Jac56 = -Mass-c3-cf3;
816 
817  if (s > 0.) { /* collegamento diretto 1->2 3->4 */
818  Jac15 = dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres12/density), p1-p2);
819  Jac25 = -Jac15;
820  Jac35 = dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres34/density), p3-p4);
821  Jac45 = -Jac35;
822  } else { /* collegamento inverso 1->3 2->4 */
823  Jac15 = -dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres13/density), p1-p3);
824  Jac25 = -dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres24/density), p2-p4);
825  Jac35 = -Jac15;
826  Jac45 = -Jac25;
827  }
828 
829  doublereal Jac65 = -1;
830  doublereal Jac66 = dCoef;
831 
832  WM.PutCoef(1, 1, Jac11);
833  WM.PutCoef(1, 2, Jac12);
834  WM.PutCoef(1, 3, Jac13);
835  WM.PutCoef(1, 4, Jac14);
836  WM.PutCoef(1, 5, Jac15);
837  WM.PutCoef(2, 1, Jac21);
838  WM.PutCoef(2, 2, Jac22);
839  WM.PutCoef(2, 3, Jac23);
840  WM.PutCoef(2, 4, Jac24);
841  WM.PutCoef(2, 5, Jac25);
842  WM.PutCoef(3, 1, Jac31);
843  WM.PutCoef(3, 2, Jac32);
844  WM.PutCoef(3, 3, Jac33);
845  WM.PutCoef(3, 4, Jac34);
846  WM.PutCoef(3, 5, Jac35);
847  WM.PutCoef(4, 1, Jac41);
848  WM.PutCoef(4, 2, Jac42);
849  WM.PutCoef(4, 3, Jac43);
850  WM.PutCoef(4, 4, Jac44);
851  WM.PutCoef(4, 5, Jac45);
852  WM.PutCoef(5, 1, Jac51);
853  WM.PutCoef(5, 2, Jac52);
854  WM.PutCoef(5, 3, Jac53);
855  WM.PutCoef(5, 5, Jac55);
856  WM.PutCoef(5, 6, Jac56);
857  WM.PutCoef(6, 5, Jac65);
858  WM.PutCoef(6, 6, Jac66);
859 
860  return WorkMat;
861 }
862 
863 SubVectorHandler&
864 Dynamic_control_valve::AssRes(SubVectorHandler& WorkVec,
865  doublereal dCoef,
866  const VectorHandler& XCurr,
867  const VectorHandler& XPrimeCurr)
868 {
869  DEBUGCOUT("Entering Dynamic_control_valve::AssRes()" << std::endl);
870 
871  WorkVec.Resize(6);
872 
873  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
874  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
875  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
876  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
877  integer iFirstIndex = iGetFirstIndex();
878 
879  doublereal p1 = pNode1->dGetX();
880  doublereal p2 = pNode2->dGetX();
881  doublereal p3 = pNode3->dGetX();
882  doublereal p4 = pNode4->dGetX();
883  doublereal density = HF->dGetDensity();
884 
885  Force = pGetDriveCaller()->dGet();
886 
887  s = XCurr(iFirstIndex+1); /* spostamento */
888  v = XCurr(iFirstIndex+2); /* velocita' */
889  sp = XPrimeCurr(iFirstIndex+1); /* velocita' */
890  vp = XPrimeCurr(iFirstIndex+2); /* accelerazione */
891 
892  doublereal jumpPres12 = fabs(p1-p2); /* salto di pressione nodo1 & nodo2 */
893  doublereal jumpPres13 = fabs(p1-p3); /* salto di pressione nodo1 & nodo3 */
894  doublereal jumpPres24 = fabs(p2-p4); /* salto di pressione nodo2 & nodo4 */
895  doublereal jumpPres34 = fabs(p3-p4); /* salto di pressione nodo3 & nodo4 */
896 
897  /* qui decido se sono a fondo od ad inizio corsa e di conseguenza calcolo
898  * i giusti coefficienti */
899 
900 
901  doublereal area_max = width*s_max; // Area massima
902  doublereal area_min = area_max*loss_area; // Area minimo
903  doublereal deltaA = s*width;
904 
905  if (s > 0.) {
906  A1 = area_min+deltaA;
907  A2 = area_min;
908  A3 = area_min+deltaA;
909  A4 = area_min;
910  } else { /* ho deltaA negativo */
911  A1 = area_min;
912  A2 = area_min-deltaA;
913  A3 = area_min;
914  A4 = area_min-deltaA;
915  }
916 
917  if (s <= -s_max) {
918  c1 = c_spost;
919 
920  if (sp < 0.) {
921  c2 = c_vel; /* ho v negativa devo smorzare */
922  } else {
923  c2 = 0.;
924  }
925 
926  if (vp < 0.) {
927  c3 = c_acc; /* ho acc negativa devo smorzare */
928  } else {
929  c3 = 0.;
930  }
931 
932  } else {
933  c1 = 0.;
934  c2 = 0.;
935  c3 = 0.;
936  }
937 
938  if (s >= s_max) {
939  cf1 = c_spost;
940 
941  if (sp > 0.) {
942  cf2 = c_vel; /* ho v positiva devo smorzare */
943  } else {
944  cf2 = 0.;
945  }
946 
947  if (vp > 0.) {
948  cf3 = c_acc; /* ho acc positiva devo smorzare */
949  } else {
950  cf3 = 0.;
951  }
952 
953  } else {
954  cf1 = 0.;
955  cf2 = 0.;
956  cf3 = 0.;
957  }
958 
959 
960  doublereal Q12 = density*Cd*A1*copysign(sqrt(2.*jumpPres12/density), p1-p2);
961  doublereal Q13 = density*Cd*A2*copysign(sqrt(2.*jumpPres13/density), p1-p3);
962  doublereal Q24 = density*Cd*A4*copysign(sqrt(2.*jumpPres24/density), p2-p4);
963  doublereal Q34 = density*Cd*A3*copysign(sqrt(2.*jumpPres34/density), p3-p4);
964 
965  doublereal Res_1 = Q12+Q13;
966  doublereal Res_2 = -Q12+Q24;
967  doublereal Res_3 = -Q13+Q34;
968  doublereal Res_4 = -Q34-Q24;
969 
970  deltaP =p1;
971 
972  if (deltaP == 0.) {
973  deltaP = 10.; /* evito di dividere per zero nello jacobiano */
974  }
975 
976  doublereal C = .4*valve_diameter*Cd*width*sqrt(density*deltaP);
977  doublereal K = .43*width*deltaP;
978  doublereal Res_5 = -Force+Mass*vp+C*sp+K*s+c1*(s+s_max)+c2*sp+c3*vp+cf1*(s-s_max)+cf2*sp+cf3*vp;
979  doublereal Res_6 = sp-v;
980 
981  flow1 = -Res_1;
982  flow2 = -Res_2;
983  flow3 = -Res_3;
984  flow4 = -Res_4;
985 
986 #ifdef HYDR_DEVEL
987  DEBUGCOUT("Force: " << Force << std::endl);
988  DEBUGCOUT("X equilibrio: " << Force/K << std::endl);
989  DEBUGCOUT("A1: " << A1 << std::endl);
990  DEBUGCOUT("A2: " << A2 << std::endl);
991  DEBUGCOUT("A3: " << A3 << std::endl);
992  DEBUGCOUT("A4: " << A4 << std::endl);
993  DEBUGCOUT("p1: " << p1 << std::endl);
994  DEBUGCOUT("p2: " << p2 << std::endl);
995  DEBUGCOUT("p3: " << p3 << std::endl);
996  DEBUGCOUT("p4: " << p4 << std::endl);
997  DEBUGCOUT("Cd: " << Cd << std::endl);
998  DEBUGCOUT("s_max : " << s_max << std::endl);
999  DEBUGCOUT("s : " << s << std::endl);
1000  DEBUGCOUT("sp: " << sp << std::endl);
1001  DEBUGCOUT("v : " << v << std::endl);
1002  DEBUGCOUT("vp: " << vp << std::endl);
1003  DEBUGCOUT("Valve_diameter:" << valve_diameter << std::endl);
1004  DEBUGCOUT("massa: " << Mass << std::endl);
1005  DEBUGCOUT("smorzatore: " << C << std::endl);
1006  DEBUGCOUT("molla: " << K << std::endl);
1007  DEBUGCOUT("density: " << density << std::endl);
1008  DEBUGCOUT("Valve_density: " << valve_density << std::endl);
1009  DEBUGCOUT("Area_max: " << area_max << std::endl);
1010  DEBUGCOUT("Width: " << width << std::endl);
1011  DEBUGCOUT("Loss_area: " << loss_area << std::endl);
1012  DEBUGCOUT("c1: " << c1 << std::endl);
1013  DEBUGCOUT("c2: " << c2 << std::endl);
1014  DEBUGCOUT("c3: " << c3 << std::endl);
1015  DEBUGCOUT("cf1: " << cf1 << std::endl);
1016  DEBUGCOUT("cf2: " << cf2 << std::endl);
1017  DEBUGCOUT("cf3: " << cf3 << std::endl);
1018  DEBUGCOUT("Q12: " << Q12 << std::endl);
1019  DEBUGCOUT("Q13: " << Q13 << std::endl);
1020  DEBUGCOUT("Q24: " << Q24 << std::endl);
1021  DEBUGCOUT("Q34: " << Q34 << std::endl);
1022  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)"<< std::endl);
1023  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
1024  DEBUGCOUT("-Res_2 (portata nodo2): " << -Res_2 << std::endl);
1025  DEBUGCOUT("-Res_3 (portata nodo3): " << -Res_3 << std::endl);
1026  DEBUGCOUT("-Res_4 (portata nodo4): " << -Res_4 << std::endl);
1027  DEBUGCOUT("-Res_5 eq.dinamica : " << -Res_5 << std::endl);
1028  DEBUGCOUT("-Res_6 sp-v : " << -Res_6 << std::endl);
1029 #endif
1030 
1031  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
1032  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
1033  WorkVec.PutItem(3, iNode3RowIndex, Res_3);
1034  WorkVec.PutItem(4, iNode4RowIndex, Res_4);
1035  WorkVec.PutItem(5, iFirstIndex+1, Res_5);
1036  WorkVec.PutItem(6, iFirstIndex+2, Res_6);
1037 
1038  return WorkVec;
1039 }
1040 
1041 
1042 void Dynamic_control_valve::Output(OutputHandler& OH) const
1043 {
1044  if (bToBeOutput()) {
1045  std::ostream& out = OH.Hydraulic();
1046  out << std::setw(8) << GetLabel()
1047  << " " << s << " " << sp << " " << vp
1048  << " " << flow1 << " " << flow2 << " " << flow3 << " " << flow4
1049  << " " << A1 << " " << A2 << " " << A3 << " " << A4 << " " << std::endl;
1050  }
1051 }
1052 
1053 
1054 void
1055 Dynamic_control_valve::SetValue(DataManager *pDM,
1056  VectorHandler& X , VectorHandler& XP,
1057  SimulationEntity::Hints *ph)
1058 {
1059  integer i = iGetFirstIndex();
1060 
1061  X.PutCoef(i+1, start);
1062  X.PutCoef(i+2, 0.);
1063  XP.PutCoef(i+1, 0.);
1064  XP.PutCoef(i+2, 0.);
1065 }
1066 
1067 /* Dynamic_control_valve - end */
1068 
1069 
1070 /* Pressure_flow_control_valve - begin */
1071 
1072 Pressure_flow_control_valve::Pressure_flow_control_valve(unsigned int uL, const DofOwner* pDO,
1073  HydraulicFluid* hf,
1074  const PressureNode* p1,
1075  const PressureNode* p2,
1076  const PressureNode* p3,
1077  const PressureNode* p4,
1078  const PressureNode* p5,
1079  const PressureNode* p6,
1080  const DriveCaller* pDC,
1081  doublereal s0, doublereal s_mx, doublereal W, doublereal Loss_A,
1082  doublereal Valve_d, doublereal Valve_rho,
1083  doublereal cs, doublereal cv,
1084  doublereal ca, flag fOut)
1085 : Elem(uL, fOut),
1086 HydraulicElem(uL, pDO, hf, fOut), DriveOwner(pDC),
1087 pNode1(p1), pNode2(p2), pNode3(p3), pNode4(p4), pNode5(p5), pNode6(p6),
1088 start(s0), s_max(s_mx), width(W), loss_area(Loss_A),
1089 valve_diameter(Valve_d), valve_density(Valve_rho),
1090 c_spost(cs), c_vel(cv), c_acc(ca)
1091 {
1092  ASSERT(pNode1 != NULL);
1093  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
1094  ASSERT(pNode2 != NULL);
1095  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
1096  ASSERT(pNode3 != NULL);
1097  ASSERT(pNode3->GetNodeType() == Node::HYDRAULIC);
1098  ASSERT(pNode4 != NULL);
1099  ASSERT(pNode4->GetNodeType() == Node::HYDRAULIC);
1100  ASSERT(pNode5 != NULL);
1101  ASSERT(pNode5->GetNodeType() == Node::HYDRAULIC);
1102  ASSERT(pNode6 != NULL);
1103  ASSERT(pNode6->GetNodeType() == Node::HYDRAULIC);
1104  ASSERT(s0 >= 0.);
1105  ASSERT(Valve_rho > std::numeric_limits<doublereal>::epsilon());
1106  ASSERT(Valve_d > std::numeric_limits<doublereal>::epsilon());
1107  ASSERT(W > std::numeric_limits<doublereal>::epsilon());
1108  ASSERT(Loss_A >= 0.);
1109  ASSERT(s_mx >= 0.);
1110 
1111  Cd = .611; /* coefficiente di perdita */
1112  Mass = 1.175*valve_diameter*valve_diameter*valve_diameter*valve_density*M_PI; // massa della valvola
1113  valve_area=0.25*valve_diameter*valve_diameter*M_PI;
1114 }
1115 
1116 
1117 Pressure_flow_control_valve::~Pressure_flow_control_valve(void)
1118 {
1119  NO_OP;
1120 }
1121 
1122 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
1123 HydraulicElem::Type Pressure_flow_control_valve::GetHydraulicType(void) const
1124 {
1125  return HydraulicElem::PRESSURE_FLOW_CONTROL_VALVE;
1126 }
1127 
1128 /* Contributo al file di restart */
1129 std::ostream& Pressure_flow_control_valve::Restart(std::ostream& out) const
1130 {
1131  return out << "Pressure_flow_control_valve not implemented yet!" << std::endl;
1132 }
1133 
1134 unsigned int Pressure_flow_control_valve::iGetNumDof(void) const
1135 {
1136  return 2;
1137 }
1138 
1139 DofOrder::Order Pressure_flow_control_valve::GetDofType(unsigned int i) const
1140 {
1141  ASSERT(i >= 0 && i <= 1);
1142  return DofOrder::DIFFERENTIAL;
1143 }
1144 
1145 
1146 void
1147 Pressure_flow_control_valve::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
1148 {
1149  *piNumRows = 8;
1150  *piNumCols = 8;
1151 }
1152 
1153 VariableSubMatrixHandler&
1154 Pressure_flow_control_valve::AssJac(VariableSubMatrixHandler& WorkMat,
1155  doublereal dCoef,
1156  const VectorHandler& XCurr,
1157  const VectorHandler& XPrimeCurr)
1158 {
1159  DEBUGCOUT("Entering Control_valve::AssJac()" << std::endl);
1160 
1161  FullSubMatrixHandler& WM = WorkMat.SetFull();
1162  WM.ResizeReset(8, 8);
1163 
1164  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1165  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1166  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
1167  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
1168  integer iNode5RowIndex = pNode5->iGetFirstRowIndex()+1;
1169  integer iNode6RowIndex = pNode6->iGetFirstRowIndex()+1;
1170  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
1171  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
1172  integer iNode3ColIndex = pNode3->iGetFirstColIndex()+1;
1173  integer iNode4ColIndex = pNode4->iGetFirstColIndex()+1;
1174  integer iNode5ColIndex = pNode5->iGetFirstColIndex()+1;
1175  integer iNode6ColIndex = pNode6->iGetFirstColIndex()+1;
1176  integer iFirstIndex = iGetFirstIndex();
1177 
1178  WM.PutRowIndex(1, iNode1RowIndex);
1179  WM.PutRowIndex(2, iNode2RowIndex);
1180  WM.PutRowIndex(3, iNode3RowIndex);
1181  WM.PutRowIndex(4, iNode4RowIndex);
1182  WM.PutRowIndex(5, iNode5RowIndex);
1183  WM.PutRowIndex(6, iNode6RowIndex);
1184  WM.PutColIndex(1, iNode1ColIndex);
1185  WM.PutColIndex(2, iNode2ColIndex);
1186  WM.PutColIndex(3, iNode3ColIndex);
1187  WM.PutColIndex(4, iNode4ColIndex);
1188  WM.PutColIndex(5, iNode5ColIndex);
1189  WM.PutColIndex(6, iNode6ColIndex);
1190 
1191  WM.PutRowIndex(7, iFirstIndex+1);
1192  WM.PutColIndex(7, iFirstIndex+1);
1193  WM.PutRowIndex(8, iFirstIndex+2);
1194  WM.PutColIndex(8, iFirstIndex+2);
1195 
1196  doublereal p1 = pNode1->dGetX();
1197  doublereal p2 = pNode2->dGetX();
1198  doublereal p3 = pNode3->dGetX();
1199  doublereal p4 = pNode4->dGetX();
1200  doublereal p5 = pNode5->dGetX();
1201  doublereal p6 = pNode6->dGetX();
1202  doublereal density = HF->dGetDensity();
1203 
1204  s = XCurr(iFirstIndex+1); /* spostamento */
1205  v = XCurr(iFirstIndex+2); /* velocita' */
1206 
1207  doublereal jumpPres12 = fabs(p1-p2); /* salto di pressione nodo1 & nodo2 */
1208  doublereal jumpPres13 = fabs(p1-p3); /* salto di pressione nodo1 & nodo3 */
1209  doublereal jumpPres24 = fabs(p2-p4); /* salto di pressione nodo2 & nodo4 */
1210  doublereal jumpPres34 = fabs(p3-p4); /* salto di pressione nodo3 & nodo4 */
1211 
1212  if (jumpPres12 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1213  jumpPres12 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1214  }
1215  if (jumpPres13 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1216  jumpPres13 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1217  }
1218  if (jumpPres24 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1219  jumpPres24 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1220  }
1221  if (jumpPres34 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1222  jumpPres34 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1223  }
1224  doublereal Jac17;
1225  doublereal Jac27;
1226  doublereal Jac37;
1227  doublereal Jac47;
1228 
1229  doublereal primo = Cd*A1/sqrt(2.*jumpPres12/density);
1230  doublereal secondo = Cd*A2/sqrt(2.*jumpPres13/density);
1231  doublereal terzo = Cd*A4/sqrt(2.*jumpPres24/density);
1232  doublereal quarto = Cd*A3/sqrt(2.*jumpPres34/density);
1233 
1234  doublereal Jac11 = -primo-secondo;
1235  doublereal Jac12 = primo;
1236  doublereal Jac13 = secondo;
1237  doublereal Jac14 = 0.;
1238  doublereal Jac21 = primo;
1239  doublereal Jac22 = -primo-terzo;
1240  doublereal Jac23 = 0.;
1241  doublereal Jac24 = terzo;
1242  doublereal Jac31 = secondo;
1243  doublereal Jac32 = 0.;
1244  doublereal Jac33 = -secondo-quarto;
1245  doublereal Jac34 = quarto;
1246  doublereal Jac41 = 0.;
1247  doublereal Jac42 = terzo;
1248  doublereal Jac43 = quarto;
1249  doublereal Jac44 = -terzo-quarto;
1250 
1251  doublereal Jac57 = -valve_area;
1252  doublereal Jac67 = +valve_area;
1253 
1254  doublereal costante = density*Cd*width*valve_diameter;
1255 
1256  doublereal Jac71 = -.2*costante/sqrt(density*deltaP)*sp-.43*width*s;
1257  doublereal Jac72 = 0.;
1258  doublereal Jac73 = 0.;
1259  doublereal Jac74 = 0.;
1260  doublereal Jac77 = -.4*valve_diameter*Cd*width*sqrt(density*deltaP)-dCoef*.43*width*deltaP-dCoef*c1-c2-dCoef*cf1-cf2;
1261  doublereal Jac78 = -Mass-c3-cf3;
1262 
1263  if (s > 0.) { /* collegamento diretto 1->2 3->4 */
1264  Jac17 = dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres12/density), p1-p2);
1265  Jac27 = -Jac17;
1266  Jac37 = dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres34/density), p3-p4);
1267  Jac47 = -Jac37;
1268  } else { /* collegamento inverso 1->3 2->4 */
1269  Jac17 = -dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres13/density), p1-p3);
1270  Jac27 = -dCoef*density*Cd*width*copysign(sqrt(2.*jumpPres24/density), p2-p4);
1271  Jac37 = -Jac17;
1272  Jac47 = -Jac27;
1273  }
1274 
1275  doublereal Jac87 = -1;
1276  doublereal Jac88 = dCoef;
1277 
1278  WM.PutCoef(1, 1, Jac11);
1279  WM.PutCoef(1, 2, Jac12);
1280  WM.PutCoef(1, 3, Jac13);
1281  WM.PutCoef(1, 4, Jac14);
1282  WM.PutCoef(1, 7, Jac17);
1283  WM.PutCoef(2, 1, Jac21);
1284  WM.PutCoef(2, 2, Jac22);
1285  WM.PutCoef(2, 3, Jac23);
1286  WM.PutCoef(2, 4, Jac24);
1287  WM.PutCoef(2, 7, Jac27);
1288  WM.PutCoef(3, 1, Jac31);
1289  WM.PutCoef(3, 2, Jac32);
1290  WM.PutCoef(3, 3, Jac33);
1291  WM.PutCoef(3, 4, Jac34);
1292  WM.PutCoef(3, 7, Jac37);
1293  WM.PutCoef(4, 1, Jac41);
1294  WM.PutCoef(4, 2, Jac42);
1295  WM.PutCoef(4, 3, Jac43);
1296  WM.PutCoef(4, 4, Jac44);
1297  WM.PutCoef(4, 7, Jac47);
1298  WM.PutCoef(5, 7, Jac67);
1299  WM.PutCoef(6, 7, Jac67);
1300  WM.PutCoef(7, 1, Jac71);
1301  WM.PutCoef(7, 2, Jac72);
1302  WM.PutCoef(7, 3, Jac73);
1303  WM.PutCoef(7, 7, Jac77);
1304  WM.PutCoef(7, 8, Jac78);
1305  WM.PutCoef(8, 7, Jac87);
1306  WM.PutCoef(8, 8, Jac88);
1307 
1308  return WorkMat;
1309 }
1310 
1311 SubVectorHandler&
1312 Pressure_flow_control_valve::AssRes(SubVectorHandler& WorkVec,
1313  doublereal dCoef,
1314  const VectorHandler& XCurr,
1315  const VectorHandler& XPrimeCurr)
1316 {
1317  DEBUGCOUT("Entering Pressure_flow_control_valve::AssRes()" << std::endl);
1318 
1319  WorkVec.Resize(8);
1320 
1321  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1322  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1323  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
1324  integer iNode4RowIndex = pNode4->iGetFirstRowIndex()+1;
1325  integer iNode5RowIndex = pNode5->iGetFirstRowIndex()+1;
1326  integer iNode6RowIndex = pNode6->iGetFirstRowIndex()+1;
1327  integer iFirstIndex = iGetFirstIndex();
1328 
1329  doublereal p1 = pNode1->dGetX();
1330  doublereal p2 = pNode2->dGetX();
1331  doublereal p3 = pNode3->dGetX();
1332  doublereal p4 = pNode4->dGetX();
1333  doublereal p5 = pNode5->dGetX();
1334  doublereal p6 = pNode6->dGetX();
1335  doublereal density = HF->dGetDensity();
1336 
1337  Force = pGetDriveCaller()->dGet();
1338 
1339  s = XCurr(iFirstIndex+1); /* spostamento */
1340  v = XCurr(iFirstIndex+2); /* velocita' */
1341  sp = XPrimeCurr(iFirstIndex+1); /* velocita' */
1342  vp = XPrimeCurr(iFirstIndex+2); /* accelerazione */
1343 
1344  doublereal jumpPres12 = fabs(p1-p2); /* salto di pressione nodo1 & nodo2 */
1345  doublereal jumpPres13 = fabs(p1-p3); /* salto di pressione nodo1 & nodo3 */
1346  doublereal jumpPres24 = fabs(p2-p4); /* salto di pressione nodo2 & nodo4 */
1347  doublereal jumpPres34 = fabs(p3-p4); /* salto di pressione nodo3 & nodo4 */
1348 
1349  /* qui decido se sono a fondo od ad inizio corsa e di conseguenza calcolo
1350  * i giusti coefficienti */
1351 
1352  doublereal x = s; /* usato per calcolare le aree */
1353 
1354  doublereal area_max = width*s_max; // Area massima
1355  doublereal area_min = area_max*loss_area; // Area minimo
1356  doublereal deltaA = x*width;
1357 
1358  if (x > 0.) {
1359  A1 = area_min+deltaA;
1360  A2 = area_min;
1361  A3 = area_min+deltaA;
1362  A4 = area_min;
1363  } else { /* ho deltaA negativo */
1364  A1 = area_min;
1365  A2 = area_min-deltaA;
1366  A3 = area_min;
1367  A4 = area_min-deltaA;
1368  }
1369 
1370  doublereal Q12 = density*Cd*A1*copysign(sqrt(2.*jumpPres12/density), p1-p2);
1371  doublereal Q13 = density*Cd*A2*copysign(sqrt(2.*jumpPres13/density), p1-p3);
1372  doublereal Q24 = density*Cd*A4*copysign(sqrt(2.*jumpPres24/density), p2-p4);
1373  doublereal Q34 = density*Cd*A3*copysign(sqrt(2.*jumpPres34/density), p3-p4);
1374 
1375  doublereal Res_1 = Q12+Q13;
1376  doublereal Res_2 = -Q12+Q24;
1377  doublereal Res_3 = -Q13+Q34;
1378  doublereal Res_4 = -Q34-Q24;
1379  doublereal Res_5 = valve_area*sp;
1380 #warning "????????????? Res_6 = -Res_6 ?"
1381  doublereal Res_6 = -Res_6;
1382 
1383  deltaP = p1;
1384  if (deltaP == 0.) {
1385  deltaP = 10.; /* evito di dividere per zero nello jacobiano */
1386  }
1387  doublereal C = .4*valve_diameter*Cd*width*sqrt(density*deltaP);
1388  doublereal K = .43*width*deltaP;
1389  doublereal Res_7 = -Force+Mass*vp+C*sp+K*s+c1*s+c2*sp+c3*vp+cf1*(s-s_max)+cf2*sp+cf3*vp;
1390  doublereal Res_8 = sp-v;
1391 
1392  flow1 = -Res_1;
1393  flow2 = -Res_2;
1394  flow3 = -Res_3;
1395  flow4 = -Res_4;
1396  flow5 = -Res_5;
1397  flow6 = -Res_6;
1398 
1399 #ifdef HYDR_DEVEL
1400  DEBUGCOUT("Force: " << Force << std::endl);
1401  DEBUGCOUT("X equilibrio: " << Force/K << std::endl);
1402  DEBUGCOUT("A1: " << A1 << std::endl);
1403  DEBUGCOUT("A2: " << A2 << std::endl);
1404  DEBUGCOUT("A3: " << A3 << std::endl);
1405  DEBUGCOUT("A4: " << A4 << std::endl);
1406  DEBUGCOUT("p1: " << p1 << std::endl);
1407  DEBUGCOUT("p2: " << p2 << std::endl);
1408  DEBUGCOUT("p3: " << p3 << std::endl);
1409  DEBUGCOUT("p4: " << p4 << std::endl);
1410  DEBUGCOUT("Cd: " << Cd << std::endl);
1411  DEBUGCOUT("s_max : " << s_max << std::endl);
1412  DEBUGCOUT("x : " << x << std::endl);
1413  DEBUGCOUT("s : " << s << std::endl);
1414  DEBUGCOUT("sp: " << sp << std::endl);
1415  DEBUGCOUT("v : " << v << std::endl);
1416  DEBUGCOUT("vp: " << vp << std::endl);
1417  DEBUGCOUT("Valve_diameter:" << valve_diameter << std::endl);
1418  DEBUGCOUT("massa: " << Mass << std::endl);
1419  DEBUGCOUT("smorzatore: " << C << std::endl);
1420  DEBUGCOUT("molla: " << K << std::endl);
1421  DEBUGCOUT("density: " << density << std::endl);
1422  DEBUGCOUT("Valve_density: " << valve_density << std::endl);
1423  DEBUGCOUT("Area_max: " << area_max << std::endl);
1424  DEBUGCOUT("Width: " << width << std::endl);
1425  DEBUGCOUT("Loss_area: " << loss_area << std::endl);
1426  DEBUGCOUT("c1: " << c1 << std::endl);
1427  DEBUGCOUT("c2: " << c2 << std::endl);
1428  DEBUGCOUT("c3: " << c3 << std::endl);
1429  DEBUGCOUT("cf1: " << cf1 << std::endl);
1430  DEBUGCOUT("cf2: " << cf2 << std::endl);
1431  DEBUGCOUT("cf3: " << cf3 << std::endl);
1432  DEBUGCOUT("Q12: " << Q12 << std::endl);
1433  DEBUGCOUT("Q13: " << Q13 << std::endl);
1434  DEBUGCOUT("Q24: " << Q24 << std::endl);
1435  DEBUGCOUT("Q34: " << Q34 << std::endl);
1436  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)"<< std::endl);
1437  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
1438  DEBUGCOUT("-Res_2 (portata nodo2): " << -Res_2 << std::endl);
1439  DEBUGCOUT("-Res_3 (portata nodo3): " << -Res_3 << std::endl);
1440  DEBUGCOUT("-Res_4 (portata nodo4): " << -Res_4 << std::endl);
1441  DEBUGCOUT("-Res_5 (portata nodo3): " << -Res_5 << std::endl);
1442  DEBUGCOUT("-Res_6 (portata nodo4): " << -Res_6 << std::endl);
1443  DEBUGCOUT("-Res_6 eq dinamica : " << -Res_7 << std::endl);
1444  DEBUGCOUT("-Res_7 sp-v : " << -Res_8 << std::endl);
1445 #endif
1446 
1447  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
1448  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
1449  WorkVec.PutItem(3, iNode3RowIndex, Res_3);
1450  WorkVec.PutItem(4, iNode4RowIndex, Res_4);
1451  WorkVec.PutItem(5, iNode5RowIndex, Res_5);
1452  WorkVec.PutItem(6, iNode6RowIndex, Res_6);
1453  WorkVec.PutItem(7, iFirstIndex+1, Res_7);
1454  WorkVec.PutItem(8, iFirstIndex+2, Res_8);
1455 
1456  return WorkVec;
1457 }
1458 
1459 
1460 void Pressure_flow_control_valve::Output(OutputHandler& OH) const
1461 {
1462  if (bToBeOutput()) {
1463  std::ostream& out = OH.Hydraulic();
1464  out << std::setw(8) << GetLabel()
1465  << " " << s << " " << sp << " " << vp
1466  << " " << flow1 << " " << flow2 << " " << flow3 << " " << flow4
1467  << " " << flow5 << " " << flow6
1468  << " " << A1 << " " << A2 << " " << A3 << " " << A4 << " " << std::endl;
1469  }
1470 }
1471 
1472 
1473 void
1474 Pressure_flow_control_valve::SetValue(DataManager *pDM,
1475  VectorHandler& X , VectorHandler& XP,
1476  SimulationEntity::Hints *ph)
1477 {
1478  integer i = iGetFirstIndex();
1479 
1480  X.PutCoef(i+1, start);
1481  X.PutCoef(i+2, 0.);
1482  XP.PutCoef(i+1, 0.);
1483  XP.PutCoef(i+2, 0.);
1484 }
1485 
1486 /* Pressure_flow_control_valve - end */
1487 
1488 
1489 /* Pressure_valve - begin */
1490 
1491 Pressure_valve::Pressure_valve(unsigned int uL, const DofOwner* pDO,
1492  HydraulicFluid* hf,
1493  const PressureNode* p1, const PressureNode* p2,
1494  doublereal A_dia,
1495  doublereal mv,
1496  doublereal A_max, doublereal s_mx, doublereal K,
1497  doublereal F0, doublereal w,
1498  doublereal cs, doublereal cv, doublereal ca,
1499  flag fOut)
1500 : Elem(uL, fOut),
1501 HydraulicElem(uL, pDO, hf, fOut),
1502 pNode1(p1), pNode2(p2), area_diaf(A_dia), mass(mv),
1503 area_max(A_max),s_max(s_mx),
1504 Kappa(K), force0(F0), width(w),c_spost(cs), c_vel(cv), c_acc(ca)
1505 {
1506  ASSERT(pNode1 != NULL);
1507  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
1508  ASSERT(pNode2 != NULL);
1509  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
1510  ASSERT(A_dia > std::numeric_limits<doublereal>::epsilon());
1511  ASSERT(mv > std::numeric_limits<doublereal>::epsilon());
1512  ASSERT(A_max > std::numeric_limits<doublereal>::epsilon());
1513  ASSERT(K >= 0. );
1514  ASSERT(F0 >= 0.);
1515  ASSERT(w > std::numeric_limits<doublereal>::epsilon());
1516  ASSERT(s_mx>=0.); // corsa massima della valvola
1517 
1518 }
1519 
1520 Pressure_valve::~Pressure_valve(void)
1521 {
1522  NO_OP;
1523 }
1524 
1525 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
1526 HydraulicElem::Type Pressure_valve::GetHydraulicType(void) const
1527 {
1528  return HydraulicElem::PRESSURE_VALVE;
1529 }
1530 
1531 
1532 /* Contributo al file di restart */
1533 std::ostream& Pressure_valve::Restart(std::ostream& out) const
1534 {
1535  return out << "Pressure_valve not implemented yet!" << std::endl;
1536 }
1537 
1538 
1539 unsigned int Pressure_valve::iGetNumDof(void) const
1540 {
1541  return 2;
1542 }
1543 
1544 
1545 DofOrder::Order Pressure_valve::GetDofType(unsigned int i) const
1546 {
1547  ASSERT(i >= 0 && i <= 1);
1548  return DofOrder::DIFFERENTIAL;
1549 }
1550 
1551 
1552 void
1553 Pressure_valve::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
1554 {
1555  *piNumRows = 4;
1556  *piNumCols = 4;
1557 
1558 }
1559 
1560 
1561 VariableSubMatrixHandler&
1562 Pressure_valve::AssJac(VariableSubMatrixHandler& WorkMat,
1563  doublereal dCoef,
1564  const VectorHandler& XCurr,
1565  const VectorHandler& XPrimeCurr)
1566 {
1567  DEBUGCOUT("Entering Pressure_valve::AssJac()" << std::endl);
1568 
1569  FullSubMatrixHandler& WM = WorkMat.SetFull();
1570  WM.Resize(4, 4);
1571 
1572  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1573  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1574  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
1575  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
1576  integer iFirstIndex = iGetFirstIndex();
1577 
1578  WM.PutRowIndex(1, iNode1RowIndex);
1579  WM.PutRowIndex(2, iNode2RowIndex);
1580  WM.PutColIndex(1, iNode1ColIndex);
1581  WM.PutColIndex(2, iNode2ColIndex);
1582 
1583  WM.PutRowIndex(3, iFirstIndex+1);
1584  WM.PutColIndex(3, iFirstIndex+1);
1585  WM.PutRowIndex(4, iFirstIndex+2);
1586  WM.PutColIndex(4, iFirstIndex+2);
1587 
1588  doublereal p1 = pNode1->dGetX();
1589  doublereal p2 = pNode2->dGetX();
1590  /* p1 = XCurr(pNode1->iGetFirstRowIndex()+1); */
1591  doublereal density = HF->dGetDensity();
1592 
1593  s = XCurr(iFirstIndex+1); /* spostamento */
1594  v = XCurr(iFirstIndex+2); /* velocita' */
1595 
1596  doublereal Cd = .6;
1597  doublereal jumpPres = fabs(p1-p2);
1598 
1599  /* evito di dividere per un numero troppo piccolo */
1600  if (jumpPres < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1601  jumpPres = 1.e8*std::numeric_limits<doublereal>::epsilon();
1602  }
1603 
1604  doublereal Jac11 = 0.;
1605  doublereal Jac12 = 0.;
1606  doublereal Jac13 = 0.;
1607  doublereal Jac14 = 0.;
1608  doublereal Jac21 = 0.;
1609  doublereal Jac22 = 0.;
1610  doublereal Jac23 = 0.;
1611  doublereal Jac24 = 0.;
1612  doublereal Jac31 = 0.;
1613  doublereal Jac32 = 0.;
1614  doublereal Jac33 = 0.;
1615  doublereal Jac34 = 0.;
1616  doublereal Jac41 = 0.;
1617  doublereal Jac42 = 0.;
1618  doublereal Jac43 = 0.;
1619  doublereal Jac44 = 0.;
1620 
1621 #ifdef HYDR_DEVEL
1622  DEBUGCOUT("Valore di c1: " << c1 << std::endl);
1623  DEBUGCOUT("Valore di c2: " << c2 << std::endl);
1624  DEBUGCOUT("Valore di c3: " << c3 << std::endl);
1625 #endif
1626 
1627  Jac11 = -density*width*s*Cd/sqrt(2*jumpPres/density)/density;
1628  Jac12 = -Jac11;
1629  Jac13 = -density*width*Cd*copysign(sqrt(2*jumpPres/density), p1-p2)*dCoef;
1630  Jac14 = 0.;
1631 
1632  Jac21 = -Jac11;
1633  Jac22 = Jac11;
1634  Jac23 = -Jac13; /* density*width*Cd*copysign(sqrt(2*jumpPres/density),(p1-p2))*dCoef; */
1635  Jac24 = 0.;
1636 
1637  Jac31 = area_max;
1638  Jac32 = 0.;
1639  Jac33 = -Kappa*dCoef-c1*dCoef-cf1*dCoef;
1640  Jac34 = -mass
1641  -density*area_max*pow(area_max/(area_diaf*Cd), 2)*fabs(v)*dCoef
1642  -c2*dCoef-c3-cf2*dCoef-cf3;
1643 
1644  Jac41 = 0.;
1645  Jac42 = 0.;
1646  Jac43 = -1.;
1647  Jac44 = dCoef;
1648 
1649 #ifdef HYDR_DEVEL
1650  DEBUGCOUT("Jac smorzatore " << density*area_max*pow(area_max/(area_diaf*Cd), 2) << std::endl);
1651 #endif
1652 
1653  WM.PutCoef(1, 1, Jac11);
1654  WM.PutCoef(1, 2, Jac12);
1655  WM.PutCoef(1, 3, Jac13);
1656  WM.PutCoef(1, 4, Jac14);
1657  WM.PutCoef(2, 1, Jac21);
1658  WM.PutCoef(2, 2, Jac22);
1659  WM.PutCoef(2, 3, Jac23);
1660  WM.PutCoef(2, 4, Jac24);
1661  WM.PutCoef(3, 1, Jac31);
1662  WM.PutCoef(3, 2, Jac32);
1663  WM.PutCoef(3, 3, Jac33);
1664  WM.PutCoef(3, 4, Jac34);
1665  WM.PutCoef(4, 1, Jac41);
1666  WM.PutCoef(4, 2, Jac42);
1667  WM.PutCoef(4, 3, Jac43);
1668  WM.PutCoef(4, 4, Jac44);
1669 
1670  return WorkMat;
1671 }
1672 
1673 SubVectorHandler&
1674 Pressure_valve::AssRes(SubVectorHandler& WorkVec,
1675  doublereal dCoef,
1676  const VectorHandler& XCurr,
1677  const VectorHandler& XPrimeCurr)
1678 
1679 {
1680  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1681  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1682 
1683  doublereal p1 = pNode1->dGetX();
1684  doublereal p2 = pNode2->dGetX();
1685 
1686  WorkVec.Resize(4);
1687 
1688  integer iFirstIndex = iGetFirstIndex();
1689 
1690  s = XCurr(iFirstIndex+1); /* spostamento */
1691  v = XCurr(iFirstIndex+2); /* velocita' */
1692  sp = XPrimeCurr(iFirstIndex+1); /* velocita' */
1693  vp = XPrimeCurr(iFirstIndex+2); /* accelerazione */
1694  doublereal density = HF->dGetDensity();
1695 
1696  doublereal Cd = .6;
1697  doublereal jumpPres = fabs(p1-p2);
1698  doublereal Res_1 = 0.;
1699  doublereal Res_2 = 0.;
1700  doublereal Res_3 = 0.;
1701  doublereal Res_4 = 0.;
1702  doublereal x0 = (p1*area_max-force0)/Kappa; /* punto di equilibrio */
1703 
1704 #ifdef HYDR_DEVEL
1705  DEBUGCOUT("s: " << s << std::endl);
1706  DEBUGCOUT("sp: " << sp << std::endl);
1707  DEBUGCOUT("v: " << v << std::endl);
1708  DEBUGCOUT("vp: " << vp << std::endl);
1709 #endif
1710 
1711  if (s <= 0.) {
1712  c1 = c_spost;
1713 
1714  if (sp < 0.) {
1715  c2 = c_vel; /* ho v negativa devo smorzare */
1716  } else {
1717  c2 = 0.;
1718  }
1719 
1720  if (vp < 0.) {
1721  c3 = c_acc; /* ho acc negativa devo smorzare */
1722  } else {
1723  c3 = 0.;
1724  }
1725 
1726  c4 = 0.; /* (force0-p1*area_max); */
1727  } else {
1728  c1 = 0.;
1729  c2 = 0.;
1730  c3 = 0.;
1731  c4 = 0.;
1732  }
1733 
1734  if (s >= s_max) {
1735  cf1 = c_spost;
1736 
1737  if (sp > 0.) {
1738  cf2 = c_vel; /* ho v positiva devo smorzare */
1739  } else {
1740  cf2 = 0.;
1741  }
1742 
1743  if (vp > 0.) {
1744  cf3 = c_acc; /* ho acc positiva devo smorzare */
1745  } else {
1746  cf3 = 0.;
1747  }
1748 
1749  cf4 = 0.; /* (force0-p1*area_max); */
1750  } else {
1751  cf1 = 0.;
1752  cf2 = 0.;
1753  cf3 = 0.;
1754  cf4 = 0.;
1755  }
1756 
1757  Res_1 = density*width*s*Cd*copysign(sqrt(2.*jumpPres/density), p1-p2);
1758  Res_2 = -Res_1;
1759  Res_3 = mass*vp-p1*area_max
1760  +copysign(.5*density*area_max*pow(sp*area_max/(area_diaf*Cd), 2), sp)
1761  +Kappa*s+force0+c1*s+c2*sp+c3*vp+c4+cf1*(s-s_max)+cf2*sp+cf3*vp+cf4;
1762  Res_4 = sp-v;
1763 
1764 
1765  flow1 = -Res_1; /* per l'output */
1766  flow2 = -Res_2; /* per l'output */
1767 
1768 #ifdef HYDR_DEVEL
1769  DEBUGCOUT("width: " << width << std::endl);
1770  DEBUGCOUT("Cd: " << Cd << std::endl);
1771  DEBUGCOUT("jumpPres: " << jumpPres << std::endl);
1772  DEBUGCOUT("density: " << density << std::endl);
1773  DEBUGCOUT("p1: " << p1 << std::endl);
1774  DEBUGCOUT("p2: " << p2 << std::endl);
1775  DEBUGCOUT("s: " << s << std::endl);
1776  DEBUGCOUT("sp: " << sp << std::endl);
1777  DEBUGCOUT("v: " << v << std::endl);
1778  DEBUGCOUT("vp: " << vp << std::endl);
1779  DEBUGCOUT("area_max: " << area_max << std::endl);
1780  DEBUGCOUT("Kappa: " << Kappa << std::endl);
1781  DEBUGCOUT("force0: " << force0 << std::endl);
1782  DEBUGCOUT("mass: " << mass << std::endl);
1783  DEBUGCOUT("x0: " << x0 << std::endl);
1784  DEBUGCOUT("s_max: " << s_max << std::endl);
1785  DEBUGCOUT("Valore dello smorzatore: "
1786  << copysign(.5*density*area_max*pow(area_max/(area_diaf*Cd), 2), sp) << std::endl);
1787  DEBUGCOUT("c1: " << c1 << std::endl);
1788  DEBUGCOUT("c2: " << c2 << std::endl);
1789  DEBUGCOUT("c3: " << c3 << std::endl);
1790  DEBUGCOUT("c4: " << c4 << std::endl);
1791  DEBUGCOUT("cf1: " << cf1 << std::endl);
1792  DEBUGCOUT("cf2: " << cf2 << std::endl);
1793  DEBUGCOUT("cf3: " << cf3 << std::endl);
1794  DEBUGCOUT("cf4: " << cf4 << std::endl);
1795  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)" << std::endl);
1796  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
1797  DEBUGCOUT("-Res_2 (portata nodo2): " << -Res_2 << std::endl);
1798  DEBUGCOUT("Res_3: " << Res_3 << std::endl);
1799  DEBUGCOUT("Res_4: " << Res_4 << std::endl);
1800 #endif
1801 
1802  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
1803  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
1804  WorkVec.PutItem(3, iFirstIndex+1, Res_3);
1805  WorkVec.PutItem(4, iFirstIndex+2, Res_4);
1806 
1807  return WorkVec;
1808 }
1809 
1810 
1811 void Pressure_valve::Output(OutputHandler& OH) const
1812 {
1813  if (bToBeOutput()) {
1814  std::ostream& out = OH.Hydraulic();
1815  out << std::setw(8) << GetLabel()
1816  << " " << s << " " << v << " "<< vp
1817  << " " << flow1 << " " << flow2 << std::endl;
1818  }
1819 }
1820 
1821 
1822 void Pressure_valve::SetValue(DataManager *pDM,
1823  VectorHandler& X, VectorHandler& XP,
1824  SimulationEntity::Hints *ph)
1825 {
1826  integer i = iGetFirstIndex();
1827 
1828  X.PutCoef(i+1, 0.);
1829  X.PutCoef(i+2, 0.);
1830  XP.PutCoef(i+1, 0.);
1831  XP.PutCoef(i+2, 0.);
1832 }
1833 
1834 /* Pressure_valve - end */
1835 
1836 
1837 
1838  /* Flow_valve - begin */
1839 
1840 Flow_valve:: Flow_valve(unsigned int uL, const DofOwner* pDO,
1841  HydraulicFluid* hf,
1842  const PressureNode* p1, const PressureNode* p2,
1843  const PressureNode* p3,
1844  doublereal A_dia,
1845  doublereal mv, doublereal A_pipe, doublereal A_max,
1846  doublereal K, doublereal F0, doublereal w,
1847  doublereal s_mx,
1848  doublereal cs, doublereal cv, doublereal ca,
1849  flag fOut)
1850 : Elem(uL, fOut),
1851 HydraulicElem(uL, pDO, hf, fOut),
1852 pNode1(p1), pNode2(p2), pNode3(p3),
1853 area_diaf(A_dia), mass(mv), area_pipe(A_pipe), area_max(A_max),
1854 Kappa(K), force0(F0), width(w), s_max(s_mx), c_spost(cs), c_vel(cv), c_acc(ca)
1855 {
1856  ASSERT(pNode1 != NULL);
1857  ASSERT(pNode1->GetNodeType() == Node::HYDRAULIC);
1858  ASSERT(pNode2 != NULL);
1859  ASSERT(pNode2->GetNodeType() == Node::HYDRAULIC);
1860  ASSERT(pNode3 != NULL);
1861  ASSERT(pNode3->GetNodeType() == Node::HYDRAULIC);
1862  ASSERT(A_dia > std::numeric_limits<doublereal>::epsilon());
1863  ASSERT(mv > std::numeric_limits<doublereal>::epsilon());
1864  ASSERT(A_pipe > std::numeric_limits<doublereal>::epsilon());
1865  ASSERT(A_max > std::numeric_limits<doublereal>::epsilon());
1866  ASSERT(K > std::numeric_limits<doublereal>::epsilon() );
1867  ASSERT(F0 >= 0.);
1868  ASSERT(w > std::numeric_limits<doublereal>::epsilon());
1869  ASSERT(s_max >= 0.);
1870 
1871  h = .02; /* coefficiente di perdita di carico concentrata tra i nodi 1 e 2 (smorza il moto della valvola) */
1872 }
1873 
1874 
1875 Flow_valve::~Flow_valve(void)
1876 {
1877  NO_OP;
1878 }
1879 
1880 
1881 /* Tipo di elemento idraulico (usato solo per debug ecc.) */
1882 HydraulicElem::Type Flow_valve::GetHydraulicType(void) const
1883 {
1884  return HydraulicElem:: FLOW_VALVE;
1885 }
1886 
1887 
1888 /* Contributo al file di restart */
1889 std::ostream& Flow_valve::Restart(std::ostream& out) const
1890 {
1891  return out << " Flow_valve not implemented yet!" << std::endl;
1892 }
1893 
1894 
1895 unsigned int Flow_valve::iGetNumDof(void) const
1896 {
1897  return 2;
1898 }
1899 
1900 
1901 DofOrder::Order Flow_valve::GetDofType(unsigned int i) const
1902 {
1903  ASSERT(i >= 0 && i <= 1);
1904  return DofOrder::DIFFERENTIAL;
1905 }
1906 
1907 
1908 void
1909 Flow_valve::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
1910 {
1911  *piNumRows = 5;
1912  *piNumCols = 5;
1913 }
1914 
1915 VariableSubMatrixHandler&
1916 Flow_valve::AssJac(VariableSubMatrixHandler& WorkMat,
1917  doublereal dCoef,
1918  const VectorHandler& XCurr,
1919  const VectorHandler& XPrimeCurr)
1920 {
1921  DEBUGCOUT("Entering Flow_valve::AssJac()" << std::endl);
1922 
1923  FullSubMatrixHandler& WM = WorkMat.SetFull();
1924  WM.Resize(5, 5);
1925 
1926  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
1927  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
1928  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
1929  integer iNode1ColIndex = pNode1->iGetFirstColIndex()+1;
1930  integer iNode2ColIndex = pNode2->iGetFirstColIndex()+1;
1931  integer iNode3ColIndex = pNode3->iGetFirstColIndex()+1;
1932  integer iFirstIndex = iGetFirstIndex();
1933 
1934  WM.PutRowIndex(1, iNode1RowIndex);
1935  WM.PutRowIndex(2, iNode2RowIndex);
1936  WM.PutRowIndex(3, iNode3RowIndex);
1937  WM.PutColIndex(1, iNode1ColIndex);
1938  WM.PutColIndex(2, iNode2ColIndex);
1939  WM.PutColIndex(3, iNode3ColIndex);
1940 
1941  WM.PutRowIndex(4, iFirstIndex+1);
1942  WM.PutColIndex(4, iFirstIndex+1);
1943  WM.PutRowIndex(5, iFirstIndex+2);
1944  WM.PutColIndex(5, iFirstIndex+2);
1945 
1946  doublereal p1 = pNode1->dGetX();
1947  doublereal p2 = pNode2->dGetX();
1948  doublereal p3 = pNode3->dGetX();
1949  s = XCurr(iFirstIndex+1); /* spostamento */
1950  v = XCurr(iFirstIndex+2); /* velocita' */
1951 
1952  doublereal Cd = .6; /* verificare */
1953  doublereal jumpPres12 = fabs(p1-p2);
1954  doublereal jumpPres23 = fabs(p2-p3);
1955  doublereal jumpPres13 = fabs(p1-p3);
1956  doublereal density = HF->dGetDensity();
1957 
1958  /* evito di dividere per un numero troppo piccolo */
1959  if (jumpPres12 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1960  jumpPres12 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1961  }
1962  if (jumpPres23 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1963  jumpPres23 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1964  }
1965  if (jumpPres13 < 1.e8*std::numeric_limits<doublereal>::epsilon()) {
1966  jumpPres13 = 1.e8*std::numeric_limits<doublereal>::epsilon();
1967  }
1968 
1969  doublereal Jac11 = 0.;
1970  doublereal Jac12 = 0.;
1971  doublereal Jac13 = 0.;
1972  doublereal Jac14 = 0.;
1973  doublereal Jac15 = 0.;
1974  doublereal Jac21 = 0.;
1975  doublereal Jac22 = 0.;
1976  doublereal Jac23 = 0.;
1977  doublereal Jac24 = 0.;
1978  doublereal Jac25 = 0.;
1979  doublereal Jac31 = 0.;
1980  doublereal Jac32 = 0.;
1981  doublereal Jac33 = 0.;
1982  doublereal Jac34 = 0.;
1983  doublereal Jac35 = 0.;
1984  doublereal Jac41 = 0.;
1985  doublereal Jac42 = 0.;
1986  doublereal Jac43 = 0.;
1987  doublereal Jac44 = 0.;
1988  doublereal Jac45 = 0.;
1989  doublereal Jac51 = 0.;
1990  doublereal Jac52 = 0.;
1991  doublereal Jac53 = 0.;
1992  doublereal Jac54 = 0.;
1993  doublereal Jac55 = 0.;
1994 
1995 #ifdef HYDR_DEVEL
1996  DEBUGCOUT("Valore di p1: " << p1 << std::endl);
1997  DEBUGCOUT("Valore di p2: " << p2 << std::endl);
1998  DEBUGCOUT("Valore di p3: " << p3 << std::endl);
1999  DEBUGCOUT("Valore di s: " << s << std::endl);
2000 #endif
2001 
2002  doublereal s12 = sqrt(2.*jumpPres12/density);
2003  doublereal s13 = sqrt(2.*jumpPres13/density);
2004 
2005  Jac11 = -width*s*Cd/s12-area_diaf*Cd/s13;
2006  Jac12 = width*s*Cd/s12;
2007  Jac13 = area_diaf*Cd/s13;
2008  Jac14 = -density*width*Cd*copysign(s12, p1-p2)*dCoef;
2009  Jac15 = 0.;
2010 
2011  Jac21 = width*s*Cd/s12;
2012  Jac22 = -width*s*Cd/s12;
2013  Jac23 = 0.;
2014  Jac24 = density*width*Cd*copysign(s12, p1-p2)*dCoef;
2015  Jac25 = 0.;
2016 
2017  Jac31 = area_diaf*Cd/s13;
2018  Jac32 = 0.;
2019  Jac33 = -area_diaf*Cd/s13;
2020  Jac34 = 0.;
2021  Jac35 = 0.;
2022 
2023  Jac41 = area_max;
2024  Jac42 = 0.;
2025  Jac43 =-area_max;
2026  Jac44 = -Kappa*dCoef-c1*dCoef-cf1*dCoef;
2027 
2028 
2029  doublereal Jac45old1;
2030  doublereal Jac45old2;
2031  Jac45old1 = -mass-c2*dCoef-c3-cf2*dCoef-cf3
2032  -density*area_max*pow(area_max/(area_diaf*Cd), 2.)*fabs(v)*dCoef;
2033  Jac45 = -mass-c2*dCoef-c3-cf2*dCoef-cf3
2034  -h*density*pow(area_max/area_pipe, 2.)*fabs(v)*dCoef;
2035  Jac45old2 = -mass-c2*dCoef-c3-cf2*dCoef-cf3-44.4*fabs(v)*dCoef;
2036 
2037  Jac51 = 0.;
2038  Jac52 = 0.;
2039  Jac53 = 0.;
2040  Jac54 = -1.;
2041  Jac55 = dCoef;
2042 
2043 #ifdef HYDR_DEVEL
2044  DEBUGCOUT("Jac smorzatore "
2045  << density*area_max*pow(area_max/(area_diaf*Cd), 2) << std::endl);
2046  DEBUGCOUT("Jac smorzatore " << h*density << std::endl);
2047  DEBUGCOUT("Jac smorzatore Jac45old1 " << Jac45old1 << std::endl);
2048  DEBUGCOUT("Jac smorzatore Jac45old2 " << Jac45old2 << std::endl);
2049  DEBUGCOUT("Jac smorzatore Jac45 " << Jac45 << std::endl);
2050 
2051  DEBUGCOUT("Jac11: " << Jac11 << std::endl);
2052  DEBUGCOUT("Jac12: " << Jac12 << std::endl);
2053  DEBUGCOUT("Jac13: " << Jac13 << std::endl);
2054  DEBUGCOUT("Jac14: " << Jac14 << std::endl);
2055  DEBUGCOUT("Jac15: " << Jac15 << std::endl);
2056  DEBUGCOUT("Jac21: " << Jac21 << std::endl);
2057  DEBUGCOUT("Jac22: " << Jac22 << std::endl);
2058  DEBUGCOUT("Jac23: " << Jac23 << std::endl);
2059  DEBUGCOUT("Jac24: " << Jac24 << std::endl);
2060  DEBUGCOUT("Jac25: " << Jac25 << std::endl);
2061  DEBUGCOUT("Jac31: " << Jac31 << std::endl);
2062  DEBUGCOUT("Jac32: " << Jac32 << std::endl);
2063  DEBUGCOUT("Jac33: " << Jac33 << std::endl);
2064  DEBUGCOUT("Jac34: " << Jac34 << std::endl);
2065  DEBUGCOUT("Jac35: " << Jac35 << std::endl);
2066  DEBUGCOUT("Jac41: " << Jac41 << std::endl);
2067  DEBUGCOUT("Jac42: " << Jac42 << std::endl);
2068  DEBUGCOUT("Jac43: " << Jac43 << std::endl);
2069  DEBUGCOUT("Jac44: " << Jac44 << std::endl);
2070  DEBUGCOUT("Jac45: " << Jac45 << std::endl);
2071  DEBUGCOUT("Jac51: " << Jac51 << std::endl);
2072  DEBUGCOUT("Jac52: " << Jac52 << std::endl);
2073  DEBUGCOUT("Jac53: " << Jac53 << std::endl);
2074  DEBUGCOUT("Jac54: " << Jac54 << std::endl);
2075  DEBUGCOUT("Jac55: " << Jac55 << std::endl);
2076 #endif
2077 
2078  WM.PutCoef(1, 1, Jac11);
2079  WM.PutCoef(1, 2, Jac12);
2080  WM.PutCoef(1, 3, Jac13);
2081  WM.PutCoef(1, 4, Jac14);
2082  WM.PutCoef(1, 5, Jac15);
2083  WM.PutCoef(2, 1, Jac21);
2084  WM.PutCoef(2, 2, Jac22);
2085  WM.PutCoef(2, 3, Jac23);
2086  WM.PutCoef(2, 4, Jac24);
2087  WM.PutCoef(2, 5, Jac25);
2088  WM.PutCoef(3, 1, Jac31);
2089  WM.PutCoef(3, 2, Jac32);
2090  WM.PutCoef(3, 3, Jac33);
2091  WM.PutCoef(3, 4, Jac34);
2092  WM.PutCoef(3, 5, Jac35);
2093  WM.PutCoef(4, 1, Jac41);
2094  WM.PutCoef(4, 2, Jac42);
2095  WM.PutCoef(4, 3, Jac43);
2096  WM.PutCoef(4, 4, Jac44);
2097  WM.PutCoef(4, 5, Jac45);
2098  WM.PutCoef(5, 1, Jac51);
2099  WM.PutCoef(5, 2, Jac52);
2100  WM.PutCoef(5, 3, Jac53);
2101  WM.PutCoef(5, 4, Jac54);
2102  WM.PutCoef(5, 5, Jac55);
2103 
2104  return WorkMat;
2105 }
2106 
2107 SubVectorHandler&
2108 Flow_valve::AssRes(SubVectorHandler& WorkVec,
2109  doublereal dCoef,
2110  const VectorHandler& XCurr,
2111  const VectorHandler& XPrimeCurr)
2112 {
2113  DEBUGCOUT("Entering Flow_valve::AssRes()" << std::endl);
2114 
2115  WorkVec.Resize(5);
2116 
2117  integer iNode1RowIndex = pNode1->iGetFirstRowIndex()+1;
2118  integer iNode2RowIndex = pNode2->iGetFirstRowIndex()+1;
2119  integer iNode3RowIndex = pNode3->iGetFirstRowIndex()+1;
2120 
2121  doublereal p1 = pNode1->dGetX();
2122  doublereal p2 = pNode2->dGetX();
2123  doublereal p3 = pNode3->dGetX();
2124  integer iFirstIndex = iGetFirstIndex();
2125 
2126  s = XCurr(iFirstIndex+1); /* spostamento */
2127  v = XCurr(iFirstIndex+2); /* velocita' */
2128  sp = XPrimeCurr(iFirstIndex+1); /* velocita' */
2129  vp = XPrimeCurr(iFirstIndex+2); /* accelerazione */
2130 
2131  doublereal Cd = .6;
2132  doublereal jumpPres12 = fabs(p1-p2);
2133  doublereal jumpPres23 = fabs(p2-p3);
2134  doublereal jumpPres13 = fabs(p1-p3);
2135 
2136  doublereal Res_1 = 0.;
2137  doublereal Res_2 = 0.;
2138  doublereal Res_3 = 0.;
2139  doublereal Res_4 = 0.;
2140  doublereal Res_5 = 0.;
2141  doublereal density = HF->dGetDensity();
2142 
2143  doublereal x0 = ((p1-p3)*area_max-force0)/Kappa;
2144 
2145  if (s <= 0.) {
2146  c1 = c_spost;
2147 
2148  if (sp < 0.) {
2149  c2 = c_vel; /* ho v negativa devo smorzare */
2150  } else {
2151  c2 = 0.;
2152  }
2153 
2154  if (vp < 0.) {
2155  c3 = c_acc; /* ho acc negativa devo smorzare */
2156  } else {
2157  c3 = 0.;
2158  }
2159 
2160  c4 = 0.; /* (force0-p1*area_max); */
2161  } else {
2162  c1 = 0.;
2163  c2 = 0.;
2164  c3 = 0.;
2165  c4 = 0.;
2166  }
2167 
2168  if (s >= s_max) {
2169  cf1 = c_spost;
2170 
2171  if (sp > 0.) {
2172  cf2 = c_vel; /* ho v positiva devo smorzare */
2173  } else {
2174  cf2 = 0.;
2175  }
2176 
2177  if (vp > 0.) {
2178  cf3 = c_acc; /* ho acc positiva devo smorzare */
2179  } else {
2180  cf3 = 0.;
2181  }
2182 
2183  cf4 = 0.; /* (force0-p1*area_max); */
2184  } else {
2185  cf1 = 0.;
2186  cf2 = 0.;
2187  cf3 = 0.;
2188  cf4 = 0.;
2189  }
2190 
2191  doublereal s12 = sqrt(2.*jumpPres12/density);
2192  doublereal s13 = sqrt(2.*jumpPres13/density);
2193 
2194  Res_1 = density*width*s*Cd*copysign(s12, p1-p2)+
2195  density*area_diaf*Cd*copysign(s13, p1-p3);
2196  Res_2 = -density*width*s*Cd*copysign(s12, p1-p2);
2197  Res_3 = -density*area_diaf*Cd*copysign(s13, p1-p3);
2198 
2199  doublereal Res_4old;
2200  Res_4old = mass*vp-area_max*p1+force0+Kappa*s+area_max*p3+c1*s+c2*sp+c3*vp
2201  +c4+cf1*(s-s_max)+cf2*sp+cf3*vp+cf4
2202  +copysign(.5*density*area_max*pow(sp*area_max/(area_diaf*Cd), 2), sp);
2203 
2204  Res_4 = mass*vp-area_max*p1+force0+Kappa*s+area_max*p3+c1*s+c2*sp+c3*vp
2205  +c4+cf1*(s-s_max)+cf2*sp+cf3*vp+cf4
2206  +copysign(.5*h*density*pow(sp*area_max/area_pipe, 2), sp);
2207 
2208  Res_5 =sp-v;
2209 
2210  flow1=-Res_1; /* per l'output */
2211  flow2=-Res_2; /* per l'output */
2212  flow3=-Res_3; /* per l'output */
2213 
2214 #ifdef HYDR_DEVEL
2215  DEBUGCOUT("Res_4: " << Res_4 << std::endl);
2216  DEBUGCOUT("Res_4old: " << Res_4old << std::endl);
2217  DEBUGCOUT("smorazatore: "
2218  << copysign(.5*density*area_max*pow(sp*area_max/(area_diaf*Cd), 2), sp) << std::endl);
2219  DEBUGCOUT("smorzatoreold: "
2220  << copysign(.5*h*density*pow(sp*area_max/area_pipe, 2), sp) << std::endl);
2221 
2222  DEBUGCOUT("width: " << width << std::endl);
2223  DEBUGCOUT("Cd: " << Cd << std::endl);
2224  DEBUGCOUT("jumpPres12: " << jumpPres12 << std::endl);
2225  DEBUGCOUT("jumpPres23: " << jumpPres23 << std::endl);
2226  DEBUGCOUT("jumpPres13: " << jumpPres13 << std::endl);
2227  DEBUGCOUT("density: " << density << std::endl);
2228  DEBUGCOUT("p1: " << p1 << std::endl);
2229  DEBUGCOUT("p2: " << p2 << std::endl);
2230  DEBUGCOUT("p3: " << p3 << std::endl);
2231  DEBUGCOUT("s: " << s << std::endl);
2232  DEBUGCOUT("sp: " << sp << std::endl);
2233  DEBUGCOUT("v: " << v << std::endl);
2234  DEBUGCOUT("vp: " << vp << std::endl);
2235  DEBUGCOUT("area_max: " << area_max << std::endl);
2236  DEBUGCOUT("Kappa: " << Kappa << std::endl);
2237  DEBUGCOUT("force0: " << force0 << std::endl);
2238  DEBUGCOUT("mass: " << mass << std::endl);
2239  DEBUGCOUT("s_max: " << s_max << std::endl);
2240  DEBUGCOUT("x0: " << x0 << std::endl);
2241  DEBUGCOUT("c1: " << c1 << std::endl);
2242  DEBUGCOUT("c2: " << c2 << std::endl);
2243  DEBUGCOUT("c3: " << c3 << std::endl);
2244  DEBUGCOUT("c4: " << c4 << std::endl);
2245  DEBUGCOUT("cf1: " << cf1 << std::endl);
2246  DEBUGCOUT("cf2: " << cf2 << std::endl);
2247  DEBUGCOUT("cf3: " << cf3 << std::endl);
2248  DEBUGCOUT("cf4: " << cf4 << std::endl);
2249  DEBUGCOUT("PORTATE AI VARI NODI (positive se entranti)" << std::endl);
2250  DEBUGCOUT("-Res_1 (portata nodo1): " << -Res_1 << std::endl);
2251  DEBUGCOUT("-Res_2 (portata nodo2): " << -Res_2 << std::endl);
2252  DEBUGCOUT("-Res_3:(portata nodo3): " << -Res_3 << std::endl);
2253  DEBUGCOUT("-Res_4: " << -Res_4 << std::endl);
2254  DEBUGCOUT("-Res_5: " << -Res_5 << std::endl);
2255 #endif
2256 
2257  WorkVec.PutItem(1, iNode1RowIndex, Res_1);
2258  WorkVec.PutItem(2, iNode2RowIndex, Res_2);
2259  WorkVec.PutItem(3, iNode3RowIndex, Res_3);
2260  WorkVec.PutItem(4, iFirstIndex+1, Res_4);
2261  WorkVec.PutItem(5, iFirstIndex+2, Res_5);
2262 
2263  return WorkVec;
2264 }
2265 
2266 void Flow_valve::Output(OutputHandler& OH) const
2267 {
2268  if (bToBeOutput()) {
2269  std::ostream& out = OH.Hydraulic();
2270  out << std::setw(8) << GetLabel()
2271  << " " << s << " " << v << " "<< vp
2272  << " " << flow1 << " "<< flow2 << " "<< flow3 << std::endl;
2273  }
2274 }
2275 
2276 void Flow_valve::SetValue(DataManager *pDM,
2277  VectorHandler& X, VectorHandler& XP,
2278  SimulationEntity::Hints *ph)
2279 {
2280  integer i = iGetFirstIndex();
2281  X.PutCoef(i+1, 0.);
2282  X.PutCoef(i+2, 0.);
2283  XP.PutCoef(i+1, 0.);
2284  XP.PutCoef(i+2, 0.);
2285 }
2286 
2287 /* Flow_valve - end */
Flow_valve::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:1882
Dynamic_control_valve::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:705
Dynamic_control_valve::sp
doublereal sp
Definition: valve.h:258
Pressure_flow_control_valve::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:1154
Dynamic_control_valve::deltaP
doublereal deltaP
Definition: valve.h:266
Dynamic_control_valve::A2
doublereal A2
Definition: valve.h:247
Pressure_flow_control_valve::s
doublereal s
Definition: valve.h:359
Control_valve::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:112
Flow_valve::s
doublereal s
Definition: valve.h:542
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Pressure_flow_control_valve::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:1134
M_PI
#define M_PI
Definition: gradienttest.cc:67
Pressure_valve::v
doublereal v
Definition: valve.h:455
Control_valve2::area_max
doublereal area_max
Definition: valve.h:164
Dynamic_control_valve::pNode1
const PressureNode * pNode1
Definition: valve.h:228
Pressure_valve::width
doublereal width
Definition: valve.h:440
Control_valve2::density
doublereal density
Definition: valve.h:160
flag
long int flag
Definition: mbdyn.h:43
Pressure_flow_control_valve::deltaP
doublereal deltaP
Definition: valve.h:369
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
Pressure_flow_control_valve::s_max
doublereal s_max
Definition: valve.h:344
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
Dynamic_control_valve::start
doublereal start
Definition: valve.h:232
Flow_valve::c1
doublereal c1
Definition: valve.h:532
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Pressure_valve::sp
doublereal sp
Definition: valve.h:456
Flow_valve::cf2
doublereal cf2
Definition: valve.h:537
Flow_valve::pNode1
const PressureNode * pNode1
Definition: valve.h:518
Control_valve2::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:385
Dynamic_control_valve::cf2
doublereal cf2
Definition: valve.h:264
Pressure_flow_control_valve::sp
doublereal sp
Definition: valve.h:361
Dynamic_control_valve::vp
doublereal vp
Definition: valve.h:259
Pressure_flow_control_valve::A1
doublereal A1
Definition: valve.h:346
Pressure_flow_control_valve::cf1
doublereal cf1
Definition: valve.h:366
Pressure_flow_control_valve::v
doublereal v
Definition: valve.h:360
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
Flow_valve::mass
doublereal mass
Definition: valve.h:522
Flow_valve::pNode2
const PressureNode * pNode2
Definition: valve.h:519
Flow_valve::c2
doublereal c2
Definition: valve.h:533
Pressure_flow_control_valve::~Pressure_flow_control_valve
~Pressure_flow_control_valve(void)
Definition: valve.cc:1117
Pressure_flow_control_valve::cf2
doublereal cf2
Definition: valve.h:367
Flow_valve::flow2
doublereal flow2
Definition: valve.h:547
Pressure_flow_control_valve::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:1129
Control_valve::A4min
doublereal A4min
Definition: valve.h:66
Pressure_flow_control_valve::pNode1
const PressureNode * pNode1
Definition: valve.h:326
Control_valve::pNode2
const PressureNode * pNode2
Definition: valve.h:47
Flow_valve::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:1901
Pressure_valve::c4
doublereal c4
Definition: valve.h:448
Pressure_flow_control_valve::loss_area
doublereal loss_area
Definition: valve.h:340
Pressure_valve::pNode2
const PressureNode * pNode2
Definition: valve.h:434
Pressure_valve::c_acc
doublereal c_acc
Definition: valve.h:444
Flow_valve::v
doublereal v
Definition: valve.h:544
Dynamic_control_valve::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:1042
Control_valve2::Cd
doublereal Cd
Definition: valve.h:163
FullSubMatrixHandler::Resize
void Resize(integer iNewRow, integer iNewCol)
Definition: submat.cc:138
Dynamic_control_valve::flow2
doublereal flow2
Definition: valve.h:253
Dynamic_control_valve::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:864
Flow_valve::~Flow_valve
~Flow_valve(void)
Definition: valve.cc:1875
Control_valve2::f
doublereal f[LAST_N]
Definition: valve.h:168
Pressure_valve::Kappa
doublereal Kappa
Definition: valve.h:438
Dynamic_control_valve::~Dynamic_control_valve
~Dynamic_control_valve(void)
Definition: valve.cc:675
Control_valve2::Q34
Definition: valve.h:145
Pressure_valve::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:1553
Pressure_valve::flow2
doublereal flow2
Definition: valve.h:459
Control_valve::A2min
doublereal A2min
Definition: valve.h:64
Pressure_flow_control_valve::pNode4
const PressureNode * pNode4
Definition: valve.h:329
Pressure_valve::cf3
doublereal cf3
Definition: valve.h:451
Control_valve::Control_valve
Control_valve(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, const PressureNode *p3, const PressureNode *p4, doublereal A_max, doublereal Loss_a, const DriveCaller *pDC, flag fOut)
Definition: valve.cc:46
Dynamic_control_valve::pNode4
const PressureNode * pNode4
Definition: valve.h:231
Flow_valve::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:2266
Dynamic_control_valve::c_acc
doublereal c_acc
Definition: valve.h:235
Pressure_flow_control_valve::pNode3
const PressureNode * pNode3
Definition: valve.h:328
Control_valve::pNode4
const PressureNode * pNode4
Definition: valve.h:49
Dynamic_control_valve::valve_diameter
doublereal valve_diameter
Definition: valve.h:241
FullSubMatrixHandler::PutCoef
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
Pressure_flow_control_valve::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:1460
Flow_valve::sp
doublereal sp
Definition: valve.h:543
Pressure_flow_control_valve::pNode6
const PressureNode * pNode6
Definition: valve.h:331
Pressure_flow_control_valve::A2
doublereal A2
Definition: valve.h:347
Control_valve2::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:404
Control_valve::A1
doublereal A1
Definition: valve.h:59
Force
Definition: force.h:46
Control_valve::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:104
Dynamic_control_valve::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:697
Pressure_flow_control_valve::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:1123
Pressure_valve::cf4
doublereal cf4
Definition: valve.h:452
Control_valve::pNode3
const PressureNode * pNode3
Definition: valve.h:48
Pressure_flow_control_valve::width
doublereal width
Definition: valve.h:339
Dynamic_control_valve::A1
doublereal A1
Definition: valve.h:246
Flow_valve::flow1
doublereal flow1
Definition: valve.h:546
NO_OP
#define NO_OP
Definition: myassert.h:74
Control_valve2::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:542
Pressure_valve::c1
doublereal c1
Definition: valve.h:445
Pressure_valve::area_max
doublereal area_max
Definition: valve.h:437
Dynamic_control_valve::Cd
doublereal Cd
Definition: valve.h:238
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
Control_valve2::N4
Definition: valve.h:140
Control_valve2::pNode
const PressureNode * pNode[LAST_N]
Definition: valve.h:157
Pressure_flow_control_valve::vp
doublereal vp
Definition: valve.h:362
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
Pressure_flow_control_valve::c3
doublereal c3
Definition: valve.h:365
SubVectorHandler::PutItem
virtual void PutItem(integer iSubRow, integer iRow, const doublereal &dCoef)
Definition: submat.h:1445
Pressure_flow_control_valve::valve_density
doublereal valve_density
Definition: valve.h:342
Control_valve2::Prepare
void Prepare(void)
Definition: valve.cc:488
Control_valve2::area_min
doublereal area_min
Definition: valve.h:166
Flow_valve::c_acc
doublereal c_acc
Definition: valve.h:531
Flow_valve::cf3
doublereal cf3
Definition: valve.h:538
Pressure_flow_control_valve::flow2
doublereal flow2
Definition: valve.h:353
Control_valve2::Q12
Definition: valve.h:144
Control_valve::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:120
Flow_valve::h
doublereal h
Definition: valve.h:541
Control_valve::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:94
Pressure_valve::cf1
doublereal cf1
Definition: valve.h:449
Control_valve2::Stato
doublereal Stato
Definition: valve.h:162
Pressure_flow_control_valve::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:1147
Dynamic_control_valve::c_spost
doublereal c_spost
Definition: valve.h:233
Flow_valve::cf1
doublereal cf1
Definition: valve.h:536
Dynamic_control_valve::pNode2
const PressureNode * pNode2
Definition: valve.h:229
Control_valve::Cd
doublereal Cd
Definition: valve.h:51
Control_valve::area_max
doublereal area_max
Definition: valve.h:52
Pressure_valve::Pressure_valve
Pressure_valve(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, doublereal A_dia, doublereal mv, doublereal A_max, doublereal s_mx, doublereal K, doublereal F0, doublereal w, doublereal cs, doublereal cv, doublereal ca, flag fOut)
Definition: valve.cc:1491
Pressure_flow_control_valve::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:1139
Pressure_flow_control_valve::valve_diameter
doublereal valve_diameter
Definition: valve.h:341
HydraulicFluid::dGetDensity
virtual doublereal dGetDensity(void) const =0
Flow_valve::force0
doublereal force0
Definition: valve.h:526
Control_valve2
Definition: valve.h:120
ScalarAlgebraicNode::dGetX
virtual const doublereal & dGetX(void) const
Definition: node.h:492
Pressure_flow_control_valve::Pressure_flow_control_valve
Pressure_flow_control_valve(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, const PressureNode *p3, const PressureNode *p4, const PressureNode *p5, const PressureNode *p6, const DriveCaller *pDC, doublereal s0, doublereal s_mx, doublereal W, doublereal Loss_A, doublereal Valve_d, doublereal Valve_rho, doublereal cs, doublereal cv, doublereal ca, flag fOut)
Definition: valve.cc:1072
Dynamic_control_valve::flow3
doublereal flow3
Definition: valve.h:254
Flow_valve::c4
doublereal c4
Definition: valve.h:535
Control_valve2::A
doublereal A[LAST_Q]
Definition: valve.h:169
Control_valve::flow4
doublereal flow4
Definition: valve.h:58
Pressure_valve::c_vel
doublereal c_vel
Definition: valve.h:443
Pressure_valve::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: valve.cc:1822
Pressure_flow_control_valve::pNode5
const PressureNode * pNode5
Definition: valve.h:330
Dynamic_control_valve::cf3
doublereal cf3
Definition: valve.h:265
Flow_valve::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:2108
Flow_valve::flow3
doublereal flow3
Definition: valve.h:548
Control_valve2::loss_area
doublereal loss_area
Definition: valve.h:165
Pressure_flow_control_valve::A3
doublereal A3
Definition: valve.h:348
Pressure_flow_control_valve::flow4
doublereal flow4
Definition: valve.h:355
Pressure_flow_control_valve::flow3
doublereal flow3
Definition: valve.h:354
Dynamic_control_valve::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:681
Flow_valve::width
doublereal width
Definition: valve.h:527
Control_valve::flow1
doublereal flow1
Definition: valve.h:55
HydraulicElem::HF
HydraulicFluid * HF
Definition: preselem.h:79
Pressure_flow_control_valve::start
doublereal start
Definition: valve.h:332
Control_valve::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:216
Dynamic_control_valve::c2
doublereal c2
Definition: valve.h:261
Pressure_valve::c3
doublereal c3
Definition: valve.h:447
Pressure_flow_control_valve::flow6
doublereal flow6
Definition: valve.h:357
Pressure_valve::s
doublereal s
Definition: valve.h:454
Pressure_valve::mass
doublereal mass
Definition: valve.h:436
Pressure_flow_control_valve::valve_area
doublereal valve_area
Definition: valve.h:343
Flow_valve::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:1909
Dynamic_control_valve::cf1
doublereal cf1
Definition: valve.h:263
Control_valve2::dp
doublereal dp[LAST_Q]
Definition: valve.h:159
Control_valve::pNode1
const PressureNode * pNode1
Definition: valve.h:46
Pressure_valve::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:1526
Flow_valve::cf4
doublereal cf4
Definition: valve.h:539
Pressure_flow_control_valve::Cd
doublereal Cd
Definition: valve.h:338
Pressure_valve::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:1562
Control_valve2::N2
Definition: valve.h:138
copysign
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
Control_valve2::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: valve.cc:397
Pressure_valve::~Pressure_valve
~Pressure_valve(void)
Definition: valve.cc:1520
Pressure_valve::cf2
doublereal cf2
Definition: valve.h:450
Control_valve::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:314
Pressure_flow_control_valve::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:1312
Control_valve2::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:592
Node::iGetFirstRowIndex
virtual integer iGetFirstRowIndex(void) const
Definition: node.cc:82
Pressure_valve::force0
doublereal force0
Definition: valve.h:439
Control_valve2::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: valve.cc:604
Dynamic_control_valve::valve_density
doublereal valve_density
Definition: valve.h:242
Pressure_flow_control_valve::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: valve.cc:1474
Control_valve::A2
doublereal A2
Definition: valve.h:60
Dynamic_control_valve::flow1
doublereal flow1
Definition: valve.h:252
Flow_valve::vp
doublereal vp
Definition: valve.h:545
Flow_valve::pNode3
const PressureNode * pNode3
Definition: valve.h:520
Pressure_valve::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:1533
Dynamic_control_valve::s
doublereal s
Definition: valve.h:256
Dynamic_control_valve::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: valve.cc:1055
Pressure_flow_control_valve::flow5
doublereal flow5
Definition: valve.h:356
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
Control_valve2::Q14
Definition: valve.h:149
Control_valve2::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:379
Control_valve2::N3
Definition: valve.h:139
Dynamic_control_valve::Dynamic_control_valve
Dynamic_control_valve(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, const PressureNode *p3, const PressureNode *p4, const DriveCaller *pDC, doublereal s0, doublereal s_mx, doublereal W, doublereal Loss_A, doublereal Valve_d, doublereal Valve_rho, doublereal cs, doublereal cv, doublereal ca, flag fOut)
Definition: valve.cc:628
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
Pressure_valve::c2
doublereal c2
Definition: valve.h:446
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
Dynamic_control_valve::A4
doublereal A4
Definition: valve.h:249
DriveOwner::pGetDriveCaller
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
Control_valve::A3min
doublereal A3min
Definition: valve.h:65
Mass
Definition: body.h:45
Control_valve2::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:390
Control_valve2::~Control_valve2
~Control_valve2(void)
Definition: valve.cc:365
Pressure_valve::flow1
doublereal flow1
Definition: valve.h:458
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
Dynamic_control_valve::s_max
doublereal s_max
Definition: valve.h:244
Dynamic_control_valve::loss_area
doublereal loss_area
Definition: valve.h:240
Dynamic_control_valve::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:712
DriveCaller::dGet
virtual doublereal dGet(const doublereal &dVar) const =0
Dynamic_control_valve::c1
doublereal c1
Definition: valve.h:260
Pressure_valve::c_spost
doublereal c_spost
Definition: valve.h:442
Control_valve2::Q23
Definition: valve.h:150
Pressure_valve::GetDofType
virtual DofOrder::Order GetDofType(unsigned int i) const
Definition: valve.cc:1545
Elem
Definition: elem.h:75
Pressure_valve::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:1539
Pressure_valve::pNode1
const PressureNode * pNode1
Definition: valve.h:433
Dynamic_control_valve::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:687
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
Dynamic_control_valve::A3
doublereal A3
Definition: valve.h:248
Dynamic_control_valve::v
doublereal v
Definition: valve.h:257
Control_valve2::Control_valve2
Control_valve2(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, const PressureNode *p3, const PressureNode *p4, doublereal A_max, doublereal Loss_a, const DriveCaller *pDC, flag fOut)
Definition: valve.cc:330
Dynamic_control_valve::width
doublereal width
Definition: valve.h:239
Control_valve2::q
doublereal q[LAST_Q]
Definition: valve.h:158
Dynamic_control_valve::pNode3
const PressureNode * pNode3
Definition: valve.h:230
Pressure_flow_control_valve::c1
doublereal c1
Definition: valve.h:363
Flow_valve::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:1895
PressureNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: presnode.h:54
Control_valve2::N1
Definition: valve.h:137
Control_valve::loss_area
doublereal loss_area
Definition: valve.h:53
Control_valve::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:88
Pressure_valve::area_diaf
doublereal area_diaf
Definition: valve.h:435
Flow_valve::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: valve.cc:2276
Dynamic_control_valve::flow4
doublereal flow4
Definition: valve.h:255
Control_valve2::Q13
Definition: valve.h:146
Control_valve2::LAST_Q
Definition: valve.h:151
a
static const doublereal a
Definition: hfluid_.h:289
Flow_valve::area_pipe
doublereal area_pipe
Definition: valve.h:523
OutputHandler::Hydraulic
std::ostream & Hydraulic(void) const
Definition: output.h:492
Flow_valve::c_spost
doublereal c_spost
Definition: valve.h:529
Flow_valve::s_max
doublereal s_max
Definition: valve.h:528
Control_valve::A1min
doublereal A1min
Definition: valve.h:63
Flow_valve::area_max
doublereal area_max
Definition: valve.h:524
Pressure_flow_control_valve::flow1
doublereal flow1
Definition: valve.h:352
Control_valve::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:99
Pressure_flow_control_valve::cf3
doublereal cf3
Definition: valve.h:368
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
Flow_valve::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: valve.cc:1889
Pressure_flow_control_valve::A4
doublereal A4
Definition: valve.h:349
doublereal
double doublereal
Definition: colamd.c:52
Control_valve2::LAST_N
Definition: valve.h:141
Control_valve::~Control_valve
~Control_valve(void)
Definition: valve.cc:82
integer
long int integer
Definition: colamd.c:51
Dynamic_control_valve::c3
doublereal c3
Definition: valve.h:262
Control_valve::Stato
doublereal Stato
Definition: valve.h:50
Flow_valve::Kappa
doublereal Kappa
Definition: valve.h:525
Dynamic_control_valve::c_vel
doublereal c_vel
Definition: valve.h:234
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Control_valve::flow2
doublereal flow2
Definition: valve.h:56
Flow_valve::Flow_valve
Flow_valve(unsigned int uL, const DofOwner *pD, HydraulicFluid *hf, const PressureNode *p1, const PressureNode *p2, const PressureNode *p3, doublereal A_dia, doublereal mv, doublereal A_pipe, doublereal A_max, doublereal K, doublereal F0, doublereal w, doublereal s_mx, doublereal cs, doublereal cv, doublereal ca, flag fOut)
Definition: valve.cc:1840
Pressure_flow_control_valve::pNode2
const PressureNode * pNode2
Definition: valve.h:327
Flow_valve::c_vel
doublereal c_vel
Definition: valve.h:530
Pressure_valve::vp
doublereal vp
Definition: valve.h:457
Pressure_valve::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:1674
Control_valve::A3
doublereal A3
Definition: valve.h:61
Control_valve2::Q24
Definition: valve.h:147
Pressure_flow_control_valve::c2
doublereal c2
Definition: valve.h:364
Flow_valve::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: valve.cc:1916
VectorHandler::Resize
virtual void Resize(integer iNewSize)=0
Control_valve::A4
doublereal A4
Definition: valve.h:62
Pressure_valve::Output
virtual void Output(OutputHandler &OH) const
Definition: valve.cc:1811
Flow_valve::area_diaf
doublereal area_diaf
Definition: valve.h:521
Flow_valve::c3
doublereal c3
Definition: valve.h:534
Dynamic_control_valve::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: valve.cc:692
Control_valve::flow3
doublereal flow3
Definition: valve.h:57
Pressure_valve::s_max
doublereal s_max
Definition: valve.h:441
Control_valve2::GetHydraulicType
virtual HydraulicElem::Type GetHydraulicType(void) const
Definition: valve.cc:372
Node::iGetFirstColIndex
virtual integer iGetFirstColIndex(void) const
Definition: node.cc:88