Message Boards Message Boards

NDSolve piecewise differential equations associated with a water hammer?

Posted 5 years ago

I am trying to solve the differential equations associated with water hammer. I have three pipe segments and I am using piecewise function for the differential equations. When I try to solve using NDSolve I get an error message " The function value {0,.....} is not a list of numbers with dimension {50} at {x,P[x,t],V[x,t]}={...}". Can someone please help me with this.

L1 = 3100; (*Length of first pipe*)
L2 = 2700;(*Length of second pipe*)

L3 = UnitConvert[Quantity[30., "ft"], "m"][[
   1]];(*Length of third pipe*)

D1 = UnitConvert[Quantity[4, "in"], "m"][[
   1]]; (*Diameter of first pipe *)

D2 = UnitConvert[Quantity[6, "in"], "m"][[
   1]];(*Diameter of second pipe *)

D3 = UnitConvert[Quantity[1.75, "in"], "m"][[
   1]];(*Diameter of third pipe *)

f = 0.002; (*Friction factor*)
\[Rho] = 1000; (*Density of fluid*)

Ef = 2.19*^9; (*Bulk modulous of fluid*)

Ep = 210*^9;(*Elastic modulous of pipe*)

ee = 4.5*^-5; (* Pipe roughness*)
\[Mu] = 
  0.001002; (* Fluid viscosity in Poise*)
Qmax = 0.1;
A1 = Pi/4 D1^2;
A2 = Pi/4 D2^2;
A3 = Pi/4 D3^2;
w = UnitConvert[Quantity[0.2, "in"], "m"][[1]];
solEe = Solve[1/EE == 1/Ef + 1/(w Ep), EE][[1]];
Ee = EE /. solEe;
c = Sqrt[Ee/\[Rho]];
pde1 = D[P[x, t], t] + \[Rho] c^2 D[V[x, t], x] == 0
sol1 = NDSolve[{pde1,
   D[V[x, t], t] + 1/\[Rho] D[P[x, t], x] == 
    Piecewise[{{-((f V[x, t] Abs[V[x, t]])/(2 D1)), 0 < x < L1},
      {-((f V[x, t] Abs[V[x, t]])/(2 D2)) == 0, L1 < x < L1 + L2},
      {-((f V[x, t] Abs[V[x, t]])/(2 D3)) == 0, 
       L1 + L2 < x < L1 + L2 + L3}}],

   V[x, 0] == 0.,
   P[x, 0] == 0.,

   P[L1 + L2 + L3, t] == \[Rho]/(2 A3^2) (A3 V[L1 + L2 + L3, t])^2,
   V[L1 + L2 + L3, t] == 10 t}, {P, V}, {x, 0, L1 + L2 + L3}, {t, 0, 
   12}]
POSTED BY: srinivas.gk
2 Replies
Posted 5 years ago

I now get the following error. How do I resolve this. I tried using Method->"StiffnessSwitching" but it doesnt help. NDSolve::ndsz: At x == 662.1269521572514`, step size is effectively zero; singularity or stiff system suspected.

L1 = 3100; (*Length of first pipe*)
L2 = 2700;(*Length of second pipe*)
L3 = UnitConvert[Quantity[30., "ft"], "m"][[
   1]];(*Length of third pipe*)
D1 = UnitConvert[Quantity[4, "in"], "m"][[
   1]]; (*Diameter of first pipe *)
D2 = UnitConvert[Quantity[6, "in"], "m"][[
   1]];(*Diameter of second pipe *)
D3 = UnitConvert[Quantity[1.75, "in"], "m"][[
   1]];(*Diameter of third pipe *)
f = 0.002; (*Friction factor*)
\[Rho] = 1000; (*Density of fluid*)
Ef = 2.19*^9; (*Bulk modulous of fluid*)
Ep = 210*^9;(*Elastic modulous of pipe*)
ee = 4.5*^-5; (* Pipe roughness*)
\[Mu] = 
  0.001002; (* Fluid viscosity in Poise*)
Qmax = 0.1;
A1 = Pi/4 D1^2;
A2 = Pi/4 D2^2;
A3 = Pi/4 D3^2;
w = UnitConvert[Quantity[0.2, "in"], "m"][[1]];
solEe = Solve[1/EE == 1/Ef + 1/(w Ep), EE][[1]];
Ee = EE /. solEe;
c = Sqrt[Ee/\[Rho]];
pde1 = D[P[x, t], t] + \[Rho] c^2 D[V[x, t], x] == 0
sol1 = NDSolve[{pde1,
   D[V[x, t], t] + 1/\[Rho] D[P[x, t], x] == 
    Piecewise[{{-((f V[x, t] Abs[V[x, t]])/(2 D1)), 0 < x < L1},
      {-((f V[x, t] Abs[V[x, t]])/(2 D2)), L1 < x < L1 + L2},
      {-((f V[x, t] Abs[V[x, t]])/(2 D3)), 
       L1 + L2 < x < L1 + L2 + L3}}],

   V[x, 0] == 0.,
   P[x, 0] == 0.,

   P[L1 + L2 + L3, t] == \[Rho]/(2 A3^2) (A3 V[L1 + L2 + L3, t])^2,
   V[L1 + L2 + L3, t] == 10 t}, {P, V}, {x, 0, L1 + L2 + L3}, {t, 0, 
   12}, Method -> "StiffnessSwitching", AccuracyGoal -> 5, 
  PrecisionGoal -> 4]
POSTED BY: srinivas.gk
Posted 5 years ago

I found the error. I have == twice in the same equation.

POSTED BY: srinivas.gk
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract