from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimationTension modulated stiff membrane
The tension is given by \[ T(u) = T_0 + T_{NL}(u) \]
We consider, \[ T_{N L}(u)=C_{N L} \frac{S(u)-S_0}{S_0} \simeq \frac{1}{2} \frac{C_{N L}}{S_0} \int_{\mathcal{S}}\|\nabla u\|^2 d \mathbf{x} \] For the rectangular plate, the Berger approximation is used for the tension \[ \begin{aligned} T_{N L}(u) \simeq & \frac{Q h}{2 L_x L_y\left(1-\nu^2\right)} \\ & \cdot \int_0^{L_x} \int_0^{L_y}\left[\left(\frac{\partial u}{\partial x}\right)^2+\left(\frac{\partial u}{\partial y}\right)^2\right] d x d y . \end{aligned} \] And therefore \[ C_{N L}=\frac{Q h}{\left(1-\nu^2\right)} \]
\[ \boxed{ \begin{aligned} \rho h \ddot{\bar{u}}_{n, m}(t) + \left(d_3 \lambda_{n, m}+d_1\right) \dot{\bar{u}}_{n, m}( t) + \left(\lambda_{n, m} \left(\lambda_{n, m}D +T_0\right)\right) \bar{u}_{n, m}(t)-\bar{f}^{(tm)}_{n, m}(u, \bar{u}) = 0 \end{aligned} } \]
Where \[ \begin{aligned} f^{(tm)}_{n, m}(u, \bar{u}) & = -\lambda_{n, m} T_{N L} (u) \bar{u}_{n, m}(t) \\ & = -\lambda_{n, m} \frac{1}{2}\frac{C_{N L}}{S_0} \left[\sum_{\tilde{n}, \tilde{m}} \frac{\lambda_{\tilde{n}, \tilde{m}}\bar{u}_{\tilde{n}, \tilde{m}}^2 (t)}{\lVert K_{\tilde{n}, \tilde{m}} \rVert_{2}^{2}} \right] \bar{u}_{n, m}(t) \end{aligned} \]
For a rectangular plate, the norm of the modes is given by \[\lVert K_{n, m} \rVert_{2}^{2} = \langle K_{n, m} , K_{n, m} \rangle = \int _0^{L_x}\int _0^{L_y}\sin ^2\left(\frac{n \pi x}{L_x}\right) \sin ^2\left(\frac{m \pi y}{L_y}\right)dydx = \frac{L_x L_y}{4} \quad \forall n, m \in \mathbb{Z} \]
Write as a system of first order ODEs:
\[ \begin{aligned} \dot{\bar{u}}(\mu, t) &= \bar{v}(\mu, t) \\ \dot{\bar{v}}(\mu, t) &= \frac{-\left(d_3 \lambda_\mu+d_1\right)}{\rho h} \bar{v}(\mu, t) - \frac{\beta_\mu}{\rho h} \bar{u}(\mu, t) + \frac{f^{(tm)}_{\mu}(u, \bar{u})}{\rho h} \end{aligned} \]
Where \[ \begin{aligned} \beta_\mu &= \lambda_\mu \left(\lambda_\mu D + T_0\right) = \lambda_\mu^2 D + \lambda_\mu T_0 \\ \bar{b}(\mu, u, \bar{u}) &= \frac{f^{(tm)}_{\mu}(u, \bar{u})}{\rho h} \\ &= -\frac{1}{\rho h} \lambda_{\mu} \frac{1}{2}\frac{Q h}{\left(1-\nu^2\right)}\frac{1}{L_x L_y} \frac{4}{L_x L_y}\left[\sum_{\eta} \lambda_{\eta}\bar{u}_{\eta}^2 (t) \right] \bar{u}_{\mu}(t) \\ &= -\lambda_{\mu} \frac{Q}{\rho \left(1-\nu^2\right)}\frac{2}{L_x^2 L_y^2}\left[\sum_{\eta} \lambda_{\eta}\bar{u}_{\eta}^2 (t) \right] \bar{u}_{\mu}(t) \\ \end{aligned} \]
${u}{n, m}(t) = {u}{}(t) $ is a matrix, but to be able to feed it to solve_ivp, we will flatten it to a vector, and index it as $= n N_y + m $ where \(N_y\) is the number of modes in the y direction.
This is a system of first-order ODEs. In matrix form, we can write this as:
\[ \boxed{ \begin{aligned} \mathbf{\dot{\bar{u}}} &= \mathbf{\bar{v}} \\ \mathbf{\dot{\bar{v}}} &= -\mathbf{M_v} \mathbf{\bar{v}} - \mathbf{M_u} \mathbf{\bar{u}} + \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{\Lambda} &= \text{diag}\left(\lambda_\mu\right) \\ \mathbf{M_v} &= \text{diag}\left(\frac{d_3 \lambda_\mu + d_1}{\rho h}\right) \\ \mathbf{M_y} &= \text{diag}\left(\frac{\beta_\mu}{\rho h}\right) \\ \mathbf{\bar{b}}(\mathbf{u}, \mathbf{\bar{u}}) &= - \frac{Q}{\rho \left(1-\nu^2\right)}\frac{2}{L_x^2 L_y^2}\mathbf{\Lambda} \left[ \mathbf{\bar{u}}^{T} \mathbf{\Lambda} \mathbf{\bar{u}} \right] \mathbf{\bar{u}} \\ \mathbf{\bar{b}}(\mathbf{u}, \mathbf{\bar{u}}) &= - \mathbf{C_b}\mathbf{\Lambda} \left[ \mathbf{\bar{u}}^{T} \mathbf{\Lambda} \mathbf{\bar{u}} \right] \mathbf{\bar{u}} \\ \end{aligned} \]
Wave2dSolverTensionModulated
Wave2dSolverTensionModulated (sampling_rate:int=16000, final_time:float=0.5, n_gridpoints_x:int=41, length_x:float=0.4, aspect_ratio:float=0.8, rho:float=1380, h:float=0.00019, E:int=3500000000.0, nu:float=0.3, d1:float=8e-05, d3:float=1.4e-05, Ts0:float=2620, n_max_modes:int=36, use_nonlinear:bool=True)
Tension modulated wave equation solver for a rectangular stiff membrane. The parameters were taken from (Fletcher, 1991, p.86) and adapted to a rectangular case.
| Type | Default | Details | |
|---|---|---|---|
| sampling_rate | int | 16000 | 1/s Temporal sampling frequency |
| final_time | float | 0.5 | s Duration of the simulation |
| n_gridpoints_x | int | 41 | pts/m Spatial sampling grid |
| length_x | float | 0.4 | m Length of x dimension |
| aspect_ratio | float | 0.8 | Aspect ratio of the membrane, Ly/Lx |
| rho | float | 1380 | kg/m**3 Density |
| h | float | 0.00019 | m Thickness |
| E | int | 3500000000.0 | Pa Young’s modulus |
| nu | float | 0.3 | Poisson’s ratio |
| d1 | float | 8e-05 | kg/(ms) Frequency independent loss |
| d3 | float | 1.4e-05 | kg m/s Frequency dependent loss |
| Ts0 | float | 2620 | N/m Tension per unit length |
| n_max_modes | int | 36 | Number of modal coordinates |
| use_nonlinear | bool | True | Use nonlinear wave equation |
# Initialize the solver
solver = Wave2dSolverTensionModulated(sampling_rate=16000, final_time=0.005, n_max_modes=625,n_gridpoints_x=161, use_nonlinear=True)
# Plot some of the modes
nx = 3
ny = 2
fig, ax = plt.subplots(nx, ny, figsize=(5, 15))
aspect_ratio_dxdy = solver.dx/solver.dy
# Print this aspect ratio and the one used in the solver
solver.print_solver_info()
print(f"Aspect ratio: {aspect_ratio_dxdy}")
print(f"Solver aspect ratio: {solver.aspect_ratio}")
for i in range(nx):
for j in range(ny):
print(f"Mode max value: {np.max(solver.modes[i, j, :, :])}")
ax[i, j].imshow(solver.modes[i, j, :, :], aspect = aspect_ratio_dxdy, origin='lower')
ax[i, j].set_title(f"Mode nx={i+1}, ny={j+1}")
ax[i, j].set_xlabel("y")
ax[i, j].set_ylabel("x")
# ax[i, j].tight_layout()dx: 0.0025 in meters
dy: 0.0025000000000000005 in meters
dt: 6.25e-05 in seconds
number of points in the x direction (n_gridpoints_x): 161
number of points in the y direction (n_gridpoints_y): 129
time in samples (nt): (80,)
number of modes in x direction (n_max_modes_x): 25
number of modes in y direction (n_max_modes_y): 25
number of modes (n_max_modes): 625
length in x direction (length_x): 0.4 in meters
length in y direction (length_y): 0.32000000000000006 in meters
grid_x shape: (161, 129)
grid_y shape: (161, 129)
wavenumbers_x shape: (25,)
wavenumbers_y shape: (25,)
modes shape: (25, 25, 161, 129)
lambdas shape: (625,)
Aspect ratio: 0.9999999999999998
Solver aspect ratio: 0.8
Mode max value: 1.0
Mode max value: 1.0
Mode max value: 1.0
Mode max value: 1.0
Mode max value: 0.9998072404820648
Mode max value: 1.0
# Create a simple initial condition, one the modes + additive noise
mode_n = 3
mode_m = 3
indn = mode_n - 1
indm = mode_m - 1
u0 = 0.0*solver.modes[indn, indm, :, :] + 0.001 * np.random.randn(solver.n_gridpoints_x, solver.n_gridpoints_y)
v0 = np.zeros_like(u0)# This was code for checking the projection to modal coordinates depending on the integration method
u0 = solver.modes[indn, indm, :, :]
v0 = np.zeros_like(u0)
# Project the initial condition to modal coordinates
bar_u_simpson, bar_v_simpson = solver.to_modal(u0,v0)
# Print the modal coordinates
print("Difference")
# Print the difference between the two methods
# print(bar_z - bar_z2)
# f = lambda y, x: np.sin(solver.wavenumbers_x[indn] * x) * np.sin(solver.wavenumbers_y[indm] * y)*np.sin(solver.wavenumbers_x[indn] * x) * np.sin(solver.wavenumbers_y[indm] * y)
# for i, x in enumerate(solver.x):
# for j, y in enumerate(solver.y):
# f_val = f(y, x)
# difference = f_val-solver.modes[indn, indm, i, j]
# if difference > 1e-15:
# print(f"Error at ({x}, {y})")
# print(f"Exact: {f_val}")
# print(f"Approx: {solver.modes[indn, indm, i, j]}")
# print(f"Difference: {difference}")
# print(solver.length_x)
# print(solver.length_y)
# print(solver.wavenumbers_x[indn])
# print(solver.wavenumbers_y[indm])
# exact_int = dblquad(f, 0, 1.0, 0, 0.5)
# print((exact_int))Differencedef gaussian_pulse(x_grid, y_grid, x0, y0, sigma):
return np.exp(-((x_grid-x0)**2 + (y_grid-y0)**2)/(2*sigma**2))
def noise(grid):
return np.random.randn(grid.shape)# Solve the wave equation
# u0 = np.zeros((solver.n_gridpoints_x, solver.n_gridpoints_y))
u0 = (solver.modes[0, 0, :, :]**2)
# Create a gaussian impulse
u0 = np.zeros((solver.n_gridpoints_x, solver.n_gridpoints_y))
u0 = np.random.randn(solver.n_gridpoints_x, solver.n_gridpoints_y)
ctr = (0.37*solver.length_x, 0.63*solver.length_y)
std = 0.005
for i in range(solver.n_gridpoints_x):
for j in range(solver.n_gridpoints_y):
u0[i, j] = np.exp(-((solver.grid_x[i, j] - ctr[0])**2 + (solver.grid_y[i, j] - ctr[1])**2)/(2*std**2))
u0_alt = gaussian_pulse(solver.grid_x, solver.grid_y, ctr[0], ctr[1], std)
# Compare the two initial conditions
print(np.allclose(u0, u0_alt))
v0 = np.zeros_like(u0)
u0 = 0.005*u0
# Add a delta impuse to the initial velocity, close to the center
# v0[solver.n_gridpoints_x//2, solver.n_gridpoints_y//2] = 20
t, u, v = solver.solve(u0=u0, v0=v0)
print(f"u shape: {u.shape}")
print(f"v shape: {v.shape}")True
bar_u shape: (625, 80)
bar_v shape: (625, 80)
u shape: (80, 161, 129)
v shape: (80, 161, 129)# Check that the boundary conditions are maintained
print(np.all(u[:, 0, :] == 0))
print(np.all(u[:, -1, :] == 0))
print(np.all(u[:, :, 0] == 0))
print(np.all(u[:, :, -1] == 0))
# Plot the sum error of the solution at the edges for all timesteps
plt.figure()
plt.plot(t, np.sum(u[:, 0, :], axis=1), label="bottom")
plt.plot(t, np.sum(u[:, -1, :], axis=1), label="top")
plt.plot(t, np.sum(u[:, :, 0], axis=1), label="left")
plt.plot(t, np.sum(u[:, :, -1], axis=1), label="right")
plt.legend()
plt.xlabel("Time")
plt.ylabel("Sum of displacement at the edges")
plt.title("Boundary conditions")
plt.show()True
False
True
False
# Plot the solution as an animation usign matplotlib
# N_plots=100
fig,ax = plt.subplots(1,2)
def animate(i):
ax[0].clear()
ax[1].clear()
ax[0].imshow(u[i,:,:], aspect = aspect_ratio_dxdy, origin='lower')
ax[1].imshow(v[i,:,:], aspect = aspect_ratio_dxdy, origin='lower')
ax[0].set_title(f"Displacement")
ax[0].set_xlabel("y")
ax[0].set_ylabel("x")
ax[1].set_title(f"Velocity")
ax[1].set_xlabel("y")
ax[1].set_ylabel("x")
ani = FuncAnimation(fig, animate, frames=20,
interval=500, repeat=False)
plt.close()
from IPython.display import HTML
HTML(ani.to_jshtml())