Tension modulated string

Step 1:

Start with the non-linear tension modulated string model. The model is defined by the following equations:

\[ \begin{align} \rho A \frac{\partial^2 y}{\partial t^2} + EI \frac{\partial^4 y}{\partial x^4} - T_{NL}(y) \frac{\partial^2 y}{\partial x^2} + d_{1} \frac{\partial y}{\partial t} + d_{3} \frac{\partial^3 y}{\partial t^3} = f(x,t) \end{align} \]

where \(T_{NL}\) is the non-linear tension function truncated at the third order:

\[ \begin{aligned} T_{N L}(y) & =T_0+T_1(y(x, t)) \\ & =T_0+\frac{E A}{2 l_0} \int_0^{l_0} y^2(x, t) d x . \end{aligned} \]

Step 2:

Perform a transformation of the PDE with respect to the spatial variable \(x\). We can do this using the transformation kernel \(K(x, \beta_{\mu})\).

\[ \begin{aligned} K\left(x, \beta_\mu\right) & =\sin \left(\frac{\mu \pi}{l_0} x\right) \\ \beta_\mu^4 & =E I\left(\frac{\mu \pi}{l_0}\right)^4+T_0\left(\frac{\mu \pi}{l_0}\right)^2 \end{aligned} \]

where \(\beta_\mu\) is the frequency variable. We see that the \(K(x, \beta_{\mu})\) is a function that generates the fourier sine series of the spatial variable \(x\).

The spatial tranformation of the non-linear term \(T_{1}\) is given as:

\[ \mathcal{T}\left\{T_1\right\} = \bar{b}(\mu,y,\bar{y}) \]

We now can write the transformed PDE as:

\[ \boxed{ \begin{aligned} \rho A \ddot{\bar{y}}(\mu, t) + \left(d_3 \eta_\mu^2+d_1\right) \dot{\bar{y}}(\mu, t) + \beta_\mu^4 \bar{y}(\mu, t)-\bar{b}(\mu, y, \bar{y}) = 0 \end{aligned} } \]

Where \(\eta_\mu=-\left(\mu \pi / l_0\right)^2\). The initial conditions are given as:

\[ \begin{aligned} \bar{y}_0(\mu) &= \mathcal{T}{\left\{y_0(x)\right\}} \\ \dot{\bar{y}}_0(\mu) &= \mathcal{T}{\left\{\dot{y}_0(x)\right\}} = 0 \end{aligned} \]

Step 3:

Discretize the non-linear term \(\bar{b}(\mu, y, \bar{y})\)

\[ \begin{aligned} \bar{b}^d(\mu) & =\eta_\mu^2 \frac{E A \bar{y}(\mu)}{2 l_0} \int_0^{l_0}\left[\sum_{\nu=1}^M \frac{1}{N_\nu} \bar{y}^d(\nu) K(\xi, \nu)\right]^{\prime 2} d x \\ & =\eta_\mu^2 \frac{E A \pi^2 \bar{y}(\mu)}{l_0^4} \sum_{\nu=1}^M \nu^2\left(\bar{y}^d(\nu)\right)^2 \end{aligned} \]

Step 4

Write as a system of first order ODEs:

\[ \begin{aligned} \dot{\bar{y}}(\mu, t) &= \bar{v}(\mu, t) \\ \dot{\bar{v}}(\mu, t) &= \frac{-\left(d_3 \eta_\mu^2+d_1\right)}{\rho A} \bar{v}(\mu, t) - \frac{\beta_\mu^4}{\rho A} \bar{y}(\mu, t) + \frac{\bar{b}(\mu, y, \bar{y})}{\rho A} \end{aligned} \]

This is a system of first-order ODEs. In matrix form, we can write this as:

