?

Prediction of seakeeping in the early stage of conventional monohull vessels design using artificial neural network

2023-11-29 11:22RomeroTelloGutirrezRomeroServCms

P.Romero-Tello , ,J.E.Gutiérrez-Romero ,B.Serván-Cms

a Departamento de Física Aplicada y Tecnología Naval,Universidad Politécnica de Cartagena (UPCT),Cartagena,Murcia,Spain

b Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE),Barcelona,Spain

ABSTRACT Nowadays seakeeping is mostly analyzed by means of model testing or numerical models.Both require a significant amount of time and the exact hull geometry,and therefore seakeeping is not taken into account at the early stages of ship design.Hence the main objective of this work is the development of a seakeeping prediction tool to be used in the early stages of ship design.This tool must be fast,accurate,and not require the exact hull shape.To this end,an artificial intelligence (AI) algorithm has been developed.This algorithm is based on Artificial Neural Networks (ANNs)and only requires a number of ship coefficients of form.The methodology developed to obtain the predictive algorithm is presented as well as the database of ships used for training the ANN.The data were generated using a frequency domain seakeeping code based on the boundary element method (BEM).Also,the AI predictions are compared to the BEM results using both,ship hulls included and not included in the database.As a result of this work it has been obtained an AI tool for seakeeping prediction of conventional monohull vessels

1.Introduction.State of art

Thanks to the advantages offered by Artificial Intelligence (AI)and Machine Learning (ML) such as generalization ability,simplicity,adaptive learning,and others,ML are being used in many different fields of engineering.But not many ML techniques have been applied in marine engineering.Recently,Ao et al.[1] developed a Deep Learning algorithm capable to predict bare hull ship resistance.A database was generated using Free Form Deformation technique to get a large number of variations of the base containership KRISO.The targeted outputs were obtained using a potential flow code.Then they applied an optimization processes to the ANNs’ hyperparameters to obtain the best ANN.The resulting algorithm achieved promising results with an acceptable error when compared to the solver used to obtain the targeted output.Zhou et al.[2] optimized the hyperparameters using different ML techniques,such as,Artificial Neural Networks (ANNs),Support Vector Regression (SVR),and Random Forest (RF),in order to develop an algorithm capable of predicting fuel consumption taking into account weather and ship conditions.They concluded that ANNs algorithms can provide a robust and accurate estimation of fuel consumption.

Other as Cepowski [3,4] developed a set of ANNs capable of predicting added resistance in waves,sway accelerations,and roll angles for Car Carriers.The author used these ANNs to optimize car carrier ships.The main advantage of this set of ANNs is the ability to make predictions without needing the exact ship’s hull geometry.They only require the length to beam (L/B) and beam to draft (B/T) ratios,and hull shape coefficients as inputs.However,these trained algorithms are only capable of make predictions for Car Carriers within a limited range of principal dimensions.

AI has also been applied to seakeeping computing.Sayli et al.[5] used Non-linear meta-models [6,7] to predict the heave and the pitch Response Amplitude Operators for fishing vessels in head waves.It was used 13 different hull shapes,combining V and U sections.And the draft was modified to obtain 39 different ships.Then,using strip theory,seakeeping loads were computed in order to get the target outputs for training the AI algorithm.They used traditional naval architecture coefficients as input for ANNs,as Cepowski in [3,4] .And similarly,to previous works,ML techniques were used but they were tailored for specific types of ship.Alarrcin et al.[8] developed an ANN algorithm for the control algorithm of a container ship dynamic roll stabilizer.They used supervised learning with the most used network in regression problems,the Multi-Layer Perceptron (MLP) [9,10],combined with backpropagation.This ANN model consists of one input layer,several hidden layers,and one output layer.The main achievement is an ANN algorithm for dynamic roll control capable of significantly reducing roll amplitudes.

Ekinci et al.[11] demonstrated that AI algorithms can be applied successfully to estimate ship design parameters such as the beam,waterline length,deadweight,and draft to improve the ship definition in the early stages of design.They analyzed different AI methods and applied ANNs to optimize a chemical tanker.Abramowski [12] analyzed the effective power and preliminary data of conventional cargo ships depending on the ship velocity,vertical position of center of gravity,displacement,and cost.This author combined ANNs with other techniques such as genetic and simulated annealing algorithms to optimize the ship design.ANNs algorithms were trained to include a wide range of ships.Abramowski concluded that the combination of those methodologies opens wide possibilities when compared to traditional approaches used in the early stages of the ship design.

