## Linearization of Differential Equations

Linearize the nonlinear differential equations and compare the linearized model simulation to the original nonlinear model response.

#### Exercise 1 - Automobile

Linearize the momentum balance for the velocity of an automobile at steady state conditions when the gas pedal is maintained at 40%.

$$m\frac{dv(t)}{dt} = F_p u(t) - \frac{1}{2} \rho \, A \, C_d v(t)^2 $$

with *u(t)* as gas pedal position (%pedal), *v(t)* as velocity (m/s), *m* as the mass of the vehicle (500 kg) plus the mass of the passengers and cargo (200 kg), *C _{d}* as the drag coefficient (0.24), `\rho` as the air density (1.225 kg/m

^{3}),

*A*as the vehicle cross-sectional area (5 m

^{2}), and

*F*as the thrust parameter (30 N/%pedal). Simulate a step change in the gas pedal position, comparing the nonlinear response to the linear response. Comment on the accuracy of the linearized model.

_{p}

**Solution Help**

Python can calculate either analytic (exact) or numeric (approximate) derivatives. This problem is not very complicated but see this video tutorial for SymPy and NumPy tutorials on derivatives.

import sympy as sp

sp.init_printing()

# define symbols

v,u = sp.symbols(['v','u'])

Fp,rho,Cd,A,m = sp.symbols(['Fp','rho','Cd','A','m'])

# define equation

eqn = Fp*u/m - rho*A*Cd*v**2/(2*m)

print(sp.diff(eqn,v))

print(sp.diff(eqn,u))

from scipy.misc import derivative

m = 700.0 # kg

Cd = 0.24

rho = 1.225 # kg/m^3

A = 5.0 # m^2

Fp = 30.0 # N/%pedal

u = 40.0 # % pedal

v = 50.0 # km/hr (change this for SS condition)

def fv(v):

return Fp*u/m - rho*A*Cd*v**2/(2*m)

def fu(u):

return Fp*u/m - rho*A*Cd*v**2/(2*m)

print('Approximate Partial Derivatives')

print(derivative(fv,v,dx=1e-4))

print(derivative(fu,u,dx=1e-4))

print('Exact Partial Derivatives')

print(-A*Cd*rho*v/m) # exact d(f(u,v))/dv

print(Fp/m) # exact d(f(u,v))/du

The differential equation in linear form becomes:

$$\frac{dv(t)}{dt} = \alpha \left(u(t) - 40\right) + \beta \left(v(t)-v_{ss}\right) $$

with the following definitions from the partial derivatives

$$\alpha = \frac{F_p}{m} = 0.042857$$

$$\beta = -\frac{A \, C_d \, \rho \, v_ss}{m} = constant$$

Complete the remainder of this exercise by calculating the steady state velocity `v_{ss}`, the constant `\beta`, and simulating a step response with the linearized and nonlinear models.

**Steady State Velocity**

The steady state velocity is determined by setting the time derivative to zero, plugging in other steady state conditions, and solving for `v_{ss}`.

$$m\frac{dv(t)}{dt} = F_p u(t) - \frac{1}{2} \rho \, A \, C_d v(t)^2$$

$$0 = F_p u_{ss} - \frac{1}{2} \rho \, A \, C_d v_{ss}^2$$

$$v_{ss} = \sqrt{\frac{2 \, F_p u_{ss}}{\rho \, A \, C_d}}$$

$$v_{ss} = \sqrt{\frac{2 \times 30 \frac{N}{\%} \times 40 \%}{1.225 \frac{kg}{m^3} \times 5 m^2 \times 0.24}} = 40.406 \frac{m}{s}$$

**Nonlinear Model Simulation**

import matplotlib.pyplot as plt

from scipy.integrate import odeint

# define model

def vehicle(v,t,u,load):

# inputs

# v = vehicle velocity (m/s)

# t = time (sec)

# u = gas pedal position (0% to 100%)

# load = passenger load + cargo (kg)

Cd = 0.24 # drag coefficient

rho = 1.225 # air density (kg/m^3)

A = 5.0 # cross-sectional area (m^2)

Fp = 30 # thrust parameter (N/%pedal)

m = 500 # vehicle mass (kg)

# calculate derivative of the velocity

dv_dt = (1.0/(m+load)) * (Fp*u - 0.5*rho*Cd*A*v**2)

return dv_dt

tf = 60.0 # final time for simulation

nsteps = 61 # number of time steps

delta_t = tf/(nsteps-1) # how long is each time step?

ts = np.linspace(0,tf,nsteps) # linearly spaced time vector

# simulate step test operation

step = np.zeros(nsteps) # u = valve % open

step[11:] = 50.0 # step up pedal position

# passenger(s) + cargo load

load = 200.0 # kg

# velocity initial condition

v0 = 0.0

# for storing the results

vs = np.zeros(nsteps)

# simulate with ODEINT

for i in range(nsteps-1):

u = step[i]

v = odeint(vehicle,v0,[0,delta_t],args=(u,load))

v0 = v[-1] # take the last value as initial condition

vs[i+1] = v0 # store the velocity for plotting

# plot results

plt.figure()

plt.subplot(2,1,1)

plt.plot(ts,vs,'b-',linewidth=3)

