HW09: Additional Equation Solving Practice#

Problem 1#

This problem reviews the use of non-linear solvers. You can use any numerical non-linear solver. Below is an image of a portion of a horizonal pipe network with pipes label P1, P2, and P3. Fluid enters the left and then exits at the right through both P2 and P3. Each of those pipes have different properties (lengths, diameters, and roughness). Your task is to solve for the flow in P2 and P3 given a flow in P1.

Action 1: Create a class called fluid using param.Parameterized that has the following properties:

  • density (default of 1000 kg/m3)

  • mu (default of 1E-3 Pa-sec for viscosity)

#Your code here

Action 2: Create a class called pipeClass using param.Parameterized that has the following properties:

  • length (m)

  • diameter (m)

  • roughness (m)

  • flow (m3/s)

  • fluidinside (a fluid object with the default being the fluid default): specify with fluidinside = param.Parameter(fluid()). You can refer to the fluidinside density or mu with self.fluidinside.density or self.fluidinside.mu, respectively.

  • a function called velocity that returns the velocity of the fluid given self.flow and self.diameter: i.e. velocity = self.flow/(np.pi/4*self.diameter**2)

  • a function called reynoldsNumber that returns the Reynolds number for the pipe flow: i.e. Re = self.velocity()*self.diameter*self.fluidinside.density/self.fluidinside.mu

    • for the Reynolds number, assume a density of 1000 kg/m3 and a viscosity of 1E-3 Pa*sec (the default fluid() parameters)

  • a function called frictionFactor that returns the friction factor for the pipe using the the Haaland Equation

  • a function to determine the pressuredrop in the horizontal pipe = \(\rho fLv^2/(2D)\). \(\rho\) is density, f the friction factor, L the length, v the velocity, D the diameter.

All of the functions in the class need only be functions of self

#Your code here

Action 3: Now create three pipe objects with the following properties:

Pipe

Length (m)

Diameter (m)

Roughness (m)

P1

100

0.5

0.0001

P2

200

0.25

0.00001

P3

50

0.75

0.0002

Other System Properties:

  • The total flow rate through the system is 0.1 m^3/s

  • The fluid is water with the properties above referenced.

  • The fluid is turbulent (Re > 4000)

#Your code here

Action 4: Setup a function to calculate the pressure drop through P1 and a second pipe object with a given flow rate in the second pipe object using the above objects and their class functions, e.g. def routepressuredrop(objectj,flowj):. You can do this step with three lines of object based code.

#Your code here

Action 5: Then setup a function to solve for the flow in each pipe. The flow in each pipe will be different such that the total pressure drop through P1 and P3 is equal to P1 and P2, e.g. routepressuredrop(P2,flow2) = routepressuredrop(P3,flow3)

#Your code here

Action 6: Report the flow rate in each pipe and verify that the flow is indeed turbulent(Re > 4000). Also report the total pressure drop through the system (P1 + P3) and verify that it is equal to the pressure drop through P1 and P2.

#Results here

Hopefully, this example shows you how using objects helps you keep track of the variables in your code. And that it helps reduce the errors made.

Problem 2#

Program the same problem above without using objects. I started to do this for the key and I thought it would be rude of me to ask you to actually do this so this problem is a free-bee. I’m giving you the result for free below. It doesn’t look as nice and organized as the code is above (especially if I don’t name my variables nicely).

from scipy.optimize import fsolve
import numpy as np
rho = 1000
mu = 1e-3
pipe1_length = 100
pipe2_length = 200
pipe3_length = 50
pipe1_dia = 0.5
pipe2_dia = 0.25
pipe3_dia = 0.75
pipe1_r = 0.0001
pipe2_r = 0.00001
pipe3_r = 0.0002
pipe1_flow = 0.1

def velocity(flow,D):
    return flow/(np.pi/4*D**2)

def Reynolds(D,v):
    return rho * D * v/mu

def frictionFactor(D,v,rough): #from Haaland equation
    return (1/(-1.8*np.log10(6.9/Reynolds(D,v)+(rough/(3.7*D))**(10/9))))**2

def pressuredrop(f,L,D,v):
    return rho*f*L/D/2*v**2

def legpressuredrop12(flow2):
    p1drop = pressuredrop(frictionFactor(pipe1_dia,velocity(pipe1_flow,pipe1_dia),pipe1_r),pipe1_length,pipe1_dia,velocity(pipe1_flow,pipe1_dia))
    p2drop = pressuredrop(frictionFactor(pipe2_dia,velocity(flow2,pipe2_dia),pipe2_r),pipe2_length,pipe2_dia,velocity(flow2,pipe2_dia))
    return p1drop + p2drop

def legpressuredrop13(flow3):
    p1drop = pressuredrop(frictionFactor(pipe1_dia,velocity(pipe1_flow,pipe1_dia),pipe1_r),pipe1_length,pipe1_dia,velocity(pipe1_flow,pipe1_dia))
    p3drop = pressuredrop(frictionFactor(pipe3_dia,velocity(flow3,pipe3_dia),pipe3_r),pipe3_length,pipe3_dia,velocity(flow3,pipe3_dia))
    return p1drop + p3drop

def eqh(flow2):
    flow3 = pipe1_flow - flow2
    return legpressuredrop12(flow2) - legpressuredrop13(flow3)

pipe2_flow = fsolve(eqh,0.03)[0]

pipe2_flow

0.0024372371378035057

By my count 32 lines of code without objects. 28 lines of code with objects. Which way do you like better? Why?

Your response here

Problem 3#

In class we used lambda for an inline function: f = lambda x: x**2 -5. I can use a lambda function to define the arguments to use in a multivariable function as follows:

Example: Solve this equation for y given: \(x^2+y^2+c =56.3\) where \(x = 3\) and \(c = -4.2\)

#solve this equation:
def func(y,x,c):
    return x**2 + y**2 + c - 56.3

Now solve this equation with x = 3 and c = -4.2 with a guess of 3. Remembering that Call like this: fsolve(function,guess).

x = 3; c = -4.2
fsolve(lambda var: func(var,x,c),3)[0]

7.176350047203662

Notice that I used the inline function definition with lambda to reduce my variable count to 1. Another way is with arguments:

x = 3; c = -4.2
fsolve(func,3,args=(x,c,))[0] #notice the variable I'm solving for is listed first in the function func

7.176350047203662

Action Solve the below set of equations given b = 3 and c = 5. You must keep b and c as parameters in the function call.

\[\begin{split} 2\cdot x^2 + y^2 = 25 \\ x + y + z = c \\ y + 4.3\cdot z = b \end{split}\]
#Your code here#