Yu and Wang [13] investigated the minimization of the ship wave resistance in calm waters using the Multi-Layer Perceptron(MLP).The original database used to train the MLP was formed by four well know hull forms: S60 with block coefficient Cb=0.6 and Cb=0.7);S175 series;and Wigley hulls.They applied transformation techniques on those four parent ships to increase the database to thousand different geometries.The algorithms inputs were parameters obtained from a principal component analysis,and the target bare hull resistance was obtained by computational analysis.

More recently Taghva et al.[14] predicted heave,roll,and pitch RAOs and added wave resistance of container ship based on the S175 hull using ANNs.As previous works,traditional naval architecture parameters were used to train the ANNs.The strip theory method was used to compute the seakeeping and added resistance in wave results for training the ANNs.

Other uses of AI are aiming at basic conceptual ship design.Gurgen et al.[15] trained ANN algorithms to predict principal dimensions for the conceptual design of chemical tankers.In this case,the database consists of 100 tanker ships.And ANNs were trained with basic ships parameters,such as freeboard,deadweight,length overall,and others.After training the ANNs,they obtained high correlation coefficients,and showing that AI can be used to determine initial ship’s particulars with high accuracy levels when compared to common regression methods.

ANNs can be sensitive to patterns between the input and the targets,difficult to reproduce with simple regression methods.However,the main drawback of ANNs is that they required large databases.For limited databases their advantages vanish and the risk of overfitting [9] increases,impeding generalization

Recently águila et al.[16] trained Recurrent Neural Networks(RNN) to predict seakeeping of the DTMB combatant ship in the time-domain.They compared Gated Recurrent Units (GRUs) and Long Short-Term Memory (LSTM) and concluded that LSTM are the best option for this specific analysis.GRUs and LSTM are types of RNN that address the ability of memory from previous points in time to help the network for future predictions.In this work,the authors used an original approach in which continuous functionals were approximated by LSTM.

Cepowski [17] continued previous works [3,4],and trained ANNS to predict added wave resistance using basic ship design parameters such as the length,beam,draft,Froude number and wavelength.The author showed the ability of AI to predict complex quantities with few basic data.While only few types of ships (passenger and ferry ships) can be predicted in previous works [3,4],in [17] Cepowski extended his algorithms to a few more typologies of ships.The dataset to train ANNs consists of computational results from 14 types of ships.Added resistance in waves coefficients were predicted by MLP,and compared with experimental results achieving a high accuracy.

Based on this literature review,several conclusions can be made:

· There are vast experimental seakeeping results,but they are scattered,limited,and non-homogenized to use in massive AI training.Then the large datasets of results required for training AI algorithms are obtained by numerical computations.

· Previous works have shown that,for many marine engineering applications,traditional form coefficient and principal dimensions are enough to define the hull geometry and feed ANNs.

· In previous works,the authors usually focused on obtaining ANN algorithms for very specific types of ships with a very limited training database.This results in algorithms that cannot be used for a wide range of ships,and might lead to overfitting.

2.Aim of this work

The main objective of this work is the development of Artificial Neural Networks (ANNs) capable of predicting the seakeeping of conventional ships.We define a conventional ship as a monohull vessel operating in displacement conditions.These ANNs will predict the Froude Krylov and wave diffraction and radiation loads as those computed using a standard frequency-domain firstorder wave diffraction-radiation solver based on the Boundary Element Method (BEM) [18,19] .These loads correspond to the incident wave load,added mass,damping,and wave diffraction loads,and hereafter will be refereed as seakeeping loads.

The purpose is to obtain a seakeeping tool capable of providing results instantly with an accuracy similar to those provided by a frequency domain BEM code.This tool will predict seakeeping based on basic hull form coefficients allowing to include seakeeping performance within the early stages of ship design.

This work is presented as follows: first it is introduced the methodology that has been used to model the seakeeping problem;second it is explained how the ANN training dataset has been obtain;third it is described how the different ANN architectures are generated,detailing the analyses carried out to select the best ANNs;fourth the results obtained from the ANNs are compared to those obtained by BEM for a number of ship shapes;and fifth results are discussed and conclusions are provided.

