Term Project Due 11/29/07

CE 5143 Civil Engineering Numerical Analysis, Dr. Chowdhury

Behavior of a Simple Structure with a Nonlinear Element

— by David Wagner 2007/11/17 12:23

Introduction
1. Background
2. Objective
3. Scope
Theory
1. Numerical Method
2. Computer Code
Numerical Analysis
Results and Discussion
Reference
Appendix: Selected Computer Output

Introduction

This paper presents a mathematical model of a simple hydraulic system with nonlinear behavior. A backward-difference approximation gives expected results only over a limited time period.

Some hydraulic systems can be modeled piecewise by mathematically combining the behavior of individual elements. This method is fairly straightforward for linear systems, but becomes more complex for nonlinear systems. Furthermore, a model created in this manner may be limited in accuracy to predictive periods too short to be of use. The simple hydraulic system analysed here demonstrates these limitations.

Objective

The objective is to model the dynamic behavior a simple hydraulic system with a nonlinear component.

Scope

This paper illustrates the analysis with one example.

Consider a simple industrial waste treatment unit, presented here with parameters chosen to illustrate the nature of the hydraulic model developed to describe it. The process unit is a plug-flow reactor designed to treat the average flow of the effluent expected from two similar industrial processes. Each of these processes is cyclic with a period near 6.28 seconds (1.00 rad/s), during which the process alternately draws some, then empties more water than it drew into the reactor. Together the two processes discharge into the reactor the average net flow rate it was designed to treat. The reactor entrance is oriented vertically and continues above the datum like a standpipe to form a chamber intended to buffer the influent flow. Also at the entrance into the process unit is a check valve to divert the wastewater through a pipe and into the sanitary sewer if influent backs up too high into the standpipe. A valve at the end of this pipe is intended to limit the flow into the sewer when this happens.

The check valve exhibits nonlinear behavior of interest. Because the function is well-behaved and applied extensively in engineering to model similar behavior, the Shockley equation is used here to model the behavior of the check valve.

Theory

Coordinate System

x :: Parallel to Fluid Flow
y :: Vertical
z :: Horizontal and Perpendicular to Fluid Flow

A periodic flow may be characterized by the sum of individual sine-wave flows.

¹⁾
$Q prime prime(t) = {{d^2}Q}/{dt^2} = -I omega^2 sin(omega t + phi)$

The pressure drop (h) across each component of a hydraulic network is related to the flow through it (Q).

Component	Pressure Drop h(Q)	dh/dt	Flow Q(h)	dQ/dt
Open Standpipe	$h = {1/{A_{xz}}}int{}{}{Q(t)}dt$	$dh/dt={1/{A_{xz}}}Q(t)$	$Q=A_{xz}dh/dt$	$dQ/dt=A_{xz}{{d^2}h}/{dt^2}$
Valve
Closed Pipe	$h = {L/{g A_{yz}}} {dQ/dt}$	$dh/dt = {L/{g A_{yz}}}{{d^2}Q}/{dt^2}$	$Q={{g A_{yz}}/L} int{}{}{h(t)dt}$	${dQ/dt}={{g A_{yz}}/L} h$
Check Valve	$h = n h_T ln({{Q}/{Q_s}}+1)$	$dh/dt = ({n h_T}/{Q+{Q_s}}) dQ/dt$	$Q = Q_s(e^{{h}/{nh_t}}-1)$	$dQ/dt = {{Q_S}/{nh_T}} e^{{h}/{nh_T}}$ ²⁾

Conservative Quantities

Conservation of Mass:
Conservation of Momentum: $sum{}{}F = rho Q (overline{u_2} - overline{u_1})$
- Fluid Inertia (mv) = ρA_pfL_pfq_pf
Conservation of Energy: $H = z + d cos theta + {{overline{u}^2}/{2g}}$
- $H = h + {{overline{u}^2}/{2g}}$ :

Fluid in a full horizontal pipe

Momentum
Energy ${1/2}m u^2 = {1/2}rho V u^2 = {1/2}rho A L u^2$

Fluid in a vertical standpipe

