An accurate prediction of ship motions, slamming loads and structural responses has long been an important research interest in the field of ship hydrodynamics. With the tendency of modern ships towards large-scale and light-weighted, the significance of hydroelastic responses in ship hydrodynamic analysis has escalated, especially for the ships with large lengths. Ships are becoming more flexible, and the hull structural deformations can significantly affect the flow field and wave loads. In this situation, fully coupled Flexible Fluid Structure Interaction (FFSI) approach is needed in the numerical simulations of ship hydroelasticity analysis. In this paper, a two-way coupled CFD-MBD solver (Zhang et al., 2025) is applied to study the seakeeping performance and hydroelastic responses of a 20,000 TEU containership. The flow field is solved by RANS (Reynolds Averaged Navier-Stokes) equations with OpenFOAM, an open source CFD software. The structural deformation of hull girder is solved by MBDyn with multi-body dynamic method. The data transformation of fluid pressure and structural deformation in the FSI procedure is realized with the help of preCICE, an open source library for FSI communication (Xiao et al., 2024). The flow field mesh is updated by solving a Laplacian’s equation. Numerical predictions of ship motions and vertical bending moment (VBM) responses are validated against experimental measurements in regular head waves with varying wavelengths and ship speeds. The effects of wave height and ship speed on hydroelastic responses and slamming loads are systematically analyzed. The results show that VBM and slamming loads can increase significantly with changes in wave parameters, with more intense nonlinear effects observed under higher wave height conditions. In higher-speed cases, the peak value of VBM shifts from the midship to the bow region. These findings provide valuable data for the evaluation of ship hydroelasticity in ocean waves, particularly for ultra-large containerships.