3.Methodology

The methodology used to generate the AI architecture is based on seven steps (see Fig.1).These steps are as follows:

Step1:Parentshipsgeometriesdatabase: Gathering of a parent ships database as broad as possible,including different types of conventional monohull ships (bulk-carriers,containerships,fishing vessels,etc.).

Step2:Shipsdatabaseaugmentation[18] : The database is augmented by parametric variations of the beam-length and draft-beam ratios.

Step 3: BEM computations and dataset generation: The wave diffraction-radiation problem is solved for every ship using a frequency-domain first-order solver based on the Boundary Element Method (BEM) [18,19] .

To carry out the comparison of results with a proper order of magnitude seakeeping loads are set dimensionless.Table 1 provides the corresponding dimensionless formulation for each hydrodynamics coefficient,taking as a reference term the diagonal term of the ship’s mass matrix (displacement and approximate inertias using the shown expressions in Table 4).Then the dataset of results to train the ANNs is generated.

Step4:ANNsgeneration: large number of ANNs are trained using parametric variations of their hyperparameters: number of layers,neurons,activation function,and optimization algorithms [9,21] .

Step5:ANNtraining: The ANNs are trained using the hull form coefficients as the input dataset.The targeted outputsDataset is composed of: added mass and damping matrices,and excitation loads.

Table 1 Dimensionless seakeeping loads.

Table 2 Range of the dimensionless form coefficients for the parent ships.

Step6:ANNscompetition: A competitive process among the different ANN architectures is carried out to select the best one(the one minimizing the errors).

Step7:Verification: Finally,a verification process is carried out comparing the ANN results against BEM results for a number of ships not included in the parent ships database,and not used neither for the training nor for the validation.

4.Data mining and simulation

4.1.Generating the dataset

A dataset of seakeeping results for a large database of ships is required for training the ANN algorithms [20] .In this work,the database of ships is generated out of 50 different parent ship geometries,considering a wide range of ship types.Fig.2 show the dimensionless cross-sectional area curves,and Tables 2 and A1 provide the form coefficients and the parent ships type.

For each parent ship,400 hulls are generated using 20 parametric variations of length-to-beam and beam-to-draft ratios.The ranges for these ratios are 1 .5 ≤≤8 .0

4.2.Seakeeping simulation particulars

Fig. 1.Methodology employed to generate the ANN architecture.

Fig. 2.Dimensionless cross-sectional area for each parent geometry.

In this work,a seakeeping BEM code,based on an improved version of Nemoh [18,22],is used to solve the first-order wave diffraction-radiation problem.A BEM code has been chosen instead of a strip code for several reasons.Based on Neumann[23] most strip theories are not valid for low frequencies.The strip theory does not have three-dimensional effects representing interactions between sections nor any forward-speed effects on the free-surface condition [24],compared with BEM formulation.Some comparison and validation about BEM code used can be found in a literature review.For instance,Parisella [25] compared the results of Nemoh with those computed by WAMIT [26],obtaining a low error.Anderson validated Nemoh with different calculation techniques [27] .For validation purposes,we compared the results obtained by BEM code used to obtain data base with those resulting from Journeé’s experimentation [28] .Fig.5 shows the comparison of heave and pitch RAO of the Wigley vessel.

For each of the 20,0 0 0 hulls,the first-order wave diffractionradiation problem has been solved for 30 waves with wave lengths between 0 .05 Lwland 1 .5Lwl,and for 7 incident wave directions between 0 and 180 °.Froude-Krylov and wave diffraction loads haven been obtained for each wave frequency and direction.And added mass and damping matrices haven been obtained for each wave frequency.

Fig. 3.Left: Amidship (Ac),midship (Am),and waterplane (Awl) sections of a ship.Right: Dataset representation based on sectional area ratios.

Fig. 4.Extreme parametric variation of a ship geometry.

Fig. 5.Comparison between BEM code and results for Wigley III Journeé’s experimentation Fr=0.2.

Fig.6 resumes the workflow strategy used to simulate the whole database of ships.Parametric variations of the parent’s geometries have been carried out directly on the mesh generated for the base geometries.