\[ \boxed{ \begin{aligned} \mathbf{\dot{\bar{y}}} &= \mathbf{\bar{v}} \\ \mathbf{\dot{\bar{v}}} &= -\mathbf{M_v} \mathbf{\bar{v}} - \mathbf{M_y} \mathbf{\bar{y}} + \mathbf{\bar{b}} \end{aligned} } \]

The Matrices \(\mathbf{M_v}\), \(\mathbf{M_y}\) are diagonal matrices whose order is equal to the number of modes \(M\).

\[ \begin{aligned} \mathbf{M_v} &= \text{diag}\left(\frac{d_3 \eta_\mu^2+d_1}{\rho A}\right) \\ \mathbf{M_y} &= \text{diag}\left(\frac{\beta_\mu^4}{\rho A}\right) \end{aligned} \]

The vector \(\mathbf{b}\) contains the non-linear terms and is given as:

\[ \mathbf{\Eta} = \text{diag}\left(\eta_\mu^2\right) \\ T_{1} = \frac{E A \pi^2}{l_0^4} \mathbf{\bar{y}}(\mu)^T \cdot \Eta \cdot \mathbf{\bar{y}} \\ \mathbf{\bar{b}} = T_{1} \text{diag}(\frac{\eta_\mu^2}{\rho A}) \mathbf{\bar{y}} \]

from matplotlib import pyplot as plt

The parameters of the model are:

  • \(N\): number of points in the string
  • \(T\): total time of simulation
  • \(f_s\): sampling rate
  • \(l\): length of string
  • \(A\): string cross-sectional area
  • \(E\): Young’s modulus
  • \(I\): string moment of inertia
  • \(\rho\): density
  • \(Ts_0\): initial tension
  • \(d_{1}\): string frequency independent damping
  • \(d_{3}\): string frequency dependent damping
  • \(M\): number of modes to include in the modal expansion

source

Wave1dSolverTensionModulated

 Wave1dSolverTensionModulated (sampling_rate:int=48000,
                               final_time:float=0.5, n_gridpoints:int=101,
                               length:float=1.0, A:float=1.9634e-07,
                               I:float=2.454e-14, rho:float=7800,
                               E:int=190000000000.0, d1:float=0.004,
                               d3:float=6e-05, Ts0:float=150,
                               n_max_modes:int=50,
                               use_nonlinear:bool=True,
                               method:str='DOP853', rtol:float=1e-12,
                               atol:float=1e-14)

*Tension modulated wave equation solver.

This solver is based on the paper “Sound Synthesis With Tension Modulated Nonlinearities” and the accompanying code to the paper “Physical Modeling using Recurrent Neural Networks with Fast Convolutional Layers”.*

Type Default Details
sampling_rate int 48000 1/s Temporal sampling frequency
final_time float 0.5 s Duration of the simulation
n_gridpoints int 101 pts/m Spatial sampling grid
length float 1.0 m Length at rest
A float 1.9634e-07 m**2 Cross section area
I float 2.454e-14 m**4 Moment of intertia
rho float 7800 kg/m**3 Density
E int 190000000000.0 Pa Young’s modulus
d1 float 0.004 kg/(ms) Frequency independent loss
d3 float 6e-05 kg m/s Frequency dependent loss
Ts0 float 150 N Tension
n_max_modes int 50 Number of modal expansion terms
use_nonlinear bool True Use nonlinear wave equation
method str DOP853 Integration method
rtol float 1e-12 Relative tolerance
atol float 1e-14 Absolute tolerance

source

nonlinear_string_func

 nonlinear_string_func (t, state, n_max_modes, use_nonlinear, pre, H, H_1,
                        M_v, M_y)
Type Details
t time
state state vector
n_max_modes number of modes
use_nonlinear use nonlinear wave equation
pre precalculated constant
H H matrix
H_1 H_1 matrix
M_v M_v matrix
M_y M_y matrix
Returns ndarray state at timestep t and position x (u(x, t))

Test the solver