plt.ylabel('Velocity (m/s)')

plt.legend(['Velocity'],loc=2)

plt.subplot(2,1,2)

plt.plot(ts,step,'r--',linewidth=3)

plt.ylabel('Gas Pedal')

plt.legend(['Gas Pedal (%)'])

plt.xlabel('Time (sec)')

plt.show()

The next step is to include the linearized model and compare to the nonlinear model response.

#### Exercise 2 - Thermocouple

Linearize an energy balance that describes the temperature of a thermocouple bead `T_t` given a change in the surrounding gas temperature `T_g`. Include the flame temperature `T_f` as a variable in the linearization. Simulate a cyclical gas temperature change and a step change in flame temperature in both the linear and nonlinear versions of the model.

Some advanced combustion systems operate in a pulsed mode where an acoustic wave resonates in the combustor in order to enhance heat and mass transfer. The frequency of the wave is approximately 100Hz (cycles/sec). As a result, the temperature at any given spatial location in the combustor fluctuates with time. The fluctuations can be approximated as a sine wave with an amplitude of `75^oC` for the system of interest. The objective is to measure the fluctuating temperature at a given point in the reactor with use of a thermocouple.

Heat is transferred to the thermocouple by convection and radiation. Heat transfer by radiation is approximated by `q = A \epsilon \sigma (T_f^4-T_t^4)` where `A` is the surface area, `\epsilon` is the emissivity, and `\sigma` is the Stefan-Boltzmann constant `(W/(m^2 K^4))`. The thermocouple bead is approximated as a sphere (e.g. ignore conduction along connecting wires).

Use the following values with subscript *t* as the thermocouple, subscript *f* as the flame, *ss* indicating steady-state, *h* as the heat transfer coefficient, *c _{p}* as the heat capacity, `\rho` as the density, and

*d*as the diameter of the thermocouple bead.

_{t}$$h=2800 \, \frac{W}{m^2 \, K} \quad \rho_{t} = 20 \, g/cm^3 \quad \sigma = 5.67e-8 \, \frac{W}{m^2\,K^4}$$ $$\epsilon=0.8 \quad T_{t,ss} = 1500 \, K \quad T_{f} = 1500 \, K$$ $$c_{p,t} = 0.4 \frac{J}{g\,K} \quad d_t = 0.01 \, cm$$

A starting Python script is available in the solution help. Add the linearized balance equation in addition to the nonlinear balance equation that is already implemented. Simulate both for the cycling gas temperature and discuss whether the linearized model is a suitable approximation of the nonlinear model. What changes could be made to the thermocouple to improve the ability to measure gas temperature?

**Solution Help**

import matplotlib.pyplot as plt

from scipy.integrate import odeint

# define thermocouple model

def thermocouple(x,t,Tg,Tf):

## States

Tt = x[0]

Tlin = x[1]

## Parameters

h = 2800.0 # W/m^2-K

rho = 20.0 # gm/cm^3

sigma = 5.67e-8 # W/m^2-K^4

eps = 0.8 #

Cp = 0.4 # J/gm-K

d = 0.01 # cm

r = d / 2.0 # radius

A = 4.0 * np.pi * (r/100.0)**2 # sphere area (m^2)

V = 4.0/3.0 * np.pi * r**3 # sphere volume (cm^3)

# acc = inlet - outlet

# acc = m * Cp * dT/dt = rho * V * Cp * dT/dt

# inlet - outlet from 2 heat transfer pathways

# q(convection) = h * A * (Tg-Tt)

# q(radiation) = A * esp * sigma * (Tf^4 - Tt^4)

q_conv = h * A * (Tg-Tt)

q_rad = A * eps * sigma * (Tf**4 - Tt**4)

dTdt = (q_conv + q_rad) / (rho * V * Cp)

dTlin_dt = 0.0 # add linearized equation

return [dTdt,dTlin_dt]

# Flame temperature

Tf0 = 1500.0 #K

# Starting thermocouple temperature

y0 = [Tf0,Tf0]

# Time Intervals (sec)

t = np.linspace(0,0.1,1000)

# Flame Temperature

Tf = np.ones(len(t))*Tf0

Tf[500:] = 1520.0 #K

# Gas temperature cycles

Tg = Tf + (150.0/2.0) * np.sin(2.0 * np.pi * 100.0 * t) #K

# Store thermocouple temperature for plotting

Tt = np.ones(len(t)) * Tf

Tlin = np.ones(len(t)) * Tf

for i in range(len(t)-1):

ts = [t[i],t[i+1]]

y = odeint(thermocouple,y0,ts,args=(Tg[i],Tf[i]))

y0 = y[-1]

Tt[i+1] = y0[0]

Tlin[i+1] = y0[1]

# Plot the results

plt.figure()

plt.plot(t,Tg,'b--',linewidth=3,label='Gas Temperature')

plt.plot(t,Tf,'g:',linewidth=3,label='Flame Temperature')

plt.plot(t,Tt,'r-',linewidth=3,label='Thermocouple')

plt.ylabel('Temperature (K)')

plt.legend(loc='best')

plt.xlim([0,t[-1]])

plt.xlabel('Time (sec)')

plt.show()