For each parent geometry,a mesh sensitivity analysis has been carried out to ensure that convergence have been reached.Examples of the test carried out are presented here.Tables 3–6 shows a convergence analysis for three hull forms (the original and two extreme parametric variations B/D=1.0 and L/B=1.5),for four types of ships.Tables 3–6 provide the error differences com puted between consecutive meshes for the normalized the diagonal terms of added mass,and damping matrices,and diffraction vector loads in percentage.From all those computed errors,the maximum is selected and shown in Tables 3–6 .The meshes are noted by M capital letter,the refinement in number of panels is increased between meshes,so Mi+1– Mishows error difference between them.Most of cases,mesh 1 has up to 200 panel,while mesh 8 reaches up to 10,0 0 0 panels.Figs.7–10 show the meshes analyzed.The error committed is less than 2% in most cases,when number of panels increases.The mesh selected is the number five.

4.3.ANN inputs and targeted outputs

In the early stages of the design of a ship,the exact hull geometry is unknown,and this is needed to predict the seakeeping particulars using numerical simulations or model testing.However,main particulars of the ship geometry,as the length,beam,draft,and form coefficient can be bounded within a range based on similar ships.

Fig. 6.Workflow for the seakeeping numerical simulation of the dataset of ships.

Table 3 Mesh sensitivity analysis for container ship hull form.

Table 4 Mesh sensitivity analysis for cruise ship hull form.

Table 5 Mesh sensitivity analysis for tanker ship hull form.

Fig. 7.Caption with three different meshes of a container parent ship.

Fig. 8.Caption with three different meshes of a cruise parent ship.

Fig. 9.Caption with three different meshes of a tanker parent ship.

Table 6 Mesh sensitivity analysis for frigate ship hull form.

Fig. 10.Caption with three different meshes of a frigate parent ship.

The objective of this work is to obtain a fast and accurate AI algorithm capable of predicting the seakeeping particulars of a conventional monohull ship,so that could be used to predict the seakeeping performance in the early stage of design.And since the exact geometry is not known,the AI should be based on the ship ′s geometry main particulars.

ANN inputs identification is a key point to enable the ANN to predict how changes in the hull geometry modify the seakeeping outputs.We define the ANN seakeeping target outputs as: the added mass and damping matrices,and the excitation loads for all six degrees of freedom.And these targets (ti) are assumed to mainly depend on the main dimensionless form factors of the hull geometry [29,30] .

where Lwlis the waterline length,Bwlis the maximum waterline beam,D is the draft,Cbis the block coefficient,Cwlis the waterplane coefficient,Cmis the midship coefficient,Ccis the amidships coefficient,andXBandZBare the longitudinal and vertical position of the buoyancy center,respectively.

For each nonzero component of the added mass and damping matrices,an ANN is built.An ANN is also built for every component of the excitation loads and wave heading.And 70% of the database of ships is used for the training,while 15% is used for the validation and the other 15% is used for the test.

5.Predictive algorithm

In this section,the process carried out to select the best ANN is presented.The selected ANN will be the one that minimize the prediction error.The Multi-Layer Perceptron (MLP) [9,10]model has been used for the regression problem.Tensorflow 2.1.0[31] along with GPU Nvidia acceleration libraries [32] have been used for the ANN development.GPU architecture are capable of speeding up the training process,achieving training times up to 10 times faster when compared to CPU architectures [9] .

Table 7 Parametric variation of ANN architecture hyper-parameters.

In order to obtain an optimal ANN architecture,a large number of ANNs with different parameters and hyper-parameters are trained and verified.Then,the errors are compared to select the best one.Eq.(1) provides a formulation for the mean absolute error (MAE) used [33–36],wheretiis the value obtained by numerical simulation,piis the ANN predicted value,andnis the total number of parameters.And Table 7 shows the parametric variation carried out to generate the different ANN architectures.

Fig. 11.Graphic representation for the ANN architectures generated by parametric variations of hyper-parameters.

Fig.11 provides a graphical representation of the different ANNs generated.The weights are initialized using a Gaussian distribution using the function Glorot Uniform function [47] .And to guarantee the tests repeatability four different seeds are used for initialization of the random functions.

