Finite element simulation of Airy waves

Question

I am attempting to solve for waves on a water surface starting with a two dimensional solution. The equations are that the water must satisfy Laplace's equation everywhere with a time dependent boundary condition on the top surface. The equations are

$$\begin{align*} &\nabla^2\phi=0\\ &\left(\frac{\partial\phi}{\partial y}+\frac1{g}\frac{\partial^2\phi}{\partial t^2}\right)_{y=y_2}=0 \end{align*}$$

where $x$ is horizontal and $y$ vertical; also, $g$ is the acceleration due to gravity. The second equation applies at the surface. The first term in this equation is the vertical gradient so I think I can use a Neumann condition to enter this equation.

So I start by making a grid.

Needs["NDSolve`FEM`"];
x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
    MaxCellMeasure -> .002, "MaxBoundaryCellMeasure" -> 0.01];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

I now make some initial conditions. I have worked out an initial velocity for the top surface and solve to give the initial value of the potential function.

a = 0.40825787026798765;
b = 0.1689925306573793;
g = 9.81;

solIC = NDSolveValue[{
    Laplacian[ϕ0[x, y], {x, y}] ==
     NeumannValue[-b (E^-(x - 1)^2 - a), 0 <= x <= x2 && y == y2],
    DirichletCondition[ϕ0[x, y] == 0, x == 0 && y == 0]
    }, ϕ0, {x, y} ∈ mesh];

The velocity on the surface and throughout the water looks reasonable.

Plot[Evaluate[(D[solIC[x, y], y]) /. y -> y2], {x, 0, x2}, 
 AspectRatio -> 1/4]
ClearAll[f];
f[x_, y_] := Evaluate[Grad[solIC[x, y], {x, y}, "Cartesian"]];
StreamPlot[f[x, y], {x, 0, x2}, {y, 0, y2}, AspectRatio -> Automatic]

Now I attempt to set up the time dependent problem.

sol = NDSolveValue[{
   Laplacian[ϕ[x, y, t], {x, y}] ==
    NeumannValue[1/g Derivative[0, 0, 2][ϕ][x, y, t], 
     0 <= x <= x2 && y == y2],
   DirichletCondition[ϕ[x, y, t] == 0, x == 0 && y == 0],
   ϕ[x, y, 0] == solIC[x, y]
   }, ϕ,
  {x, y} ∈ mesh, {t, 0, 0.1}
  ]

This does not work and gives me the error

CoefficientArrays::poly: ϕ$1664+ϕ$1665-NeumannValue[0.101937 ϕ$1666,0<=x<=4&&y==1] 
is not a polynomial. >>

Now I am stuck. This does not help. Any suggestions?

Suggestions from comments

1. If I put in a time derivative in the equation this gives the same error.

2. If I take the time derivative out of the NeumannValue and put a time derivative into the equation then it solves but is meaningless. Does this suggest we can't have time derivatives in the boundary conditions?

Can you suggest how we might set up an alternative equation who's solution might go into the boundary condition?

sol = NDSolveValue[{
    Laplacian[ϕ[x, y, t], {x, y}] == 0.1 D[ϕ[x, y, t], t] +
      NeumannValue[ϕ[x, y, t], 0 <= x <= x2 && y == y2],
    DirichletCondition[ϕ[x, y, t] == 0, x == 0 && y == 0],
    ϕ[x, y, 0] == solIC[x, y]
    }, ϕ,
   {x, y} ∈ mesh, {t, 0, 0.1}
   ];
Plot3D[sol[x, y2, t], {x, 0, x2}, {t, 0, 0.1}]

So the problem could be the derivative in the NeumannValue.

Answer 1

Update: When I wrote this answer it was not clear to me that this question is about Airy wave theory which I do not know how to solve.

Old answer:

I am not sure if this is what you are looking for but here is a time dependent solution. We can take it from here if it's not what you need.

sol = NDSolveValue[{1/g D[ϕ[x, y, t], {t, 1}] - 
      Laplacian[ϕ[x, y, t], {x, y}] == 0, 
    DirichletCondition[ϕ[x, y, t] == 0, 
     x == 0 && y == 0], ϕ[x, y, 0] == solIC[x, y]}, ϕ, {x, 
     y} ∈ mesh, {t, 0, 0.1}];
Plot3D[Evaluate[{sol[x, y, 0], sol[x, y, 0.1]}], {x, y} ∈ 
  mesh, PlotLegends -> Automatic]

You can have derivatives of time in the NeumannValue - but that needs to be a wave type of equation:

sol = NDSolveValue[{1/g D[ϕ[x, y, t], {t, 2}] - 
      Laplacian[ϕ[x, y, t], {x, y}] == 
     NeumannValue[-Derivative[0, 0, 1][ϕ][x, y, t], 
      0 <= x <= x2 && y == y2], 
    DirichletCondition[ϕ[x, y, t] == 0, 
     x == 0 && y == 0], ϕ[x, y, 0] == solIC[x, y], 
    Derivative[0, 0, 1][ϕ][x, y, 0] == 0}, ϕ, {x, 
     y} ∈ mesh, {t, 0, 0.1}];

For this example, however, you'd need to set up an initial velocity that makes sense. Setting it to zero is not correct.

Here is a simplified example of how to set time derivatives within the NeumannValue. This is in essence an absorbing boundary condition (Note, that in this 1D case there is a small bug that will be fixed in a upcoming version, the sign in the NeumannValue is not correct. This has no effect on the 2D case you are looking at):

eqn = D[u[t, x], {t, 2}] == 
   D[u[t, x], {x, 2}] + 
    NeumannValue[
     If[$VersionNumber <= 10.4, 1, -1]*Derivative[1, 0][u][t, x], 
     x == 0 || x == 1];
u0[x_] := Evaluate[D[0.125 Erf[(x - 0.5)/0.125], x]];
ic = {u[0, x] == u0[x], Derivative[1, 0][u][0, x] == 0};
ufun = NDSolveValue[{eqn, ic}, u, {t, 0, 1}, {x, 0, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement"}}];
list = Table[
   Plot[ufun[t, x], {x, 0, 1}, PlotRange -> {-0.1, 1.3}], {t, 0, 1, 
    0.1}];
ListAnimate[list]