Linearization of a Nonlinear Differentials Equations

Linearization is the process of taking the gradient of a nonlinear function with respect to all variables and creating a linear representation at that point. It is required for certain types of analysis such as stability analysis, solution with a Laplace transform, and to put the model into linear state-space form. Consider a nonlinear differential equation model that is derived from balance equations with input u and output y.

$$\frac{dy}{dt} = f(y,u)$$

The right hand side of the equation is linearized by a Taylor series expansion, using only the first two terms.

$$\frac{dy}{dt} = f(y,u) \approx f \left(\bar y, \bar u\right) + \frac{\partial f}{\partial y}\bigg|_{\bar y,\bar u} \left(y-\bar y\right) + \frac{\partial f}{\partial u}\bigg|_{\bar y,\bar u} \left(u-\bar u\right)$$

If the values of `\bar u` and `\bar y` are chosen at steady state conditions then `f(\bar y, \bar u)=0` because the derivative term `{dy}/{du}=0` at steady state. To simplify the final linearized expression, deviation variables are defined as `y' = y-\bar y` and `u' = u - \bar u`. A deviation variable is a change from the nominal steady state conditions. The derivatives of the deviation variable is defined as `{dy'}/{dt} = {dy}/{dt}` because `{d\bar y}/{dt} = 0` in `{dy'}/{dt} = {d(y-\bar y)}/{dt} = {dy}/{dt} - \cancel{{d\bar y}/{dt}}`. If there are additional variables such as a disturbance variable `d` then it is added as another term in deviation variable form `d' = d - \bar d`.

$$\frac{dy'}{dt} = \alpha y' + \beta u' + \gamma d'$$

The values of the constants `\alpha`, `\beta`, and `\gamma` are the partial derivatives of `f(y,u,d)` evaluated at steady state conditions.

$$\alpha = \frac{\partial f}{\partial y}\bigg|_{\bar y,\bar u,\bar d} \quad \quad \beta = \frac{\partial f}{\partial u}\bigg|_{\bar y,\bar u,\bar d} \quad \quad \gamma = \frac{\partial f}{\partial d}\bigg|_{\bar y,\bar u,\bar d}$$

Exercise 1

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), Cd as the drag coefficient (0.24), `\rho` as the air density (1.225 kg/m3), A as the vehicle cross-sectional area (5 m2), and Fp 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.

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.

# analytic solution
import sympy as sp
# 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)

# numeric solution
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('Exact Partial Derivatives')
print(-A*Cd*rho*v/m) # exact d(f(u,v))/dv
print(Fp/m) # exact d(f(u,v))/du

Exercise 2

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, cp as the heat capacity, `\rho` as the density, and dt as the diameter of the thermocouple bead.

$$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 below. 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?

import numpy as np
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.plot(t,Tg,'b--',linewidth=3,label='Gas Temperature')
plt.plot(t,Tf,'g:',linewidth=3,label='Flame Temperature')
plt.ylabel('Temperature (K)')
plt.xlabel('Time (sec)')