Fig.12 shows the competition process applied to select best ANN to predict the added mass in heave (a33) .Each ANN (individual point in Fig.12) is trained with a combination of the hyperparameters from 1 to 5 (see Table 7).The best combination comes from the minimization of the MAE.Pareto front shows the minimum MAEs in the first step of the competition process applied to different number of neurons and layers.Then,overfitting prevent techniques (dropout),batch normalization,and regularization are applied,in next two septs,to those points having minimum MAE.In second step,dropout techniques are applied in those winner architectures from the first step.Then,in step three,the combinations of hyperparameters with lowest MAE are tested with batch normalization and regularization techniques.

After applying the three training steps shown in Fig.12,the ANN with the lowest MAE is selected.It is worth mentioning that steps 2 and 3 do not always improve on the training performed in step 1.Fig.12 shows an example of this,where the best ANN had been trained in the first step of the process.

From the training process it can be concluded that one layer is not enough to reach fitted predictions.Increasing the number of layers from one to two,the curves fitting improves up to 25%,and the increasing from two to three the gain is up to 15% additional.It should be noted,that the greater number of layers the greater risk of overfitting.

As for the number of neurons per layer,using 30 neurons fitted predictions are achieved with MAE errors below 5%.And it is observed that the Adam optimizer provides better results than the RMSprop.Furthermore,the results improve when the two activation functions are included in the ANN.

6.Numerical assessment and computational performance

6.1.Numerical assessment

This work focuses on predicting almost instantly the seakeeping loads for a wide range of conventional monohull ships.This range includes crude carriers,container ships,supply vessels,and many others as reported in Section 6 .

To demonstrate the validity ANNs developed,five assessment cases are computed and compared to the results obtained by the BEM code used to generate the target outputs of the ANNs.This inter-code comparison is carried for five hull geometries different to those of the parent ships database.Fig.13 compares the dimensionless form factor for the parent ships database and the assessment cases.Table 8 provides the main hull form parameters as well as their typology and Fig.14 shows sections and 3D views for the assessment cases.

Appendices B -E show graphical comparisons of the results obtained for ANNs versus BEM.The dimensionless diagonal values of the added mass and damping matrices are compared for different wavelengths.In addition,results of excitation loads for 7 headings and RAOs for 3 headings (0,30,and 60 °) are also compared.

6.2.Computational performance

One of the objectives of this work is to develop an AI capable of predicting seakeeping wave loads so that this AI can be used in the early stages of the ship design.Then,this tool should be capable of making fast predictions that can be used iteratively within the ship design process.

One of the main advantages of ANNs is that,once they have been trained and calibrated,they can compute complex problems almost immediately.In the case of this work,it is expected that seakeeping wave loads for a conventional monohull ship can be obtained much faster than using a traditional BEM code.And equally important,without the need of the exact hull geometry and the corresponding computational mesh.

Fig 15 shows the computational time required for computing all the seakeeping loads for 60 0 0 hulls,including the computational time for reading the input file,carrying out the ANN computation,and writing all the outputs.For each of those cases,30 wave frequencies and 7 wave directions are considered.This experiment was carried out with a computer with the following particulars:

-CPU: AMD Ryzen 7 3700 × 8 core Processor.

-GPU: Nvidia GeForce RTX2060.

-RAM: 32 GB.

It is observed that the computational velocity is about 200 cases/second.Given this computing velocity,and optimization algorithm based on a systematic variation of the ship form parameters could be used in the early design of a ship.

7.Results discussion

The fitting of an ANNs curve to a BEM curve is measured by the Mean Normalized Relative Error (MNRE) defined as follows,

beingNREithe normalized error in a pointi,the dimensionless target,the dimensionless prediction,andnthe number of curve points.

Fig.B.1 shows dimensionless values of added mass diagonal terms.Fig.C.1 compares the dimensionless diagonal values of damping matrices for the assessment ships.It can be observed that the ANNs underestimate the peak values of the damping curves.

Tables 9 and 10 shows MNRE corresponding to the prediction of the diagonal terms of the added mass and damping matrices.The average value of MNRE for the added masses is around 3.56%while the average value for the damping is 3.43%.

Fig. 13.Dimensionless ship form parameters for the parent ships used for the trainings and the five assessment ships.

Table 8 Main dimensionless hull parameters for the five assessment ships.

Table 9 Normalized relative mean error for predicted dimensionless added masses.

Table 10 Normalized relative mean error for predicted dimensionless dampings.

