I'm trying to model some two-dimensional waves, and am unsure how to combine my boundary conditions with spectral methods. The PDE I'd like to explore resembles the equations for damped, driven bending-mode waves in a thin disk or beam:
$$\ddot w(\boldsymbol x,t) = -c \nabla^2 \nabla^2 w(\boldsymbol x,t) - d \dot w(\boldsymbol x,t) + u(\boldsymbol x, t),$$
where $w(\boldsymbol x,t)$ is the spatiotemporal response of the system (analogous to vertical displacement of a horizontally oriented thin plate), $\boldsymbol x = \{x_1,x_2\}$ are 2D spatial coordiantes, $u(\boldsymbol x, t) = u(\boldsymbol x) u(t)$ is a spatiotemporally separable external input, and $- d \dot w(\boldsymbol x,t)$ is a damping/dissipation term with unknown $d$.
I've settled on this equation for my model based on an empirical dispersion relationship seen in experimental data; There may be other equivalent models. For example, the Schrödinger equation has a similar dispersion relationship and I'd be equally happy with methods employing that if its any easier.
The boundary I'd like to use looks like this:
The actual region I want to simulate isn't circular (or defined by any simple function), this was just easier to draw.
This is what I've done so far:
- I define a square 2D grid with resolution sufficiently smaller than the wavelengths of interest.
- I mask out a sub-domain of this grid with the irregular shape I'd like to model.
- I construct a graph adjacency matrix for the pixels/cells in this grid.
- I calculate the eigenvectors and eigenvalues of this matrix, and keep the top 16 components as a low-rank basis on which to approximate my solutions.
- These 16 components represent the top 16 eigenmodes of the Laplacian operator $\nabla^2$ with Dirichlet boundary conditions set to zero.
From here, numerical solutions are relatively straightforward. I prefer to think of the functions $w$ and $u$ as vectors in a Hilbert space, and write the spatial and time dependances as subscript indecies:
$$\ddot {\boldsymbol w}_{\boldsymbol x,t} = - a {\nabla^2 \nabla^2} {{\boldsymbol w}_{\boldsymbol x,t}} -d \dot {\boldsymbol w}_{\boldsymbol x,t} + {\boldsymbol u}_{\boldsymbol x} {\boldsymbol u}_t. $$
Take the Fourier transform in time , and change to the eigenbasis of the Laplacian operator in space. Denote temporal angular frequency as $\nu$. Denote spatial eigenmode number as $k$, with associated eigenvector $\boldsymbol v_{k}$ and eigenvalue $-\lambda_k$. We'll denote the thusly transformed components as $\boldsymbol\xi_{k,\nu}$. Note that the Fourier transform of the first and second (time) derivatives are $i\nu$ and $-\nu^2$, respectively. In this basis, the PDE becomes:
$$-\nu^2 {\boldsymbol\xi}_{k,\nu} = -a \lambda_k^2 \boldsymbol\xi_{k,\nu} - i\nu d {\boldsymbol\xi}_{k,\nu} + \hat{\boldsymbol u}_k \tilde{\boldsymbol u}_\nu,$$
where $\hat{\boldsymbol u}_k$ is the projection of $\boldsymbol u_{\boldsymbol x}$ onto the $k^{\text{th}}$ spatial mode and $\tilde{\boldsymbol u}_\nu$ is the (temporal) Fourier transform of ${\boldsymbol u}_t$. Collect terms and solve for ${\boldsymbol\xi}_{k,\nu}$:
$${\boldsymbol\xi}_{k,\nu} = \frac {\hat{\boldsymbol u}_k \tilde{\boldsymbol u}_\nu} {a \lambda_k^2 + i\nu d -\nu^2}$$
The spatiotemporal solution can then be obtained by applying the inverse Fourier transform to ${\boldsymbol\xi}_{k,\nu}$ and projecting back onto our 16 spatial components. This approach is much faster than simulating and waiting for said simulation to reach steady-state, and lets me explore parameter variations efficiently.
The issue, however is: how can one apply the same approach using the mixed ("Robin"?) boundary conditions described earlier? I suspect these boundary conditions should change the eigendecomposition, but I'm not sure how to construct/describe an operator encoding them?
Many thanks in advance! Disclaimer: I don't normally do this sort of thing, so apologies if I've mis-used some of the terminology.