Momentum
Energy $m g h_{c} + {1/2}m u^2 = rho V g {1/2}h + {1/2}rho V u^2 = rho A h g {1/2}h + {1/2}rho A h u^2 = {1/2}rho g A h^2 + {1/2}rho A h u^2$

Resonant Frequency of a full pipe.

$omega_{pf,0} = sqrt{g/L_{pf}} = sqrt{9.8/9.8}=$ 1.00 rad/s
$omega_{op,0} = sqrt{g/L_{op}} = sqrt{9.8/98}=$ 0.1 rad/s

Mechanical Impedance: Z=F/v

Statics

$F = ma = {mv^2}/Delta$

Numerical Method

Since the objective is to show how this system behaves in general, a finite-backward-difference approximation should be sufficient to visualize the nature of its behavior.

As both outflows discharge to free surfaces at the datum, conservation of energy requires the net headloss through each pathway to be zero.

$H(t) + h_{sp}=0$
$H(t) + h_{pf}=0$
$H(t) + h_{cv} + h_{op} + h_{ov}=0$
→ $h_{sp} = h_{pf} = h_{cv} + h_{op} + h_{ov}$

The sum of their derivatives must also be zero.

$dH/dt + dh_{sp}/dt=0$
$dH/dt + dh_{pf}/dt=0$
$dH/dt + dh_{cv}/dt + dh_{op}/dt + dh_{ov}/dt=0$
→ $dh_{sp}/dt = dh_{pf}/dt = dh_{cv}/dt + dh_{op}/dt + dh_{ov}/dt$

Although H(t) is not known (the pressure is whatever is necessary to maintain the given flow rate Q(t)), the following relations hold regardless of H(t).

$h_{sp} = h_{pf} = h_{cv} + h_{op} + h_{ov}$
$dh_{sp}/dt = dh_{pf}/dt = dh_{cv}/dt + dh_{op}/dt + dh_{ov}/dt$

Since the fluid is considered incompressible, the flow into any enclosed point must equal the flow out of it, thus the net flow at each connecting node is zero.

$Q_{cv} = Q_{op} = Q_{ov}$
$Q(t) = Q_{sp} + Q_{pf} + Q_{cv}$

Differentiate these.

${dQ_{cv}}/dt = {dQ_{op}}/dt = {dQ_{ov}}/dt$
$dQ/dt = {dQ_{sp}}/dt + {dQ_{pf}}/dt + {dQ_{cv}}/dt$

These equations can be solved in terms of head or flow. Solving for head is conventional, so put everything in terms of h.

${{Q_S}/{nh_T}} e^{{h_cv}/{nh_T}} = {{g A_{op}}/L} h_op = {1/R_{ov}}{{dh_{ov}}/{dt}}$
$dQ/dt = A_{sp}{{d^2}{h_sp}}/{dt^2} + {{g A_{pf}}/L} h_pf + {{Q_S}/{nh_T}} e^{{h_cv}/{nh_T}} =$ $A_{sp}{{d^2}{h_sp}}/{dt^2} + {{g A_{pf}}/L} h_sp + {{Q_S}/{nh_T}} e^{{h_cv}/{nh_T}}$

