from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
import time
Tension modulated stiff membrane
This notebook is a solver using the functional transform method for a stiff membrane with tension modulation. The membrane is described by the following PDE:
\[ \begin{align} D \nabla^4 u +\rho h \frac{\partial^2 u}{\partial t^2} - T(u) \nabla^2 u + d_{1} \frac{\partial u}{\partial t} + d_{3} \frac{\partial \nabla^2 u}{\partial t} = f^{(ext)} \end{align} \] Assumming bending stiffness to be small \(D/T \ll 1\) Where for a rectangular plate \(x \in [0, L_x], y \in [0, L_y]\) of thickness \(h\) , \(D = Q h^3 / 12 (1 - \nu ^2)\) \(\rho\) is the mass density \((kg/m^3)\), T(u) is the tension, \(d_1\) and \(d_3\) are the damping coefficients, and \(f^{(ext)}\) is the external force.
The modes for simply supported BC are \[ K_{n,m}\left(x, y \right) =\sin \left(\frac{n \pi x}{L_x} \right) sin \left(\frac{m \pi y}{L_y} \right) \] Where \(n, m\) are integers. \[ \begin{aligned} & \nabla^2 K_{n, m}(x, y)=-\lambda_{n, m} K_{n, m}(x, y), \\ & \text { with } \quad \lambda_{n, m}=\pi^2\left[\left(\frac{n}{L_x}\right)^2+\left(\frac{m}{L_y}\right)^2\right] . \end{aligned} \]
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} \]
::: {#cell-9 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
import numpy as np
from scipy.integrate import solve_ivp, simpson
:::
::: {#cell-11 .cell 0=‘e’ 1=‘x’ 2=‘p’ 3=‘o’ 4=‘r’ 5=‘t’}
class Wave2dSolverTensionModulated:
"""
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.
"""
def __init__(
self,
int = 16000, # 1/s Temporal sampling frequency
sampling_rate: float = 0.5, # s Duration of the simulation
final_time: int = 41, # pts/m Spatial sampling grid
n_gridpoints_x: float = 0.4, # m Length of x dimension
length_x: float = 0.9, # Aspect ratio of the membrane, Ly/Lx
aspect_ratio: float = 1380, # kg/m**3 Density
rho: float = 1.9e-4, # m Thickness
h: int = 3.5e9, # Pa Young's modulus
E: float = 0.3, # Poisson's ratio
nu: float = 8e-5, # kg/(ms) Frequency independent loss
d1: float = 1.4e-5, # kg m/s Frequency dependent loss
d3: float = 2620, # N/m Tension per unit length
Ts0: int = 36, # Number of modal coordinates
n_max_modes: bool = True, # Use nonlinear wave equation
use_nonlinear: bool = False, # Calculate the energy of the system
include_energy:
):# TODO: use a ratio for the side lengths, and use the same ratio for gridpoints
# Attributes intrinsic to the Wave Equation PDE
self.pde_num_variables = 1
self.pde_num_spatial_dimensions = 2
self.pde_order_time_derivatives = 2
# Attributes of the simulation
self.sampling_rate = sampling_rate
self.final_time = final_time
self.n_gridpoints_x = n_gridpoints_x
self.length_x = length_x
self.aspect_ratio = aspect_ratio
self.use_nonlinear = use_nonlinear
self.include_energy = include_energy
# Attributes of the membrane
self.rho = rho
self.h = h
self.E = E
self.nu = nu
self.d1 = d1
self.d3 = d3
self.Ts0 = Ts0
# Use the same number of modes in x and y, this might not be a good idea for aspect ratios far from 1
self.n_max_modes_x = int(np.floor(np.sqrt(n_max_modes)))
self.n_max_modes_y = int(np.floor(np.sqrt(n_max_modes)))
self.n_max_modes = self.n_max_modes_x * self.n_max_modes_y
# if n_max_modes > n_gridpoints - 2:
# self.n_max_modes = n_gridpoints - 2
# print(f"n_max_modes too high, setting to {self.n_max_modes}")
# else:
# self.n_max_modes = n_max_modes
= np.linspace(0, self.length_x, self.n_gridpoints_x)
x self.dx = x[1] - x[0] # m spatial sampling interval
# calculate the gridpoints and length in y using the aspect ratio
# We want dy to be as close as possible to dx
self.length_y = self.length_x * self.aspect_ratio
self.n_gridpoints_y = (
int(np.floor((self.n_gridpoints_x - 1) * self.aspect_ratio)) + 1
)
= np.linspace(0, self.length_y, self.n_gridpoints_y)
y self.x = x
self.y = y
self.grid_x, self.grid_y = np.meshgrid(x, y, indexing="ij")
# spatial grid
# self.grid_x = np.m
self.dy = y[1] - y[0] # m spatial sampling interval
# temporal grid use arange to make sure that the timestep corresponds exactly to the sampling frequency
self.dt = 1 / self.sampling_rate # s temporal sampling interval
self.timesteps = np.arange(0, self.final_time, self.dt)
self.mu_x = np.arange(1, self.n_max_modes_x + 1)
self.mu_y = np.arange(1, self.n_max_modes_y + 1)
self.wavenumbers_x = self.mu_x * np.pi / self.length_x
self.wavenumbers_y = self.mu_y * np.pi / self.length_y
self.grid_wavenumber_x, self.grid_wavenumber_y = np.meshgrid(
self.wavenumbers_x, self.wavenumbers_y
)self.modes = np.zeros(
(self.n_max_modes_x,
self.n_max_modes_y,
self.n_gridpoints_x,
self.n_gridpoints_y,
)
)self.lambdas = np.zeros((self.n_max_modes_x * self.n_max_modes_y))
# modes and eigenvalues of the modes (lanbdas)
for i, wx in enumerate(self.wavenumbers_x):
for j, wy in enumerate(self.wavenumbers_y):
self.modes[i, j, :, :] = np.sin(wx * self.grid_x) * np.sin(
* self.grid_y
wy
)self.lambdas[i * self.n_max_modes_y + j] = wx**2 + wy**2
self.norm_factor = 4 / (
self.length_x * self.length_y
# normalization factor for the modes
) self.D = self.E * (self.h) ** 3 / (12 * (1 - self.nu**2))
= self.E * self.h / (1 - self.nu**2)
CNL
= self.D * self.lambdas**2 + self.Ts0 * self.lambdas
beta_mu
# calculate the matrices
self.Lambdadiag = np.diag(self.lambdas)
# self.H_1 = np.diag(self.lambdas) / (self.rho * self.h)
self.M_v = np.diag(self.d1 + self.d3 * self.lambdas) / (self.rho * self.h)
self.M_u = np.diag(beta_mu) / (self.rho * self.h)
# coeffficient for the nonlinear term
self.Cb = (
2
* self.E
/ (self.rho * self.length_x**2 * self.length_y**2 * (1 - self.nu**2))
)
def print_matrices(self):
print(f"Lambdadiag shape: {self.Lambdadiag.shape}")
# print(f"H_1 shape: {self.H_1.shape}")
print(f"M_v shape: {self.M_v.shape}")
print(f"M_u shape: {self.M_u.shape}")
return
def print_solver_info(self):
# Print some information
print(f"dx: {self.dx} in meters")
print(f"dy: {self.dy} in meters")
print(f"dt: {self.dt} in seconds")
print(
f"number of points in the x direction (n_gridpoints_x): {self.n_gridpoints_x}"
)print(
f"number of points in the y direction (n_gridpoints_y): {self.n_gridpoints_y}"
)print(f"time in samples (nt): {self.timesteps.shape}")
print(f"number of modes in x direction (n_max_modes_x): {self.n_max_modes_x}")
print(f"number of modes in y direction (n_max_modes_y): {self.n_max_modes_y}")
print(
f"number of modes (n_max_modes): {self.n_max_modes_x * self.n_max_modes_y}"
)print(f"length in x direction (length_x): {self.length_x} in meters")
print(f"length in y direction (length_y): {self.length_y} in meters")
# Print the shapes of the grids and the wavenumbers
print(f"grid_x shape: {self.grid_x.shape}")
print(f"grid_y shape: {self.grid_y.shape}")
print(f"wavenumbers_x shape: {self.wavenumbers_x.shape}")
print(f"wavenumbers_y shape: {self.wavenumbers_y.shape}")
print(f"modes shape: {self.modes.shape}")
print(f"lambdas shape: {self.lambdas.shape}")
return
def to_modal(
self,
u,
v,str = "simpson", # displacement # velocity
integrator: -> tuple[np.ndarray, np.ndarray]:
) """Project the displacement and velocity to modal coordinates.
Also flatten the arrays to be 1D.
"""
= np.zeros(self.n_max_modes)
bar_u # bar_z2 = np.zeros(self.n_max_modes)
= np.zeros(self.n_max_modes)
bar_v
if integrator == "simpson":
for i in range(self.n_max_modes_x):
for j in range(self.n_max_modes_y):
= u * self.modes[i, j, :, :]
uu * self.n_max_modes_y + j] = simpson(
bar_u[i =self.y) for uu_y in uu], x=self.x
[simpson(uu_y, x
)= v * self.modes[i, j, :, :]
vv * self.n_max_modes_y + j] = simpson(
bar_v[i =self.y) for vv_y in vv], x=self.x
[simpson(vv_y, x
)elif integrator == "trapz":
# This is unverified, use simpson for now
raise NotImplementedError
# for i in range(self.n_max_modes_x):
# for j in range(self.n_max_modes_y):
# bar_u[i * self.n_max_modes_y + j] = self.dx*self.dy*np.sum(self.modes[i, j, :, :] * u)
# # bar_z_dot[i * self.n_max_modes_y + j] = np.sum(
# self.modes[i, j, :, :] * z_dot
# )
else:
raise ValueError(f"Integrator {integrator} not recognised")
return bar_u, bar_v
def to_modal_vecs(
self,
u,
v,str = "simpson", # displacement # velocity
integrator: -> tuple[np.ndarray, np.ndarray]:
) """Project the displacement and velocity to modal coordinates.
Also flatten the arrays to be 1D.
"""
= np.zeros((u.shape[0], self.n_max_modes))
bar_u = np.zeros((u.shape[0], self.n_max_modes))
bar_v if integrator == "simpson":
for i in range(self.n_max_modes_x):
for j in range(self.n_max_modes_y):
= u * self.modes[i, j, :, :]
uu * self.n_max_modes_y + j] = simpson(
bar_u[:, i =self.y) for uu_y in uu], x=self.x
[simpson(uu_y, x
)= v * self.modes[i, j, :, :]
vv * self.n_max_modes_y + j] = simpson(
bar_v[:, i =self.y) for vv_y in vv], x=self.x
[simpson(vv_y, x
)elif integrator == "trapz":
raise NotImplementedError
else:
raise ValueError(f"Integrator {integrator} not recognised")
return bar_u, bar_v
def to_displacement(
self,
# modal displacement
bar_u, # modal velocity
bar_v, -> tuple[np.ndarray, np.ndarray]:
) """Sum the modal displacements and velocities to get the displacement and velocity"""
= np.zeros((self.n_gridpoints_x, self.n_gridpoints_y))
u = np.zeros((self.n_gridpoints_x, self.n_gridpoints_y))
v = self.norm_factor
norm_factor for i in range(self.n_max_modes_x):
for j in range(self.n_max_modes_y):
+= bar_u[i * self.n_max_modes_y + j] * self.modes[i, j, :, :]
u += bar_v[i * self.n_max_modes_y + j] * self.modes[i, j, :, :]
v *= norm_factor
u *= norm_factor
v return u, v
def enforce_boundary_conditions(self, u, v):
"""Enforce dirichlet boundary conditions.
ATTN: remember v is respect to time, not x,y
"""
0, :] = 0
u[-1, :] = 0
u[0] = 0
u[:, -1] = 0
u[:, 0, :] = 0
v[-1, :] = 0
v[0] = 0
v[:, -1] = 0
v[:, return u, v
def nonlinear_membrane(
self,
# time
t, # state vector
state, -> np.ndarray: # state at timestep t and position x (u(x, t))
) # unpack the state vector
= state[: self.n_max_modes] # displacement in modal coordinates
bar_u = state[self.n_max_modes :] # velocity in modal coordinates
bar_v
# add axis for calculation
= bar_u[..., None]
bar_u
# calculate the $\bar{b}$ vector
# ATTN: this b doesn't include the minus sign
if self.use_nonlinear:
= bar_u.T @ self.Lambdadiag @ bar_u
deformation = self.Cb * deformation * self.Lambdadiag @ bar_u
b
# update the state vector
# Why .squeeze()? Is it really needed?
= bar_v
bar_udot = -self.M_v @ bar_v - self.M_u @ bar_u.squeeze()
bar_vdot if self.use_nonlinear:
-= b.squeeze()
bar_vdot
# return the state derivatives
return np.concatenate([bar_udot, bar_vdot])
def calculate_energy(
self,
# modal coordinates of the displacement
bar_u: np.ndarray, # modal coordinates of the velocity
bar_v: np.ndarray, -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
) """Calculate the kinetic, elastic potential, stiffness potential and nonlinear potential energy of the system.
Args:
bar_u (np.ndarray): modal coordinates of the displacement. Shape (nt, n_max_modes)
bar_v (np.ndarray): modal coordinates of the velocity. Shape (nt, n_max_modes)
Returns:
np.ndarray: _description_
"""
# calculate the kinetic energy
= (
energy_kinetic 0.5 * self.rho * self.h * self.norm_factor * np.sum(bar_v**2, axis=-1)
)# calculate the potential energy
= bar_u**2
bar_uu = np.sum(bar_uu * self.lambdas, axis=-1)
summation = 0.5 * self.Ts0 * self.norm_factor * summation
energy_potential_elastic_linear = (
energy_potential_stiffness_linear 0.5 * self.D * self.norm_factor * np.sum(bar_uu * self.lambdas**2, axis=-1)
)if self.use_nonlinear:
= (
TNL
summation* self.norm_factor
* self.E
* self.h
/ (2 * (1 - self.nu**2) * self.length_x * self.length_y)
)= (
energy_potential_elastic_nonlinear 0.5 * TNL * self.norm_factor * summation
)else:
= np.zeros_like(
energy_potential_elastic_nonlinear
energy_potential_elastic_linear
)
return (
energy_kinetic,
energy_potential_elastic_linear,
energy_potential_stiffness_linear,
energy_potential_elastic_nonlinear,
)
def solve(
self,
= None, # initial displacement (default: None)
u0: np.ndarray = None, # initial velocity (default: None)
v0: np.ndarray -> tuple[np.ndarray, np.ndarray, np.ndarray]:
) """Solve the wave equation with the given initial conditions."""
# set the initial conditions if not given
= u0 if u0 is not None else np.zeros_like(self.grid_x)
u0 = v0 if v0 is not None else np.zeros_like(self.grid_x)
v0
# Enforce the simply suppported boundary conditions a tht e edges of the membrane
= self.enforce_boundary_conditions(u0, v0)
u0, v0 # transform the initial conditions to modal coordinates
= self.to_modal(u0, v0)
bar_u0, bar_v0 # solve the wave equation in modal coordinates
= solve_ivp(
sol =self.nonlinear_membrane,
fun=[0, self.final_time],
t_span=np.concatenate([bar_u0, bar_v0], axis=0),
y0=self.timesteps,
t_eval="DOP853",
method=1e-12,
rtol=1e-14,
atol
)# sol has dimensions (n_max_modes*2, nt)
# unpack the solution
= (sol.y[: self.n_max_modes]).T
bar_u = (sol.y[self.n_max_modes :]).T
bar_v
# Print the shapes of the solutions
print(f"bar_u shape: {bar_u.shape}")
print(f"bar_v shape: {bar_v.shape}")
# transform back to the physical domain
# This loop is probably extremely slow but is only done once
= np.zeros(
u
(len(sol.t),
self.n_gridpoints_x,
self.n_gridpoints_y,
)
)= np.zeros(
v
(len(sol.t),
self.n_gridpoints_x,
self.n_gridpoints_y,
)
)for i in range(len(sol.t)):
= self.to_displacement(bar_u[i, :], bar_v[i, :])
u[i, :, :], v[i, :, :]
# calculate the energy of the system
if self.include_energy:
(
energy_kinetic,
energy_potential_elastic_linear,
energy_potential_stiffness_linear,
energy_potential_elastic_nonlinear,= self.calculate_energy(bar_u, bar_v)
) return (
sol.t,
u,
v,
energy_kinetic,
energy_potential_elastic_linear,
energy_potential_stiffness_linear,
energy_potential_elastic_nonlinear,
)return (
sol.t,
u,
v, )
:::
= 16000
sampling_rate # Initialize the solver
= Wave2dSolverTensionModulated(
solver =sampling_rate,
sampling_rate=0.1,
final_time=225,
n_max_modes=41,
n_gridpoints_x=800,
Ts0=0.5,
d1=5e-3,
d3=True,
use_nonlinear=True,
include_energy
)# Plot some of the modes
= 3
nx = 2
ny = plt.subplots(nx, ny, figsize=(5, 15))
fig, ax = solver.dx / solver.dy
aspect_ratio_dxdy # 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(=aspect_ratio_dxdy, origin="lower"
solver.modes[i, j, :, :], aspect
)f"Mode nx={i + 1}, ny={j + 1}")
ax[i, j].set_title("y")
ax[i, j].set_xlabel("x")
ax[i, j].set_ylabel(# ax[i, j].tight_layout()
# Compare the performance of to_modal and to_modal_vecs
# Generate a tensor of states
= 10
n_states = np.random.randn(n_states, solver.n_gridpoints_x, solver.n_gridpoints_y)
u = np.random.randn(n_states, solver.n_gridpoints_x, solver.n_gridpoints_y)
v
# Vectorised version
= time.time()
start = solver.to_modal_vecs(u, v)
bar_u_vec, bar_v_vec = time.time()
end print(f"Vectorised version took {end - start} seconds")
# Non-vectorised version
= time.time()
start = np.zeros((n_states, solver.n_max_modes))
bar_u = np.zeros((n_states, solver.n_max_modes))
bar_v for i in range(n_states):
= solver.to_modal(u[i], v[i])
bar_u[i], bar_v[i] = time.time()
end print(f"Non-vectorised version took {end - start} seconds")
# Compare the results
print(f"bar_u_vec shape: {bar_u_vec.shape}")
print(f"bar_u shape: {bar_u.shape}")
print(f"bar_v_vec shape: {bar_v_vec.shape}")
print(f"bar_v shape: {bar_v.shape}")
print(f"bar_u_vec - bar_u: {np.max(np.abs(bar_u_vec - bar_u))}")
print(f"bar_v_vec - bar_v: {np.max(np.abs(bar_v_vec - bar_v))}")
def 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))
= solver.modes[0, 0, :, :]
u0 # Create a gaussian impulse
# u0 = np.zeros((solver.n_gridpoints_x, solver.n_gridpoints_y))
= np.random.randn(solver.n_gridpoints_x, solver.n_gridpoints_y)
u0 = (0.2 * solver.length_x, 0.6 * solver.length_y)
ctr = 0.05 * solver.length_x
std for i in range(solver.n_gridpoints_x):
for j in range(solver.n_gridpoints_y):
= np.exp(
u0[i, j] -((solver.grid_x[i, j] - ctr[0]) ** 2 + (solver.grid_y[i, j] - ctr[1]) ** 2)
/ (2 * std**2)
)
= gaussian_pulse(solver.grid_x, solver.grid_y, ctr[0], ctr[1], std)
u0_alt # Compare the two initial conditions
print(np.allclose(u0, u0_alt))
# v0 = np.zeros_like(u0)
# u0 = 0.01*u0
= 20 * u0
v0 = np.zeros_like(u0)
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
= solver.solve(u0=u0, v0=v0)
t, u, v, e_k, e_ple, e_pls, e_pne print(f"t shape: {t.shape}")
print(f"u shape: {u.shape}")
print(f"v shape: {v.shape}")
# 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()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.plot(t, np.
plt.legend()"Time")
plt.xlabel("Sum of displacement at the edges")
plt.ylabel("Boundary conditions")
plt.title( plt.show()
= np.max(np.abs(u))
vmax_u = np.max(np.abs(v))
vmax_v
# N_plots=100
= plt.subplots(1, 2)
fig, ax
def animate(i):
0].clear()
ax[1].clear()
ax[0].imshow(
ax[=aspect_ratio_dxdy, origin="lower", vmin=-vmax_u, vmax=vmax_u
u[i, :, :], aspect
)1].imshow(
ax[=aspect_ratio_dxdy, origin="lower", vmin=-vmax_v, vmax=vmax_v
v[i, :, :], aspect
)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")
ax[
5)
fig.set_figheight(10)
fig.set_figwidth(= FuncAnimation(fig, animate, frames=range(0, 320, 50), interval=300, repeat=False)
ani
plt.close()from IPython.display import HTML
HTML(ani.to_jshtml())
# Max amplitude after t_late seconds of the simulation
= 0.04
t_late = int(t_late * sampling_rate)
sample_ind = np.max(np.abs(u[sample_ind, :, :]))
max_amplitude print(f"Max amplitude after {t_late}s : {max_amplitude}")
# Max speed in the last 1 seconds of the simulation
= np.max(np.abs(v[sample_ind, :, :]))
max_speed print(f"Max speed after {t_late}s : {max_speed}")
# Plot state at the beginning and after 2 seconds of the simulation
= plt.subplots(2, 2)
fig, ax 0, 0].imshow(
ax[0, :, :], aspect=aspect_ratio_dxdy, origin="lower", vmin=-vmax_u, vmax=vmax_u
u[
)0, 0].set_title("Initial displacement")
ax[0, 1].imshow(
ax[0, :, :], aspect=aspect_ratio_dxdy, origin="lower", vmin=-vmax_v, vmax=vmax_v
v[
)0, 1].set_title("Initial velocity")
ax[1, 0].imshow(
ax[
u[sample_ind, :, :],=aspect_ratio_dxdy,
aspect="lower",
origin=-vmax_u,
vmin=vmax_u,
vmax
)1, 0].set_title("Final displacement")
ax[1, 1].imshow(
ax[
v[sample_ind, :, :],=aspect_ratio_dxdy,
aspect="lower",
origin=-vmax_v,
vmin=vmax_v,
vmax
)1, 1].set_title("Final velocity")
ax[ plt.show()
# PLot the spectrogram for the output of a single point o the membrane
import matplotlib.pyplot as plt
from scipy.signal import spectrogram
import IPython.display as ipd
# print the shape of u and v
print("u shape:", u.shape)
print("v shape:", v.shape)
# Print the maximum value of u and v
print("Max u:", np.max(np.abs(u)))
print("Max v:", np.max(np.abs(v)))
= (int(u.shape[1] * np.random.random()), int(u.shape[2] * np.random.random()))
out_point # Check the that the output point is within the bounds of the data
assert out_point[0] < u.shape[1]
assert out_point[1] < u.shape[2]
= u[:, out_point[0], out_point[1]]
u_out
# plot the time domain output of a single point on the membrane
plt.figure()
plt.plot(u_out)f"u at point {out_point}")
plt.title(
plt.figure()
plt.show()
=sampling_rate))
display(ipd.Audio(u_out, rate
# PLot the spectrogram for the output of a single point on the membrane
= 512
win_size = win_size // 8
win_over = spectrogram(
f, time_spec, Sxx
u_out,=sampling_rate,
fs=win_size,
nperseg=win_size,
nfft=win_size - win_over,
noverlap
)= Sxx / Sxx.max()
Sxx_norm = 10 * np.log10(Sxx_norm)
Sxx_norm_db # Make all values below -60 dB to -1000 dB
< -60] = -60
Sxx_norm_db[Sxx_norm_db
plt.figure()
plt.pcolormesh(time_spec, f, Sxx_norm_db)"Frequency [Hz]")
plt.ylabel(
plt.semilogy()# plt.yscale('log')
100, 1000)
plt.ylim(
plt.yticks(10, 100, 200, 300, 400, 600, 800, 2000], [10, 100, 200, 300, 400, 600, 800, 2000]
[
)# plt.xlim(0, 0.25)
"Time [sec]")
plt.xlabel(f"Spectrogram of u at point {out_point}")
plt.title(="dB")
plt.colorbar(label plt.show()
from scipy.sparse import diags, kron, eye
# Functions for calculating the energy of the system
def calculate_energy_kinetic(solver, v):
= v**2
vv = simpson([simpson(vv_y, x=solver.y) for vv_y in vv], x=solver.x)
integral = solver.dx * solver.dy * np.sum(vv)
integral_trapz # print(f"Integral simpson: {integral}")
# print(f"Integral trapz: {integral_trapz}")
return solver.rho * solver.h * integral / 2
def calculate_energy_potential_elastic_linear(solver, u):
= solver.n_gridpoints_x
Nx = solver.n_gridpoints_y
Ny = solver.dx
dx = solver.dy
dy
# create finite difference matrices
# Dxx
= diags(
Dxx 1, -2, 1], [-1, 0, 1], shape=(Nx, Nx)
[# Create sparse matrix Dxx
).toarray() /= dx**2 # Divide matrix Dxx by h**2
Dxx # Dyy
= diags(
Dyy 1, -2, 1], [-1, 0, 1], shape=(Ny, Ny)
[# Create sparse matrix Dyy
).toarray() /= dy**2 # Divide matrix Dyy by h**2
Dyy # Implement dxxu0 = 0 boundary conditions (second order, centered)
0, 0] = 0
Dxx[-1, 0] = 0
Dxx[0, 0] = 0
Dyy[-1, 0] = 0
Dyy[
# Implement dirichlet boundary conditions
0, :] = 0
Dxx[-1, :] = 0
Dxx[0, :] = 0
Dyy[-1, :] = 0
Dyy[= kron(Dxx, eye(Ny))
Dxx2 = kron(eye(Nx), Dyy)
Dyy2 = Dxx2 + Dyy2
nabla2
# Flatten the displacement
= u.reshape(-1, 1)
u_flat
= (nabla2 @ u_flat) * u_flat
integrand_flat = integrand_flat.reshape(Nx, Ny)
integrand = simpson(
integral =solver.y) for integrand_y in integrand], x=solver.x
[simpson(integrand_y, x
)return -solver.Ts0 * integral / 2
def calculate_energy_potential_stiffness_linear(solver, u):
= solver.n_gridpoints_x
Nx = solver.n_gridpoints_y
Ny = solver.dx
dx = solver.dy
dy
# create finite difference matrices
# Dxx
= diags(
Dxx 1, -2, 1], [-1, 0, 1], shape=(Nx, Nx)
[# Create sparse matrix Dxx
).toarray() /= dx**2 # Divide matrix Dxx by h**2
Dxx # Dyy
= diags(
Dyy 1, -2, 1], [-1, 0, 1], shape=(Ny, Ny)
[# Create sparse matrix Dyy
).toarray() /= dy**2 # Divide matrix Dyy by h**2
Dyy # Implement dxxu0 = 0 boundary conditions (second order, centered)
0, 0] = 0
Dxx[-1, 0] = 0
Dxx[0, 0] = 0
Dyy[-1, 0] = 0
Dyy[
# Implement dirichlet boundary conditions
0, :] = 0
Dxx[-1, :] = 0
Dxx[0, :] = 0
Dyy[-1, :] = 0
Dyy[= kron(Dxx, eye(Ny))
Dxx2 = kron(eye(Nx), Dyy)
Dyy2 = Dxx2 + Dyy2
nabla2 = nabla2 @ nabla2
nabla4
# Flatten the displacement
= u.reshape(-1, 1)
u_flat
= (nabla4 @ u_flat) * u_flat
integrand_flat = integrand_flat.reshape(Nx, Ny)
integrand = simpson(
integral =solver.y) for integrand_y in integrand], x=solver.x
[simpson(integrand_y, x
)return solver.D * integral / 2
def calculate_energy_potential_elastic_nonlinear(solver, u):
= solver.n_gridpoints_x
Nx = solver.n_gridpoints_y
Ny = solver.dx
dx = solver.dy
dy
# create finite difference matrices
# Dxx
= diags(
Dx -1, 0, 1], [-1, 0, 1], shape=(Nx, Nx)
[# Create sparse matrix Dxx
).toarray() /= 2 * dx # Divide matrix Dxx by h**2
Dx # Dyy
= diags(
Dy -1, 0, 1], [-1, 0, 1], shape=(Ny, Ny)
[# Create sparse matrix Dyy
).toarray() /= 2 * dy # Divide matrix Dyy by h**2
Dy # Implement dxxu0 = 0 boundary conditions (second order, centered)
0, 0] = 0
Dx[-1, 0] = 0
Dx[0, 0] = 0
Dy[-1, 0] = 0
Dy[
# Implement dirichlet boundary conditions
0, :] = 0
Dx[-1, :] = 0
Dx[0, :] = 0
Dy[-1, :] = 0
Dy[= kron(Dx, eye(Ny))
Dx2 = kron(eye(Nx), Dy) # Flatten the displacement
Dy2 = u.reshape(-1, 1)
u_flat
# Calculate the gradient of the displacement
= (Dx2 @ u_flat) ** 2
Dx2u_sq = (Dy2 @ u_flat) ** 2
Dy2u_sq
= Dx2u_sq + Dy2u_sq
integrand_flat_T = integrand_flat_T.reshape(Nx, Ny)
integrand_T = simpson(
integral_T =solver.y) for integrand_y in integrand_T], x=solver.x
[simpson(integrand_y, x
)
= (
TNL
solver.E* solver.h
/ (2 * (1 - solver.nu**2) * solver.length_x * solver.length_y)
* integral_T
)
# create finite difference matrices
# Dxx
= diags(
Dxx 1, -2, 1], [-1, 0, 1], shape=(Nx, Nx)
[# Create sparse matrix Dxx
).toarray() /= dx**2 # Divide matrix Dxx by h**2
Dxx # Dyy
= diags(
Dyy 1, -2, 1], [-1, 0, 1], shape=(Ny, Ny)
[# Create sparse matrix Dyy
).toarray() /= dy**2 # Divide matrix Dyy by h**2
Dyy # Implement dxxu0 = 0 boundary conditions (second order, centered)
0, 0] = 0
Dxx[-1, 0] = 0
Dxx[0, 0] = 0
Dyy[-1, 0] = 0
Dyy[
# Implement dirichlet boundary conditions
0, :] = 0
Dxx[-1, :] = 0
Dxx[0, :] = 0
Dyy[-1, :] = 0
Dyy[= kron(Dxx, eye(Ny))
Dxx2 = kron(eye(Nx), Dyy)
Dyy2 = Dxx2 + Dyy2
nabla2
= (nabla2 @ u_flat) * u_flat
integrand_flat = integrand_flat.reshape(Nx, Ny)
integrand = simpson(
integral =solver.y) for integrand_y in integrand], x=solver.x
[simpson(integrand_y, x
)
return -TNL * integral / 2
# Calculate the energy of the system over time, plot the kinetic, linear potential and total energy
= len(t) - 1
block_size = 0
start_t_samples = t[start_t_samples : start_t_samples + block_size]
t_energy = np.zeros_like(t_energy)
energy_kinetic = np.zeros_like(t_energy)
energy_potential_elastic_linear = np.zeros_like(t_energy)
energy_potential_stiffness_linear = np.zeros_like(t_energy)
energy_potential_elastic_nonlinear = np.zeros_like(t_energy)
energy_total for i in range(block_size):
= calculate_energy_kinetic(solver, v[start_t_samples + i, :, :])
energy_kinetic[i] = calculate_energy_potential_elastic_linear(
energy_potential_elastic_linear[i] + i, :, :]
solver, u[start_t_samples
)= calculate_energy_potential_stiffness_linear(
energy_potential_stiffness_linear[i] + i, :, :]
solver, u[start_t_samples
)= (
energy_potential_elastic_nonlinear[i]
calculate_energy_potential_elastic_nonlinear(+ i, :, :]
solver, u[start_t_samples
)
)
if solver.use_nonlinear:
= (
energy_total
energy_kinetic+ energy_potential_elastic_linear
+ energy_potential_stiffness_linear
+ energy_potential_elastic_nonlinear
)else:
= (
energy_total
energy_kinetic+ energy_potential_elastic_linear
+ energy_potential_stiffness_linear
)
= (
energy_linear + energy_potential_elastic_linear + energy_potential_stiffness_linear
energy_kinetic )
# Calculate the energy of the system using the modal coordinates
from matplotlib.pylab import solve
from matplotlib.pyplot import bar
def calculate_energy_kinetic_modal(solver, bar_v):
= bar_v**2
bar_vv = np.sum(bar_vv)
summation return solver.rho * solver.h * solver.norm_factor * summation / 2
def calculate_energy_potential_elastic_linear_modal(solver, bar_u):
= bar_u**2
bar_uu = np.sum(bar_uu * solver.lambdas)
summation return solver.Ts0 * solver.norm_factor * summation / 2
def calculate_energy_potential_stiffness_linear_modal(solver, bar_u):
= bar_u**2
bar_uu = np.sum(bar_uu * solver.lambdas**2)
summation return solver.D * solver.norm_factor * summation / 2
def calculate_energy_potential_elastic_nonlinear_modal(solver, bar_u):
= bar_u**2
bar_uu = np.sum(bar_uu * solver.lambdas)
summation = (
TNL
summation* solver.norm_factor
* solver.E
* solver.h
/ (2 * (1 - solver.nu**2) * solver.length_x * solver.length_y)
)return TNL * solver.norm_factor * summation / 2
# Calculate the energy of the system over time, plot the kinetic, linear potential and total energy
= np.zeros_like(t_energy)
energy_kinetic_modal = np.zeros_like(t_energy)
energy_potential_elastic_linear_modal = np.zeros_like(t_energy)
energy_potential_stiffness_linear_modal = np.zeros_like(t_energy)
energy_potential_elastic_nonlinear_modal
# Time these calculations
= time.time()
start_time
= solver.to_modal_vecs(
bar_u_vec, bar_v_vec + block_size],
u[start_t_samples : start_t_samples + block_size],
v[start_t_samples : start_t_samples
)= time.time()
end_time print(f"Time to calculate modal coordinates: {end_time - start_time}")
= time.time()
start_time for i in range(block_size):
= bar_u_vec[i], bar_v_vec[i]
bar_u, bar_v = calculate_energy_kinetic_modal(solver, bar_v)
energy_kinetic_modal[i] = (
energy_potential_elastic_linear_modal[i]
calculate_energy_potential_elastic_linear_modal(solver, bar_u)
)= (
energy_potential_stiffness_linear_modal[i]
calculate_energy_potential_stiffness_linear_modal(solver, bar_u)
)= (
energy_potential_elastic_nonlinear_modal[i]
calculate_energy_potential_elastic_nonlinear_modal(solver, bar_u)
)
= time.time()
end_time print(f"Time to calculate energy in modal coordinates: {end_time - start_time}")
if solver.use_nonlinear:
= (
energy_total_modal
energy_kinetic_modal+ energy_potential_elastic_linear_modal
+ energy_potential_stiffness_linear_modal
+ energy_potential_elastic_nonlinear_modal
)else:
= (
energy_total_modal
energy_kinetic_modal+ energy_potential_elastic_linear_modal
+ energy_potential_stiffness_linear_modal
)
= (
energy_linear_modal
energy_kinetic_modal+ energy_potential_elastic_linear_modal
+ energy_potential_stiffness_linear_modal
)
# Compare the energy of the system using the modal coordinates and the physical coordinates
# Kinetic energy
plt.figure()="Kinetic physical")
plt.plot(t_energy, energy_kinetic, label="--", label="Kinetic modal")
plt.plot(t_energy, energy_kinetic_modal, linestyle
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Kinetic energy of the system")
plt.title(
plt.show()
# Linear potential elastic energy
plt.figure()
plt.plot(="Linear potential elastic physical"
t_energy, energy_potential_elastic_linear, label
)
plt.plot(
t_energy,
energy_potential_elastic_linear_modal,="--",
linestyle="Linear potential elastic modal",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Linear potential elastic energy of the system")
plt.title(
plt.show()
# Linear potential stiffness energy
plt.figure()
plt.plot(
t_energy,
energy_potential_stiffness_linear,="Linear potential stiffness physical",
label
)
plt.plot(
t_energy,
energy_potential_stiffness_linear_modal,="--",
linestyle="Linear potential stiffness modal",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Linear potential stiffness energy of the system")
plt.title(
plt.show()
# Nonlinear potential elastic energy
plt.figure()
plt.plot(
t_energy,
energy_potential_elastic_nonlinear,="Nonlinear potential elastic physical",
label
)
plt.plot(
t_energy,
energy_potential_elastic_nonlinear_modal,="--",
linestyle="Nonlinear potential elastic modal",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Nonlinear potential elastic energy of the system")
plt.title(
plt.show()
# Total energy
plt.figure()="Total physical")
plt.plot(t_energy, energy_total, label="--", label="Total modal")
plt.plot(t_energy, energy_total_modal, linestyle
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Total energy of the system")
plt.title(
plt.show()
# Linear energy
plt.figure()="Linear physical")
plt.plot(t_energy, energy_linear, label="--", label="Linear modal")
plt.plot(t_energy, energy_linear_modal, linestyle
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Total linear energy of the system")
plt.title( plt.show()
# Plot the difference between the two energies
# Total energy
plt.figure()- energy_total_modal, label="Total")
plt.plot(t_energy, energy_total
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two total energies")
plt.title(
plt.show()
# Linear energy
plt.figure()- energy_linear_modal, label="Linear")
plt.plot(t_energy, energy_linear
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two linear energies")
plt.title(
plt.show()
# Kinetic energy
plt.figure()- energy_kinetic_modal, label="Kinetic")
plt.plot(t_energy, energy_kinetic
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two kinetic energies")
plt.title(
plt.show()
# Linear potential elastic energy
plt.figure()
plt.plot(
t_energy,- energy_potential_elastic_linear_modal,
energy_potential_elastic_linear ="Linear potential elastic",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two linear potential elastic energies")
plt.title(
plt.show()
# Linear potential stiffness energy
plt.figure()
plt.plot(
t_energy,- energy_potential_stiffness_linear_modal,
energy_potential_stiffness_linear ="Linear potential stiffness",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two linear potential stiffness energies")
plt.title(
plt.show()
# Nonlinear potential elastic energy
plt.figure()
plt.plot(
t_energy,- energy_potential_elastic_nonlinear_modal,
energy_potential_elastic_nonlinear ="Nonlinear potential elastic",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel("Difference between the two nonlinear potential elastic energies")
plt.title( plt.show()
# The enerygy of the system is not conserved, should be the same as the Kinetic energy initially added
# Plot the change in energy over time
assert np.allclose(energy_total[0], energy_kinetic[0])
assert np.allclose(energy_total_modal[0], energy_kinetic_modal[0])
plt.figure()# plt.plot(t_energy, energy_total/energy_total[0] -1, label="Total")
/ energy_total_modal[0] - 1, label="Total modal")
plt.plot(t_energy, energy_total_modal
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel(0.0, 0.02)
plt.xlim("Change in total energy of the system")
plt.title( plt.show()
# Plot the nonlinear energy as a fraction of the total energy
plt.figure()
plt.plot(
t_energy,/ energy_total,
energy_potential_elastic_nonlinear ="Nonlinear physical",
label
)
plt.plot(
t_energy,/ energy_total_modal,
energy_potential_elastic_nonlinear_modal ="--",
linestyle="Nonlinear modal",
label
)
plt.legend()"Time")
plt.xlabel("Energy")
plt.ylabel(0.0, 0.02)
plt.xlim("Nonlinear potential elastic energy as a fraction of the total energy")
plt.title( plt.show()
# Check that the energies calculated internally are the same as the ones calculated externally
assert np.allclose(energy_kinetic_modal[:300], e_k[:300])
assert np.allclose(energy_potential_elastic_linear_modal[:300], e_ple[:300])
assert np.allclose(energy_potential_stiffness_linear_modal[:300], e_pls[:300])
assert np.allclose(energy_potential_elastic_nonlinear_modal[:300], e_pne[:300])