Computation of wave-induced force

FORCE is the wave-driven stress, i.e. the force per unit surface driving the wave-driven current, expressed in N/m$^2$, is defined as the gradient of the radiation stresses:

  $$S_{xx} = \rho g \int \lfloor n \cos^2\theta + n - \frac{1}{2} \rfloor E d\sigma d\theta$$ (8.30)

  $$S_{xy} = S_{yx} = \rho g \int n \sin\theta \cos\theta E d\sigma d\theta$$ (8.31)

  $$S_{yy} = \rho g \int \lfloor n \sin^2\theta + n - \frac{1}{2} \rfloor E d\sigma d\theta$$ (8.32)

with $n$ the ratio of group velocity and phase velocity. The force is then

  $$F_x = -\frac{\partial S_{xx}}{\partial x} - \frac{\partial S_{xy}}{\partial y}$$ (8.33)

and

  $$F_y = -\frac{\partial S_{yx}}{\partial x} - \frac{\partial S_{yy}}{\partial y}$$ (8.34)

In order to compute the force in all internal vertices of the unstructured mesh, we consider a control volume (CV) as depicted in Figure 8.4.

Figure 8.4: Control volume (centroid dual) of vertex is shaded. Some notation is introduced.

This CV is called centroid dual and is constructed by joining the centroids neighbouring the vertex under consideration. The set of CVs must fill the whole computational domain and must also be non-overlapping. In the following we use the numbering from Figure 8.4. Let $\varphi$ be one of the radiation stresses $S_{xx}$, $S_{xy}$ and $S_{yy}$. The gradient of $\varphi$ is computed as follows

  $$\nabla \varphi \approx \frac{1}{A_{\rm CV}} \sum_e \varphi_e \vec{n}_e$$ (8.35)

where $A_{\rm CV}$ is the area of the CV and the summation runs over the associated edges $e$ of this CV. The values $\varphi_e$ at edges of the centroid dual are taken as averages, i.e. $(\varphi_0+\varphi_1)/2$, $(\varphi_1+\varphi_2)/2$, etc. Moreover, the value of the radiation stresses inside each triangle is simply the average of the radiation stresses in the associated vertices of the cell. Now, the derivatives of $\varphi$ inside CV are

  $$\frac{\partial \varphi}{\partial x} = \frac{1}{2A_{\rm CV}} \sum_{i=0}^{n-1} (\varphi_i+\varphi_{i+1})\,(y_{i+1}-y_i)$$ (8.36)

and

  $$\frac{\partial \varphi}{\partial y} = \frac{1}{2A_{\rm CV}} \sum_{i=0}^{n-1} (\varphi_i+\varphi_{i+1})\,(x_i-x_{i+1})$$ (8.37)

with $n$ the number of surrounding cells of the considered vertex and $\varphi_n = \varphi_0$, $x_n = x_0$ and $y_n = y_0$. The area of the CV is given by

  $$A_{\rm CV} = \frac{1}{2} \sum_{i=0}^{n-1} (x_i\,y_{i+1} - x_{i+1}\,y_i)$$ (8.38)