from physmodjax.solver.generator import (
    Gaussian,
    generate_initial_condition,
)
length = 0.65
A = 0.5188e-6
I = 0.141e-12
rho = 1140
E = 5.4e9 
d1 = 8e-5
d3 = 1.4e-5
Ts0 = 60.97
n_gridpoints = 101
n_max_modes = 99
sr = 96000


solver = Wave1dSolverTensionModulated(
    sampling_rate=sr,
    final_time=0.1,
    A=A,
    I=I,
    rho=rho,
    E=E,
    d1=d1,
    d3=d3,
    Ts0=Ts0,
    length=length,
    n_gridpoints=n_gridpoints,
    n_max_modes=n_max_modes,
    use_nonlinear=True,
    # d3=1e-7,
)

# initial conditions
ic_max_amplitude = 0.02
u0, v0 = generate_initial_condition(
    rng=np.random.default_rng(354),
    generator=Gaussian(num_points=n_gridpoints),
    ic_max_amplitude=ic_max_amplitude,
    ic_amplitude_random=False,
)

t, u, v = solver.solve(u0, v0)
print(f"u shape: {u.shape}")
dx: 0.006500000000000001 in meters
dt: 1.0416666666666666e-05 in seconds
number of points (n_gridpoints): (101,)
time in samples (nt): (9600,)
n_samples_plot = 7000
fig, ax = plt.subplots(1, 2, figsize=(5, 6))
ax[0].pcolormesh(solver.grid, t[-n_samples_plot:], u[-n_samples_plot:, :])
ax[0].set_title("u")
ax[0].set_xlabel("x")
ax[0].set_ylabel("t")
ax[1].pcolormesh(solver.grid, t[-n_samples_plot:], v[-n_samples_plot:, :])
ax[1].set_title("v")
ax[1].set_xlabel("x")
ax[1].set_ylabel("t")
fig.tight_layout()
plt.show()
# plot the solution at the endpoints
plt.figure(figsize=(10, 5))
plt.plot(t, u[:, 0], label="left")
plt.plot(t, u[:, -1], label="right", alpha=0.5)
plt.legend()
plt.show()

# Plot FFT of the solution for a point on the string, in db scale
plt.figure(figsize=(10, 5))
freq = np.fft.rfftfreq(u.shape[-1], d=solver.dt)
plt.plot(freq, 20 * np.log10(np.abs(np.fft.rfft(u[17, :]))), label="fft")
# plt.figure(figsize=(10, 5))
# freq = np.fft.rfftfreq(u.shape[0], d=solver.dt)
# plt.plot(freq, np.abs(np.fft.rfft(u[:, 17])), label="fft")
# plt.xlim(3000, 8000)
plt.xlim(100, 5000)
plt.legend()
plt.show()

# Plot spectrogram of the solution for a point on the string
plt.figure(figsize=(10, 5))
plt.specgram(u[:, 17], Fs=solver.sampling_rate, NFFT=2048, noverlap=1024)
plt.ylim(0, 20000)
plt.show()

# Plot FFT of the initial condition
plt.figure(figsize=(10, 5))
freq = np.fft.rfftfreq(u.shape[-1], d=solver.dx)
# plt.plot(freq, np.abs(np.fft.rfft(u0)), label="fft")
plt.plot(np.abs(np.fft.rfft(u0)), label="fft")
# plt.xlim(0, 1)
plt.legend()

# plot the initial conditions
plt.figure(figsize=(10, 5))
plt.plot(solver.grid, u0, label="u0")
plt.plot(solver.grid, u[0, :], label="u0_sol", alpha=0.5)
plt.legend()
plt.show()

# Plot the difference between the initial condition given and the one obtained from the solution
plt.figure(figsize=(10, 5))
plt.plot(solver.grid, u0 - u[0], label="u0 - u0_sol")
plt.legend()
plt.show()

from IPython import display as ipd
# Play the solution as audio
ipd.Audio(u[:, 27], rate=solver.sampling_rate)