## 2 Numerical Implementation  
### Freak Wave Based on Superposition Model  
The superposition model, commonly employed as a freak-wave model, has been shown to be effective in generating freak waves experimentally and numerically. In the present study, a two-wave-train superposition model is adopted, which combines a random wave train and a convergent wave train to form the wave surface elevation. By adjusting the two trains of the waves' energy ratio and the phase of the linear cosine component waves of the convergent wave train, a freak wave is produced at a particular time and location [33]. The wave surface elevation of the two-wave-train superposition method can be expressed as:  
\[\zeta=\sum_{1i=1}^{N_{1}}a_{1i}\cos(o_{1i}x-w_{1i}t+\theta_{1i})+\sum_{2i=1}^{ N_{2}}a_{2i}\cos[o_{2i}(x-x_{c})-w_{2i}(t-t_{c})] \tag{1}\]  
In Equation (1), the first item creates the background random waves, while the second item generates the freak-wave peaks. \(N_{1}\) and \(N_{2}\) represent the respective quantities of component waves present in the two wave trains. \(o_{1i}\) and \(o_{2i}\) represent the wave numbers  
\[\mathbf{\sigma}^{s}=\mu^{s}\Big{[}\
abla\mathbf{z}+(\
abla\mathbf{z})^{T}\Big{]}+\lambda^{s}( \
abla\cdot\mathbf{z})\mathbf{I} \tag{8}\]  
where \(\mu^{s}\) and \(\lambda^{s}\) are the Lame Constants, \(\mathbf{I}\) is the identity tensor.  
The structural responses of the deck and breakwater structures are calculated by discretizing the momentum conservation equation for the structure, using FEM for spatial discretization and the \(Newmark-\beta\) method [38] for temporal discretization. The discrete form of the structural dynamic equation is:  
\[\mathbf{M}\cdot\ddot{\mathbf{z}}+\mathbf{C}\cdot\dot{\mathbf{z}}+\mathbf{K}\cdot\mathbf{z}=\mathbf{f} \tag{9}\]  
where \(\mathbf{M}\) is the mass, \(\mathbf{C}\) is the damping and \(\mathbf{K}\) is stiffness matrix matrices of the structure. \(\mathbf{f}\) is the force vector.  
Time can be discretized using the \(Newmark-\beta\) method [38], which involves the calculation of displacement, velocity, and acceleration within a specific time interval.  
\[\dot{\mathbf{z}}^{\mathbf{n+1}}=\dot{\mathbf{z}}^{\mathbf{n}}+\Big{[}(1-\beta)\ddot{\mathbf{z}}^{ \mathbf{n+1}}+\beta\cdot\ddot{\mathbf{z}}^{\mathbf{n+1}}\Big{]}\Delta t \tag{10}\]  
\[\mathbf{z}^{\mathbf{n+1}}=\mathbf{z}^{\mathbf{n}}+\dot{\mathbf{z}}^{\mathbf{n}}\Delta t+\Big{[}\Big{(} \frac{1}{2}-\gamma\Big{)}\ddot{\mathbf{z}}^{\mathbf{n}}+\gamma\cdot\ddot{\mathbf{z}}^{ \mathbf{n+1}}\Big{]}\Delta t^{2} \tag{11}\]  
where \(\Delta t\) is time step. \(\beta\) and \(\gamma\) are parameters that determine the accuracy and stability of the equation. Using Nastran as the structural solver in this paper, which is a widely used finite element analysis software, various dynamic loads  
Figure 2 illustrates the comparison between the heave and pitch motions of the container ship obtained using the present numerical method and experimental and numerical measurements conducted by [PERSON] et al. [43]. Due to the reflection of the beach and wave maker [44], as well as the shallow-water squatting phenomenon, there are some discrepancies between the numerical simulation results and the experimental results. However, good consistency was found when comparing them with the numerical results. The simulation results obtained from the present method were in strong agreement with the numerical simulations conducted by [PERSON] et al. [43]. This indicates that both the present solver and meshing have the ability to produce precise simulations of a ship's heave and pitch motions in waves.