Numerical integration of damped wave equation

Question

I am trying to solve the equation

$$ \frac{d^2u}{dt^2}-\frac{d^2}{dx^2}\left(c_s^2u+\nu\frac{du}{dt}\right)=0 $$

with initial conditions

$$u(x, 0)=0$$

$$\frac{du}{dt}|_{t=0}=0$$

and boundary conditions

$$\frac{du}{dx}|_{x=0,L}=\frac{\rho_0\omega_d v_a}{c^4_s+\nu^2\omega_d^2}[c_s^2\sin{(\omega_d t)}-\nu\omega_d\cos{(\omega_d t)}] $$

My code is

ClearAll[\[ScriptC], \[Nu], \[Rho], \[ScriptL], \[ScriptN], \[Kappa], \[Omega], \[ScriptCapitalV], A, u, v];

\[ScriptC] = 1489; (*m/s - speed of sound*)
\[Nu] = 1.004;(*m^2/s - kinematic viscosity*)
\[Rho] = 998.21; (*kg/m^3 - density*)
\[ScriptL] = 2*10^-2; (*m - piezo distance*)
\[ScriptN] = 8; (*harmonic order*)
\[Kappa] = 2*Pi*\[ScriptN]/\[ScriptL]; (*1/m - wave number*)
\[Omega] = \[ScriptC]*\[Kappa]; (*rad/s - circular frequency*)
\[ScriptCapitalV] = \[Omega]*10^-8; (*m/s - piezo velocity*)

A = (\[Rho]*\[Omega]*\[ScriptCapitalV]/(\[ScriptC]^4+\[Nu]^4*\[Omega]^2));

\[Tau] = 10^-9;
(*f[t_] := (1-Exp[-t/\[Tau]]);*)
b[t_] := A*((\[ScriptC]^2)*Sin[\[Omega]*t] - (\[Nu]*\[Omega])*Cos[\[Omega]*t]);
b1[t_]:= \[Omega]*A*((\[ScriptC]^2)*Cos[\[Omega]*t] + (\[Nu]*\[Omega])*Sin[\[Omega]*t]);

f[t_] := If[t < 10^-9, 0, 1];
pde = {D[v[t, x], t] - (\[ScriptC]^2)*D[u[t, x], x, x] - \[Nu]*D[v[t, x], x, x] == 0, 
       v[t, x] == D[u[t, x], t]};
ics = {u[0, x] == 0, v[0, x] == 0};
bcs = {(D[u[t, x], x] /. x -> 0) == A*((\[ScriptC]^2)*Sin[\[Omega]*t] - (\[Nu]*\[Omega])*Cos[\[Omega]*t]) f[t],
       (D[u[t, x], x] /. x -> \[ScriptL]) == A*((\[ScriptC]^2)*Sin[\[Omega]*t] - (\[Nu]*\[Omega])*Cos[\[Omega]*t]) f[t]};
bcs1 = {(D[v[t, x], x] /. x -> 0) == \[Omega]*A*((\[ScriptC]^2)*Cos[\[Omega]*t] + (\[Nu]*\[Omega])*Sin[\[Omega]*t]) f[t],
        (D[v[t, x], x] /. x -> \[ScriptL]) == \[Omega]*A*((\[ScriptC]^2)*Cos[\[Omega]*t] + (\[Nu]*\[Omega])*Sin[\[Omega]*t]) f[t]};

tmax=0.0001;

{U, V} = NDSolveValue[{pde, ics, bcs, bcs1}, {u, v}, {x, 0, \[ScriptL]}, {t, 0, tmax}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 4000, "MinPoints" -> 4000, "DifferenceOrder" -> 4}}];

{DensityPlot[U[t, x], {x, 0, \[ScriptL]}, {t, 0, tmax}, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> Automatic, PlotLabel -> "u"], 
 DensityPlot[V[t, x], {x, 0, \[ScriptL]}, {t, 0, tmax}, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> Automatic, PlotLabel -> "v"]}

Manipulate[Plot[U[t, x], {x, 0, \[ScriptL]}], {t, 0, tmax}]

where I had to multiply to the boundary conditions a sigmoid (or step) function, otherwise these are inconsistent with the initial conditions.

The numerical solution I get with this code is inconsistent with what I expect from the analytical solution: an harmonic oscillation. I don't understand where the problem is.

Note: This is the follow-up of my previous question Solving a damped wave equation. I am now trying to use values of parameters with a physical meaning.