Fig. 14.Hull shapes below waterline for the five tests.Left: Body plan;Right: 3D view.

Fig. 15.Computational time in seconds versus number of cases computed.

Fig.D.1 shows the dimensionless excitation loads (sum of the diffraction and the Froude-Krylov loads) for surge () and pitch().It is compared for seven headings and four wavelengths.The MNRE values for the excitation forces are shown in Tables F.1 to F.5 .The average value for predicter forces goes from 1.78 to 3.28%.The largest deviations are observed for pitch moments with wave headings closed to the ship’s transversal direction.

Figs.E.1–E.5 show the RAOs for heave,roll and pitch for 0,30,and 60 ° headings.Tables G.1–G.5 provide MNRE values for the heave,roll,and pitch RAO curves.The average difference between curves ranges from 1.34% to 3.38%.Let’s point out that RAO curves are not directly predicted by ANN,but obtained from postprocessing the ANN results along with ships particulars such as position of the gravity center,and the mass and inertia matrix.

In general,the prediction of dimensionless added masses and damping values is better than those for the excitation loads.As a first approach,in this work it was assumed that optimizing an ANN for the added masses in heave (a33) would predict well the other outputs.This was assumed to reduce the large number of combinations given to solve full ANN optimization problem,leading to train millions of ANNs.Also,it is observed that trained ANNs naturally remove the irregular frequencies typically obtained from BEM codes,improving the targeted outputs.

8.Conclusions

This work demonstrate that algorithms based on ANNs could be suitable to predict the seakeeping loads needed to include seakeeping in early design stages of conventional monohulls ships.A large database and widespread has been used to achieve an accurate ANN.Starting from a database of 50 parent ships (see Table A1),data augmentation technique was used to increase the number of hulls up to 20,0 0 0.This augmentation was based on systematic geometric variations of the parent ships.Then,a BEM solver was used to compute a large output for training the ANNs.

One of the main achievements of this work is to obtain an AI algorithm capable of predict the seakeeping loads acting on a ship with no need of the exact geometry of the hull.Instead basic hull form coefficients,used in the early design,are required.And the generated algorithm will be able to predict the seakeeping particulars of any conventional monohull ship whose shape coefficients are within the limits shown in Table 1 .

A competitive process has been proposed to select best combination of hyper-parameters.Then a set of ANNs capable of predict different com ponents of the seakeeping loads have been obtained.

In this work it has been demonstrated the potential of ANNs to fast computation of seakeeping loads achieving more than 200 ships per second.Moreover,the ANNs can naturally remove irregular output data computed by BEM solvers.

DeclarationofCompetingInterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

AppendixA

Table A1 Parent ships.

AppendixB

Fig. B.1. Added mass comparison between ANN and BEM for the five assessment ships.

Fig. B.1. Continued

AppendixC

Fig. C.1.Damping comparison between ANN and BEM for the five assessment ships.

AppendixD

Fig. D.1.Wave excitation loads comparison between ANN and BEM for the five assessment ships.

Fig. D.1.Continued

Fig. D.1.Continued

AppendixE

Fig. E.1.Heave and Pitch RAO comparisons between ANN and BEM for the container assessment ship.

Fig. E.3.Heave and Pitch RAO comparisons between ANN and BEM for the Trawler assessment ship.

Fig. E.5.Heave and Pitch RAO comparisons between ANN and BEM for the anchor handling assessment ship.

AppendixF

Table F.1 Assessment ship 1 (Anchor Handling): MNRE for excitation loads.

Table F.2 Assessment ship 2 (Landing craft): MNRE for excitation loads.

Table F.3 Assessment ship 3 (Trawler): MNRE for excitation loads.

Table F.4 Assessment ship 4 (Yatch): MNRE for excitation loads.

Table F.5 Assessment ship 5 (Container ship): MNRE for excitation loads.

AppendixG

Table G.1 Assessment ship 1 (Anchor Handling): MNRE for RAO.

Table G.2 Assessment ship 2 (Landing craft): MNRE for RAO.

Table G.3 Assessment ship 3 (Trawler): MNRE for RAO.

Table G.4 Assessment ship 4 (Yatch): MNRE for RAO.

Table G.5 Assessment ship 5 (Container ship): MNRE for RAO.

91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合