Following the flow, first consider the head in the standpipe using this last equation and substituting backward-difference formulas (Taylor's series expansions) for the unknown terms. Note the rate of change of the total influent flow (dQ/dt) is known or can be approximated.

$dQ/dt approx A_{sp}( {h_{sp,i}-2h_{sp,i-1}+h_{sp,i-2}} /{(Delta t)^2} ) + {{g A_{pf}}/L} h_sp + {{Q_S}/{nh_T}} e^{{h_{cv,i-1}}/{nh_T}}$
${dQ}/dt approx {A_{sp}/{(Delta t)^2}} ( {h_{sp,i}-2h_{sp,i-1}+h_{sp,i-2}} ) + {{g A_{pf}}/L} h_{sp,i} + {{Q_S}/{nh_T}} e^{{ h_{cv,i-1} }/{nh_T}}$
$( {{A_{sp}}/{(Delta t)^2}} + {{g A_{pf}}/L}) h_{sp,i} approx dQ/dt + {A_{sp}/{(Delta t)^2}} ( 2h_{sp,i-1}-h_{sp,i-2} ) - {{Q_S}/{nh_T}} e^{{ h_{cv,i-1} }/{nh_T}}$
$h_{sp,i} approx {dQ/dt + {A_{sp}/{(Delta t)^2}} ( 2h_{sp,i-1}-h_{sp,i-2} ) - {{Q_S}/{nh_T}} e^{{ h_{cv,i-1} }/{nh_T}}} / { {{A_{sp}}/{(Delta t)^2}} + {{g A_{pf}}/L}}$

Now solve for Q_ov by considering the head losses in the overflow line and using a Taylor's series backward-difference approximation.

$dh_{sp}/dt = dh_{cv}/dt + dh_{op}/dt + dh_{ov}/dt$
$dh_{sp}/dt = {{nh_T}/{Q_ov+Q_s}}{dQ_{ov}/dt} + {L/gA_{op}} {{d^2}{Q_ov}}/{dt^2} +R {dQ_{ov}/dt}$
Let $Delta h_sp = h_{sp,i}-h_{sp,i-1}$ and $Delta Q_ov = Q_{ov,i}-Q_{ov,i-1}$
${Delta h_sp}/{Delta t} = {{nh_T}/{Q_{ov,i}+Q_s}} ({Delta Q_ov}/{Delta t}) + {L/gA_{op}} ({Q_{ov,i}-2Q_{ov,i-1}+Q_{ov,i-2}} /{(Delta t)^2}) + R({Delta Q_ov}/{Delta t})$
${Delta h_sp}{Delta t} = {{nh_T Delta t Delta Q_ov }/{Q_{ov,i}+Q_s}} + {L/gA_{op}} ({Q_{ov,i}-2Q_{ov,i-1}+Q_{ov,i-2}} ) + R Delta t Delta Q_ov$

Organize this a bit.

0=
+ ${L/gA_{op}} ({Q_{ov,i}+Q_s})({Q_{ov,i}-2Q_{ov,i-1}+Q_{ov,i-2}} )$
+ $R Delta t Delta Q_ov (Q_{ov,i}+Q_s)$
- ${Delta h_sp}{Delta t}({Q_{ov,i}+Q_s})$

Expand it out.

0= $nh_T Delta t (Q_{ov,i}-Q_{ov,i-1})$
+ ${L/gA_{op}} ({Q_{ov,i}+Q_s})({Q_{ov,i}-2Q_{ov,i-1}+Q_{ov,i-2}} )$
+ $R Delta t (Q_{ov,i}-Q_{ov,i-1}) (Q_{ov,i}+Q_s)$
- ${Delta h_sp}{Delta t}({Q_{ov,i}+Q_s})$

Expand it further.

0= $nh_T Delta t Q_{ov,i}$
- $nh_T Delta t Q_{ov,i-1}$
+ ${L/gA_{op}} {Q_{ov,i}}^2$
- ${{2L}/gA_{op}} Q_{ov,i-1} {Q_{ov,i}}$
+ ${L/gA_{op}} Q_{ov,i-2} Q_{ov,i}$
+ ${L/gA_{op}} Q_s Q_{ov,i}$
- ${{2L}/gA_{op}} Q_{ov,i-1} Q_s$
+ ${L/gA_{op}} Q_s Q_{ov,i-2}$
+ $R Delta t {Q_{ov,i}}^2$
+ $R Delta t Q_s Q_{ov,i}$
- $R Delta t Q_{ov,i-1} Q_{ov,i}$
+ $R Delta t {Q_s}^2$
- $h_{sp,i} {Delta t} Q_{ov,i}$
+ $h_{sp,i-1} {Delta t} Q_{ov,i}$
- $h_{sp,i} {Delta t} Q_s$
+ $h_{sp,i-1}{Delta t} Q_s$

Put this is in terms of Q_ov,i, the variable of interest.

0= $(L/gA_{op} +R Delta t) {Q_{ov,i}}^2$ + $(nh_T Delta t - {{2L}/gA_{op}} Q_{ov,i-1} + {L/gA_{op}} Q_{ov,i-2} + {L/gA_{op}} Q_s + R Delta t Q_s - R Delta t Q_{ov,i-1} - h_{sp,i}{Delta t} + h_{sp,i-1}{Delta t}) Q_{ov,i}$ $+ ( - nh_T Delta t Q_{ov,i-1} - {{2L}/gA_{op}} Q_{ov,i-1} Q_s + {L/gA_{op}} Q_s Q_{ov,i-2} + R Delta t {Q_s}^2 - h_{sp,i} {Delta t} Q_s + h_{sp,i-1} {Delta t} Q_s)$ .

The outflow at each time step can be found by solving for Q_ov,i.

Computer Code

This model is written as a script file for Octave, an open-source mathematical analysis program similar to MatLab.

# CE 5143 Civil Engineering Numerical Analysis, Dr. Chowdhury
# Fall 2007 Term Project
# By David Wagner

clear;
imax    = 1000;
i       = 1:imax     ;# Index counter
Delta_t = 1       ;# (s) :: Time slice
t       = Delta_t*i  ;
g       = 9.8        ;# (m/s²) :: Acceleration due to Gravity
rho     = 1000       ;# (kg/m³) :: Density of wastewater (assumed incompressible)

#Influent Standpipe
A_sp = 10      ;# (m²) :: Average Cross-sectional Area

#Plug-flow Reactor
L_pf = 9.8     ;# (m) :: Hydraulic Path Length
A_pf = 10      ;# (m²) :: Average Cross-sectional Area

#Overflow Check Valve
h_T  = 2.0     ;# (m) :: Minimum Pressure Threshold
n    = 1.0     ;# :: Emission Coefficient or Ideality Factor
Q_s  = 0.01    ;# :: Leakage Flow

#Overflow Pipe
L_op = 9.8*sqrt(10)     ;# (m) :: Hydraulic Path Length
A_op = 1.0     ;# (m²) :: Average Cross-sectional Area

#Overflow Valve
R_ov = 0.01    ;# (m/unit flow) :: Valve Headloss (Resistance)

#Influent Flow
I_1 = 1   ;omega_1 = 1     ;phi_1 = 0   ;
I_2 = 1   ;omega_2 = 1.1   ;phi_2 = 0   ;

dQ_dt(i) = I_1*omega_1*cos(omega_1*t(i) + phi_1)  \
 + I_2*omega_2*cos(omega_2*t(i) + phi_2) ;

h_sp(1:2) = 0  ;# (m) :: Standpipe Fluid Level, initial values
Q_ov(i)   = 0  ;

i = 3:imax;

h_cv(i-1)=0;

for i = 3:imax

h_sp(i) = ( dQ_dt(i)                               \
 + (A_sp/(Delta_t^2)) .* (2*h_sp(i-1)-h_sp(i-2))   \
 - (Q_s/(n*h_T)) * exp(h_cv(i-1)/(n*h_T)) )        \
 / ( (A_sp/(Delta_t^2)) + ((g*A_pf)/L_pf) );

a= (L_op/(g*A_op) + R_ov*Delta_t) ; # {Q_{ov,i}}^2

b=(0
+ n*h_T*Delta_t
- (2*L_op/(g*A_op)) * Q_ov(i-1)
+ (  L_op/(g*A_op)) * Q_ov(i-2)
+ (  L_op/(g*A_op)) * Q_s
+ R_ov * Delta_t * Q_s
- R_ov * Delta_t * Q_ov(i-1)
- h_sp(i)  * Delta_t
+ h_sp(i-1)* Delta_t ); # Q_ov(i);

c= (0
- n * h_T * Delta_t * Q_ov(i-1)
- (2*L_op/(g*A_op)) * Q_ov(i-1) * Q_s
+ (  L_op/(g*A_op)) * Q_s * Q_ov(i-2)
+ R_ov * Delta_t * Q_s^2
- h_sp(i)   * Delta_t * Q_s
+ h_sp(i-1) * Delta_t * Q_s );


Q_ov(i)=max(real(roots([a,b,c])));
h_cv(i)=n * h_T * log( (Q_ov(i) / Q_s) +1 );

endfor

plot(t,Q_ov, t,h_sp)
print(strcat("project_1_",num2str(imax),".png"));

Numerical Analysis

The system is assumed to be flowing full at all times, and both discharge into water bodies with free surfaces at the same elevation (0 m).

g = 9.8 (m/s²) :: Acceleration due to Gravity
ρ = 1000 (kg/m³) :: Density of wastewater (assumed incompressible)
Q (m³/sec) :: Fluid Flow Through the System

Influent Standpipe
- A_sp = 10 (m²) :: Average Cross-sectional Area
- R_sp = 0 :: Negligible Flow Resistance
- h_sp (m) :: Standpipe Fluid Level
Plug-flow Reactor
- L_pf = 9.8 (m) :: Hydraulic Path Length
- A_pf = 10 (m²) :: Average Cross-sectional Area
- R_pf = 0 :: Negligible Flow Resistance
Overflow Check Valve
- h_T = 2.5 (m) :: Minimum Pressure Threshold
- n = 1.0 :: Emission Coefficient or Ideality Factor
- Q_s = 0.01 :: Leakage Flow
- h_cv (m) :: Pressure Head across the Check Valve
Overflow Pipe
- L_op = 9.8√10 (m) :: Hydraulic Path Length
- A_op = 1 (m²) :: Average Cross-sectional Area
- h_op (m) :: Fluid Level
Outflow Valve
- R_ov = 0.01 (m/unit flow) :: Valve Friction

Free water surfaces are at the datum level.

h_sp,in = h_pf,out = h_ov,out = 0 (m)

Continuity leaves only three unknowns.

h_sp,out = h_pf,in = h_cv,in
h_cv,out = h_op,in
h_op,out = h_ov,in

The desired value Q_ov = h_ov,in/R_ov

For this exercise, two flows of nearly the same frequency are considered.

I₁ = I₂ = 1.0 m³/sec
ω₁ = 1.0 rad/sec; ω₂ = 1.1 ω₁ = 1.1 rad/sec
φ₁ = φ₂ = 0.0
${dQ/dt} = I_1 omega_1 cos(omega_1 t + phi_1) + I_2 omega_2 cos(omega_2 t + phi_2) = 1.0 cos(1.0 t) + 1.1 cos(1.1 t)$
${{d^2}Q}/{dt^2} = -I_1 {omega_1}^2 sin(omega_1 t + phi_1) - I_2 {omega_2}^2 sin(omega_2 t + phi_2) = -1.0 sin(1.0 t) -1.21 sin(1.1 t)$

Results and Discussion

The results were unexpected and appear to show the limitations of the model more clearly than the behavior of this system. The continually increasing flow rate past a few thousand iterations of the model should violate the continuity equations, but it is there nonetheless, likely a result of cumulative arithmetic error. The exponential term often exacerbates such errors.

It is also entirely possible there is a mistake in this model, perhaps more than one. A more systematic modeling method would reduce this kind of error. Unfortunately I am not aware of a published system or open-source software package available to guide physical system modeling efforts. Most modeling software is either very specific to an application (SewerCAD, Dr. Frame, HEC-RAS, and such) or overly general (Mathematica, MatLab, and others). These are useful for practicing in an engineering specialty, or for detailed research, respectively. But there appears to be little available for one-off projects or exploring nonstandard modeling techniques.

Reference

* Applied Nomerical Methods for Engineers and Scientists, 2002, Singiresu S. Rao, Prentice Hall

Fundamentals of Engineering Supplied-reference Handbook, Seventh Edition, 2006, National Council of Examiners for Engineering and Surveying
Non-linear DC Analysis, 2007-06-17, Stefan Jahn, Qucs Technical Papers, http://qucs.sourceforge.net/tech/node16.html, accessed 2007-11-29

In addition, Wikipedia, http://wikipedia.org, was of great use for checking basic concepts.

Appendix: Selected Computer Output

The overflow seems to increase without bounds while the average standpipe level gradually drops.

¹⁾

And, similarly for the next equation.

²⁾ Non-linear DC Analysis, Stefan Jahn, 2007