• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Multi-objective design optimization of composite submerged cylindrical pressure hull for minimum buoyancy and maximum buckling load capacity

    2021-09-02 05:39:00MuhmmdImrnDongynShiLiliTongAhsnElhiHfizMuhmmdWqsMuqeemUddin
    Defence Technology 2021年4期

    Muhmmd Imrn , Dong-yn Shi , , Li-li Tong , Ahsn Elhi , Hfiz Muhmmd Wqs , Muqeem Uddin

    a College of Mechanical and Electr ical Eng ineering Harbin Engineeri ng Universi ty, Harbin 150001, China

    b College of Aerospace and Civil Engineeri ng Harbin En gineering University, H arbin 1 50001, China

    c Pakhtunkhwa Energy Development Organization (PEDO), Khybe r P akhtunkhwa, Pakist an

    ABSTRACT This paper presents the design optimization of composite submersible cylindrical pressure hull subjected to 3 MPa hydrostatic pressure. The design optimization study is conducted for cross-ply layups [0s/90t/0u], [0s/90t/0u]s, [0s/90t]s and [90s/0t]s considering three uni-directional composites, i.e. Carbon/Epoxy, Glass/Epoxy, and Boron/Epoxy. The optimization study is performed by coupling a Multi-Objective Genetic Algorithm (MOGA) and Analytical Analysis. Minimizing the buoyancy factor and maximizing the buckling load factor are considered as the objectives of the optimization study. The objectives of the optimization are achieved under constraints on the Tsai-Wu, Tsai-Hill and Maximum Stress composite failure criteria and on buckling load factor. To verify the optimization approach, optimization of one particular layup configuration is also conducted in ANSYS with the same objectives and constraints.

    1. Introduction

    Fiber reinforced plastic (FRP) composites have very useful properties such as increased specific rigidity, increased elastic modulus, better fatigue properties, low weight, better corrosion resistance, and low magnetic and acoustic signatures [1]. Due to the aforementioned useful properties, FRP composites are used in the fabrication of high-speed vessels, mine countermeasure platforms, submerged pressure hulls and in numerous other parts of big naval ships. Also, as the strength to weight ratio of FRP is higher than that of metals and alloys, the composite pressure hulls shall have low weight for a specified working depth than other hulls constructed with steel and alloys [2]. The design of composite structures is very sensitive to fiber angle, the number of ply’s layers and the material system. To completely utilize the benefits of using composite materials in the construction of composite structures, it is very important to carry out optimization study of composite structures before finalization of the design for manufacturing or for further investigations such as nonlinear buckling analysis, impact analysis, and underwater explosion analysis, etc. Many studies on the optimization of structures constructed with composites carried out using a genetic algorithm (GA) and numerical techniques can be found in the open literature. A detailed review of the optimization of different composite structures has been published in Ref [3]. Several studies have focused on the optimization of composite plates utilizing a GA and other optimization methods linked to numerical techniques [4–13]. A lot of research has also been done on the optimization of composite submersible pressure hulls employing a GA and other optimization methods, finite element analysis (FEA) and analytical analysis. A brief review of the literature on the optimization of pressure hull constructed with composites is put forth in the following paragraph. Messager et al. [14] combined a GA with the analytical formulation for shell buckling to carry out the optimization of Carbon/Epoxy, and Glass/Epoxy cylindrical pressure vessels. The objective of the optimization was to maximize the buckling pressure. Liang et al. [15] used a Hybrid Genetic Algorithm (HGA) for optimization of Poly Vinyl Chloride (PVC) sandwich composite submerged pressure vessel. They used materials’ failure and shell buckling as constraints. Three different composites i.e. Graphite/Epoxy, Boron/Epoxy, and Glass/Epoxy were utilized for modeling of the face skin. Pelletier and Vel [5] optimized a Graphite/Epoxy cylindrical shell using a GA. The maximization of axial and hoop rigidities and minimization of the mass were used as objectives. Tsai-Wu composite material failure criterion was used as the material failure constraint. Keshtegar et al. [8] carried out study for buckling load optimization of laminated composites plate subjected to compressive loads. They combined the adaptive Kriging with enhanced partial swarm optimization algorithm (IPSO) and developed a hybrid optimization technique known as A-Kriging-IPSO. Ebrahimi et al [16] performed buckling study of nanocomposite shells by using First Order Shear Deformation Theory (FOSDT). The analytical model for the nanocomposite shell was developed by employing the principle of virtual work. The resultant stability equations were solved with the help of Galerkin’s method. Topal [17] used the Modified Feasible Direction (MFD) method to obtain optimized designs of a Graphite/Epoxy cylindrical shell. Maximizing the fundamental natural frequency and buckling pressure combined into a single objective and the ply orientation angle were selected as the objective and design variable of the optimization respectively. Fathallah et al. [18–20] used the Subproblem Approximation (SPA) technique in ANSYS to optimize the stiffened elliptical composite submerged pressure hull for minimum buoyancy factor. Tsai-Wu and Maximum Stress composite failure criteria and buckling load factor were used as constraints in the optimization process. Helal et al. [21] used the SPA technique in ANSYS for conducting the multi-objective design optimization of the sandwich core stiffened composite elliptical pressure hull. The objectives were to minimize the buoyancy factor and maximize the deck area and buckling load factor under similar constraints as used in Ref [18–20]. Craven et al. [22] carried out a conceptual design analysis of submerged Carbon/Epoxy and Glass/Epoxy composite pressure hull using ABAQUS. Their study was intended to reduce the weights of the composite pressure hulls subjected to constraints on material failure and buckling instability. Li et al. [23] proposed a methodology for the optimization of composite pressure hull utilizing the Ellipsoidal Basis Function Neural Network instead of the FEA in the layout optimization tool. They used the Tsai-Wu composite material failure criterion for evaluating composite material failure. Shen et al. [24,25] used a GA coupled with FEA analysis in ANSYS for optimization of composite submerged pressure hulls for achieving maximum buckling pressure subjected to restrictions on material failure and geometric instability . Imran et al. [26] conducted an optimization study of composite submersible pressure hull making use of a GA and FEA for minimum buoyancy factor. They used the Tsai-Wu and Tsai-Hill material failure criteria and buckling load factor as constraints on the optimization. Wei et al. [27] performed optimization of symmetrical laminated cylindrical pressure hull to achieve maximum buckling pressure. First, they combined a GA with FEA under no constraints to determine optimum lamination for the pressure hull and determined its ratios of extensional and bending stiffnesses. In the second step, they coupled a stiffness coefficient-based method (SCBM) using CLT with a GA to determine the extensional and bending stiffness ratios close to those previously determined in the first part of the study. Honda et al. [28] used NSGA-II for carrying out multi-objective optimization of shapes of curvilinear fiber for composite plates. Two types of conflicting objective functions, i.e. (1) maximizing the fundamental frequency and minimizing the average curvature value of the composite plate and (2) minimizing the failure index for Tsai-Wu failure criterion and minimizing the average curvature value of the composite plate, were used in the multi-objective optimization. Almeida Jr. et al. [29] carried out optimization of composite cylinder subjected to internal pressure by using a GA coupled with FEA. To predict damage and failure of the composite laminates in progressive failure analysis, a proposed damage model [30] was implemented through User Material Subroutine (UMAT) in ABAQUS. The fitness function was obtained for constrained and unconstrained optimization. Almeida Jr. et al. [31] in another study, carried out design optimization of a topologically-optimized composite structure by using a GA. The optimization was carried out for maximizing the specific stiffness subjected to the manufacturing constraints of the fiber arrangements of the composite structure. A 22% increase in the specific stiffness of the composite structure was obtained during optimization as compared to the topologically-optimized structure. Rouhi et al. [32] carried out optimization of a variable stiffness laminated cylinder under bending loads at the ends. To reduce the computation time the FEA was replaced by a surrogate model. First, two objective functions (structure performance at the two opposite ends) were defined and then a unified objective function was defined by using different weight factors for each objective function. Blom et al. [33] performed optimization study of a variable stiffness composite cylinder to obtain maximum buckling pressure under bending load. The design procedure was executed in ABAQUS for carrying out FEA. To reduce the computation duration, the FEA was replaced by a surrogate model in the optimization of the composite cylinder subjected to material failure constraint. Khani et al. [34] employed a semi-analytical model linked to a multi-level optimization method for optimization of stiffened composite cylinders under bending moments. The objective of the optimization was to achieve maximum buckling capacity subjected to material failure constraints. The aim of the present study is to find out the decrease in the buoyancy factor and increase in the buckling strength factor of the composite pressure hull with the use of composite materials instead of steel. The optimized design can be used for carrying out further detailed investigations such as non-linear buckling analysis and underwater explosion analysis etc. In the present study, the analytical formulations developed for material failure and buckling instability are based on linear elastic analyses. Linear elastic analysis has also been considered in various previously conducted conceptual design and design optimization studies of the composite pressure hulls and composite shells [14,15,18–20, 22,23,25,35–37]. In the linear analysis, the material is assumed to be linear elastic and free of imperfection. However, in actual circumstances, the structure may not remain linear elastic and free of imperfection due to many reasons. Therefore, in the present study, a minimum factor of safety of 2 is selected for both material and instability constraints. Moreover, the present study is based on the analytical modeling and the inclusion of imperfection in the analysis requires separate procedures such as axisymmetric imperfection (ASI), single load perturbation imperfection (SPLI), single boundary perturbation imperfection (SBPI), midsurface imperfection (MSI), geometric dimple-shaped imperfection (GDI), worst multiple perturbation load imperfection (WMPI), and linear buckling mode-shaped imperfection (LBMI) [38,39]. Therefore, the incorporation of imperfection in the analytical model is not possible during the optimization simulation process presented in this paper. However, the optimized designs obtained in this study, can further be investigated for imperfection sensitivity by using the aforementioned methods in finite element softwares such as ABAQUS [40]. The approach developed in the present study greatly reduces the computational time required to carry out the design optimization as compared to the time required for performing the same optimization by using FEA in commercially available software. The present paper considers four types of composite layup formations, including [0s/90t/0u], [0s/90t/0u]s, [0s/90t]s and [90s/0t]s, for design optimization study of cylindrical submersible pressure hull. These layup formations are considered in this study because the analytical modeling developed in the optimization of the composite pressure hull is based on the cross-ply layup formations. The same layup formations have alsobeen used in earlier studies on the design optimization of composite pressure hulls [19,22,26]. The optimization process is performed utilizing a Multi-Objective Genetic Algorithm of ISIGHT and Analytical Analysis performed in MATLAB. The optimization study uses two objectives: one is the minimization of the buoyancy factor, and the other is the maximization of the buckling load factor. Similarly, Tsai-Wu, Tsai-Hill, and Maximum Stress criteria and buckling load factor are considered as constraints on the optimization.

    2. Materials and Methods

    2.1. Analytical formulation for cylindrical composite pressure hull

    The Classical Lamination Theory (CLT) and Classical Shell Theory (CST) are employed to analyze the mechanical characteristics of composite laminated cylindrical pressure hull under hydrostatic pressure, depicted in Fig. 1. According to CLT, the stress-strain relation for orthotropic material according to the coordinate system presented in Fig. 1 can be written as [41].

    Fig. 1. Geometry of the composite submersible pressure hull and coordinates of plies in the laminates.

    where k=1,2,3,···,N represents the ply number from the bottom of the laminate andrepresents the elements of the transformed reduced stiffness matrixand are calculated by Eq. (2).

    where C = Cosθ, S = Sosθ and [Qij] are the reduced stiffnesses and can be calculated as follows.

    As stated by the First Order Shear Deformation Theory (FOSDT), the displacement relationships of the cylindrical shell can be represented as follows [16,41–45].

    where u0, v0and w0are the displacements at the midplane (z= 0) in the x, θ and z directions respectively.

    The kinematic relations for cylindrical shell in terms of strains and curvatures can be written as [46]:

    According to Classical Shell Theory (CST) for cylindrical shell [47],

    Substituting the above relations for rotations of the middle surface in the kinematic relations, we get

    Qatu [42] performed a number of analyses to investigate the outcome ofon the stress resultants for different types of composite and isotropic shells. The results showed that neglecting the termfor the cross-ply laminated thin deep cylindrical shell gave a 1.9% error. Therefore, it was concluded that, for the thin laminated composite shell, the relations of the composite laminated plate could be used for thin laminated shells. Using this finding, the resultant forces and moments at the middle surface of the shell are obtained as given below.

    where h represents the overall thickness of the shell and zk, zk-1represent the coordinates of the kthply calculated from the middle surface.

    The laminate strain–displacement relationships in terms of strains of the middle surface and curvatures of the shell are given as follows:

    Putting Eq. (10) in Eq. (1) we get,

    Putting Eq. (11) in Eqs. (8) and (9) and combining the resultant equations, we get,

    where Aij, Bijand Dijare extensional, bending extensional and bending stiffnesses respectively and are calculated as follows:

    In order to find the matrix, Eq. (12) is written in the fully inverted form as follows [48].

    where

    and

    2.2. Composite failure criteria

    The Tsai-Hill and Tsai-Wu failure criteria are the most commonly used interactive failure criteria for assessing composite materials’ failure [49]. Several earlier studies on the design optimization of composite pressure hull used Tsai-Wu [14,15,18–20,22,23,25,36] and Tsai-Hill [22,37] failure criteria for composites materials. The Tsai-Hill failure criterion is founded on the Von Mises’ distortional energy yield criterion when used for anisotropic materials. Failure index of the Tsai-Hill failure criterion (FITH) is determined by using Eq. (15). Contrary to other failure criteria, Tsai-Hill failure criterion does not differentiate between tensile and compressive stresses which can lead to underestimation of the load [50].

    The factor of safety of Tsai-Hill failure criterion (FSTH) is calculated by using the following formula [51].

    Tsai-Hill failure criterion can be modified to incorporate tensile and compressive stresses.

    where σ11andσ22are the stresses along the fiber direction and along perpendicular direction respectively and τ12is the in-plane shear stress. Xtand Xcare the tensile and compressive strengths in the X direction, Ytand Ycare the tensile and compressive strengths in the Y direction and S is the shear strength of unidirectional lamina.

    Tsai-Wu [52] failure criterion is based on the total strain energy theory of Beltrami [50]. Failure index of Tsai-Wu criterion is given by Eq. (17). This theory distinguishes between tensile and compressive stresses in its equation.

    where

    The factor of safety of Tsai-Wu failure criterion is calculated as follows [51].

    where

    As per the Maximum stress criterion, failure occurs when the stress developed along the principal direction surpasses its respective strength. Failure index of Maximum Stress criterion is given by the following Eq. (19).

    Similarly, the factor of safety (FSMS) for Maximum Stress failure criteria is calculated as follows.

    In addition, the principal and shear stresses in the composite shell can be calculated from the following equation.

    When the cylindrical pressure hull is subjected to external hydrostatic pressure, the stresses and moments resultants on the cylindrical shell are given by the following relationships [15,53,54].

    By inserting these values in Eq.14, we get the value ofwhich is substituted in Eq. (11) to get the value of. Then f r om Eq. (21), we getand finally, the failure indices for each layer of the laminate are calculated from Eqs.(15), (17) and (19). Similarly, from Eqs. (16), (18) and (20), the safety factors for each layer of the laminate are calculated.

    2.3. Analytical formulation for buckling of the composite pressure hull

    The relationship for critical pressure for composite cylindrical pressure hull subjected to external hydrostatic pressure is obtained by using the stationary form of Hamilton principle. The stationary form of Hamilton principle is represented by Eq. (23) [16,46,55].

    The strain energy (U), in the absence of transvers shear deformation, can be calculated as follows

    The work done by the external pressure on cylindrical can be calculated with the help of Eq. (25).

    where V is the change in volume per unit length enclosed by the mid-surface of the shell and is given by Eq.(26) [56].

    Where η is equal to 0 and 1 for dead and live pressure respectively. Similarly, f can take the values of 0 and 1 for Donnel’s theory and Sander’s theory respectively. In the case of dead loading, the applied pressure is considered to be centrally directed i.e. it remains perpendicular to the mid-surface of the shell in the undeformed state (original state), and in the case of live pressure, the applied pressure remains perpendicular to the mid-surface of the shell in the deformed state [56,57]. Now, using Sander’s theory, Eq.(25) for closed cylindrical pressure hull and live pressure can be written as follows.

    Putting Eqs. (24) and (27) in Eq. (23), we get

    Putting the strain relationships in Eq. (28) we get,

    Simplifying Eq.(29), the following equation is obtained.

    Integrating by parts Eq. (30), we get the following equation.

    Integrating by parts the remaining terms in Eq. (31), we get the following equation.

    Rearranging Eq. (32), we get the final simplified equation as follows.

    Simplifying and using Gauss-Ostrogradsky theorem, the second part of the Eq. (33) can be written as integrals over the area’s boundary (Γ).

    The stability equations are obtained by equating the coefficients under the double integralto zero.

    These equations of motion represent stability of the cylindrical pressure hull from buckling state to post buckling state. The pre-buckling or equilibrium state solution of the equations of motion is given by the following relations [55].

    where terms with bar sign indicate equilibrium or pre-buckling state.

    To find out the solution of the equations of motion in the buckling state, the equilibrium or pre-buckling state is disturbed by adding small increments to all quantities representing the pre-buckling state.

    where the terms with subscript (*)represent quantities in the disturbed state.

    Putting the quantities representing the buckling state in the governi ng Eqs. (35-37) and negl ecti ng the termsand, which represent the strains in the fundamental states, as well as the products of two terms in the disturbed state, we get the stability equations for cylindrical pressure hull subjected to hydrostatic pressure as follows.

    The stress and momentum resultant relations are rewritten for buckling state as follows.

    Similarly, for the buckling state, the strains and curvatures changes are linearized and are listed as below.

    Putting stresses and moments resultants from Eq. (41) and the linearized strains and curvatures changes in the stability Eqs. (38-40), we get the following system of equations for cross ply laminates composite submerged pressure hull subjected to hydrostatic pressure.

    The shell buckling governing Eqs. (43-45) are solved employing the subsequent buckling displacements relations, which satisfy simply-supported boundary conditions at the two ends of the pressure hull [15].

    where A, B, C are the buckling amplitudes’ coefficients and, where m and n are the number of half waves in the longitudinal direction and number of complete waves in the circumferential direction [27].

    Putting the above solution in shell buckling governing Eqs. (43-45), the following system of equations is obtained in matrix form.

    Eq. (47) represents a general Eigenvalue problem and the value of buckling load factor can be easily calculated in MATLAB for different values of m and n. The elements of matrix [C] and [L] are listed as follows.

    2.4. Problem’s solution procedure

    A MATLAB program is written for calculation of composite failure criteria and Eigenvalues. The MATLAB is linked to the design optimization module in ISIGHT [58]. In ISIGHT parameters, constraints and objectives are defined. The multi-objectives genetic algorithm NSGA-II is used for optimization of the composite cylindrical pressure hull subjected to 3 MPa hydrostatic pressure. The material properties of all composite materials considered in this study are given in Table 1. The length (l) and radius(R) of the pressure hull are 10.127 m and 1.370 m respectively. The length of the pressure hull is chosen based on a similar volume of cylindrical pressure hull with hemispherical ends used in our previous study of composite pressure hull [26].

    Table 1 Material properties of Carbon/Epoxy, Glass/Epoxy, and Boron/Epoxy [19,62].

    To verify the analytical model and design optimization procedure, separate studies are conducted for verification of composite material failure criteria, buckling load factor and design optimization procedure. The results of composite material failure criteria obtained from the present study are compared with similar results of the previous study [59]. The buckling load factor calculated by the present analytical procedure is compared with the result obtained from ANSYS. For verification of the optimization procedure, optimization analysis of one chosen layup formation is also performed in ANSYS with the same objectives and constraints. The detail of the finite element analysis and design optimization study performed in ANSYS for the verification of the analytical modeling and optimization approach developed in the present paper, is given bellow.

    2.5. Composite layup generation

    In the ANSYS Composite PrePost (ACP) geometric model of the pressure hull is generated by using the revolved feature of the design modeling tool. The geometric model is meshed with shell element of approximate size 0.001 m. The composite layup for the pressure hull is then generated by using setup of ACP. A Fabric is created and the Orientation and Reference Direction of fiber are defined by using corresponding features. The composite plies for each layer are defined with help of Modeling Groups.

    2.6. Boundary and loading conditions

    To prevent rigid body motion, simply supported boundary conditions are applied to the both ends of the pressure hull. In ANSYS, simply supported boundary conditions can be applied to the shell model. These boundary conditions are similar to the boundary conditions used in the analytical modeling of the pressure performed in the present paper. A uniform pressure of 3 MPa is applied to the external surface of the pressure hull.

    2.7. Procedure for finite element model solution

    After completing the composite layup in the ACP, parameters are created for the number of composite ply layers. The data comprising the composite shell model is then transferred to the Static Structural analysis. The finite element model is solved and parameters are created for the factors of safety of composite failure criteria. The data from static structural analysis is then transferred to linear Buckling analysis. The modeled is solved and the parameter for the buckling load factor is created. The data form the Parameters Tab in the ANSYS project is then transferred to the Direct Optimization component in Design Explorer. In the Design Explorer input parameters, constraints and objectives are defined. A parametric design analysis of a reference HY100 steel hull is also carried out in ANSYS for calculation of the percent decrease in the buoyancy factor (DBF) of the pressure hull. The yield strength and density of the HY100 steel are 790 MPa and 7828 kg/m3respectively [60,61]. The shell element of ANSYS is used for modeling of the composite pressure hull. The design optimization procedure flowchart for MATLAB-ISIGHT and ANSYS design optimization is shown in Fig. 2.

    Fig. 2. Optimization procedure flow chart.

    3. Statement of design optimization

    The NSGA-II of ISIGHT is used for optimization of the composite cylindrical pressure hull. The main features of the NSGA-II used in this study are given below.

    Initial population:The size of the initial population and the number of generations are selected based on the number of input parameters. In this study the size of the initial population is equal to 10 times the number of input parameters and the number of generations is equal to 40, so the size of the total population remains 1600 in case of all layup formations studied in the present study.

    Crossover:During Crossover, two parent chromosomes are combined to produce new chromosomes. In this study, the crossover probability of 0.9 is used with the crossover distribution index of 10. The crossover probability of 0.9 has been used in an earlier study on the design optimization of composite structures and the crossover distribution index of 10 is normally sufficient [28,63,64].

    Mutation:The mutation al ters one or more gene val ues in a chromosome from its ori ginal value. When new gene values are added to the chromosomes, then the GA can generate the best outcome as compared to those generated earlier. In this study, a mutation i ndex of 20 is used in accordance with the recommendation of the reference [64].

    Convergence Criteria:In the MATLAB and ISIGHT optimization, the program carries out calculations and computes optimum results based on design feasibility, objective and penalty functions. The optimization tool stops calculations after producing 40 generations.

    3.1. Objectives of optimization

    The objective function F (x) of the optimization used in this study is composed of two objective functions. The two objectives are: to minimize the buoyancy factor and to maximize the buckling load factor of the pressure hull subjected to hydrostatic pressure. The buoyancy factor (BF) has been utilized as an optimization objective in earlier research [17–20,25]. Buoyancy factor (BF) of the pressure hull is defined by Eq. (48).

    The buckling load factor (λ) is defined as the ratio of critical pressure (P ) to the applied pressure (P0) as given in (Eq. (49)).

    The two objectives functions F1(x) and F2(x) are combined into a single objective function F(x) by using weight factors of 0.5 for each objective.

    where Ω and μ are weight factors for objectives functions F1(x) and F2(x) respectively.

    3.2. Variables of design optimization

    The number of ply layers of the composite laminates are employed as variables of the design optimization.

    3.3. Constraints

    The following three types of constraints on the optimization are considered.

    Material constraints:To prevent material failure, the factors of safety for Tsai-Hill (FSTH) and Tsai-Wu (FSTW) and the Maximum stress (FSMS) composite failure criteria must be larger than or equal to 1. Also, in the case of steel hull, Von-Mises stress (σv) must be smaller than the yield stress (σy) of the steel. The material constraints are denoted by Eqs. (51-54).

    Structure instability constraint:Buckling load factor (λ) is considered as structure instability constraint. It must be greater than or equal to 1.

    Side constraints:The lower and upper limits for the number of layers of the laminate are considered as side constraints.

    where N, NUand NLrepresent the number of layers and its upper and lower limits for composite laminate.

    Constraint on buoyancy factor (BF):The buoyancy factor (BF) of composite pressure hull can have a value between 0 and 1. In course of this study, it has been observed that when no constraint is used for the buoyancy factor (BF) and buckling load factor, the optimization tool always produced results with increased buckling load factor and buoyancy factor (BF) due to the non-dominated sorting. In the multi-objective optimization, it is impossible to fulfill one objective without sacrificing the other objective. High buoyancy factor (BF) means increased number of composite layers and hence more computational time. Therefore, a strategy needs to be developed to fulfill all the objectives and obtain optimum results in less time and hence at a low computational cost. One possible strategy may be to use an upper limit on buckling load factor based on the safety factor. Using an upper limit on buckling load factor, the algorithm will search for inputs which may result in a buckling load factor which lies within the limits. But due to the orthotropic nature of composite, it is possible that inputs which result in high buckling load factor should also result in low buoyancy factor (BF). The other possible strategy is to use an upper limit on buoyancy factor (BF) for each composite material. Using an upper limit on buoyancy factor (BF) will force the algorithm to search for inputs which may result in high buckling load factor and low buoyancy factor (BF). In an earlier study on optimization of similar sized composite cylindrical pressure hull [26], the buoyancy factor for Carbon/Epoxy and Glass/Epoxy has been obtained in the range from 0.19 to 0.27, while for Glass/Epoxy, the buoyancy factor has been obtained in the range from 0.30 to 0.39 for different layup formations. Therefore, in the present study an upper limit of 0.25 on the buoyancy factor (BF) for both Carbon/Epoxy and Boron/Epoxy is selected, while an upper limit of 0.42 on the Buoyancy factor (BF) is selected for Glass/Epoxy.

    4. Results and discussions

    The results of parametric design analysis of steel hull are shown in Table 2. The buoyancy factor (BF) for the steel hull is 0.5. The optimum thickness of the hull is calculated as 46 mm and the maximum Von-Mises stress (σv) is 101.9 MPa. The buckling load factor (λ) is calculated as 2.05. The results of composite material failure criteria obtained from the present study are compared with similar results of the previous study [59] and are listed in Table 3. The results of composite failure criteria from both studies are almost similar to each other.

    Table 2 Parametric analysis of cylindrical submersible steel hull.

    Table 3 Results of composite failure criteria.

    To further verify the analytical formulations for composite failure criteria and buckling instability, comparison is also carried out between MATLAB and ANSYS calculations. The results for buckling load factor and safety factors for composite failure criteria for layup [01/9075/027], obtained from MATLAB and ANSYS (Fig. 3) are given in Table 4. As can be seen from Table 4, results of both studies closely follow each other. Similarly, to verify the overall MATLAB-ISIGHT and ANSYS optimization studies, the design optimization of one Carbon/Epoxy layup [0s/90t/0u] is also carried out in ANSYS. The results of MATLAB-ISIGHT and ANSYS optimization studies are presented in Table 5. In the case of ANSYS optimization analysis, the optimum layup is [05/9075/011]. The buoyancy factor (BF) for this layup is 0.2. The values of FSTW, FSTH, FSMS, and λ are 9.6, 11.8 and 12.3 and 2.0 respectively.

    Table 4 MATLAB and ANSYS calculations for optimized Carbon/Epoxy layup [01/9075/027].

    Table 5 ANSYS and MATLAB-ISIGHT optimization.

    Optimization history charts of objectives, constraints, design parameters and objective function for Boron/Epoxy layup [90m/0n]s are shown in Figs. 4-11. In the optimization history plots, the green points represent the optimum points as determined by the optimization tool, while the blue points represent the other optimum points. Furthermore, the black points in the history plots represent points in the design space which are feasible but not optimum, while the redpoints represent the infeasible points i.e. for these points the design feasibility is low and they violate the constraints limits as well. Similarly, pareto chart of objectives for Boron/Epoxy layup [90m/0n]s is shown in Fig. 12. In Fig. 12, the green point represents the optimum value of the buoyancy factor (BF) and buckling load factor (λ) determined by the optimization tool, while the blue points also represent other optimum values of the buoyancy factor (BF) and buckling load factor (λ) i.e. pareto optimal set. Moreover, the red points in Fig. 15 represent the infeasible points. The optimization analysis for one selected layup formation through ANSYS is run on a computer unit with specifications as: Intel (R) Xeon (R) CPU E5-2620 V2 @ 2.10 GHz (24 CPU), 32 GB RAM and 2 TB HD. The required time for this optimization analysis is about 10 h. On the other hand, the optimization simulations performed using the analytical approach developed in the present study required very less amount of computational time. The average time for one complete optimization simulation is about 15 min on a core i5, 8GB RAM personal computer. So, using the analytical approach developed in the present study considerably reduces the computational time and hence the computational cost.

    Fig. 3. Contour plots for Eigenvalue buckling (a) and safety factors for Tsai-Wu (b) Tsai-Hill (c) and Maximum Stress failure criteria (d) obtained from ANSYS.

    Fig. 4. Optimization history plot for buoyancy factor.

    Fig. 5. Optimization history plot for buckling load factor.

    Fig. 6. Optimization history plot for Maximum Stress failure criterion.

    Fig. 7. Optimization history plot for Tsai-Hill failure criterion.

    Fig. 8. Optimization history plot for Tsai-Wu failure criterion.

    Fig. 9. Optimization history plot for number of layers (s).

    Fig. 10. Optimization history plot for number of layers (t).

    Fig. 11. Optimization history plot for objective function.

    Fig. 12. Pareto chart of buoyancy factor and buckling load factor.

    The results of the design optimization of all layup formations used in the present paper are given in Tables 6-9 and discussed in the following paragraphs.

    4.1 Layup [0s/90t/0u]

    In the case of all layup formations and materials considered in the present paper, optimization is performed for achieving minimum buoyancy factor (BF) and maximum buckling load factor (λ). During optimization, the number of layers for all layup formations varies from 1 to 75. The size of the initial population is selected to be 10 times the number of input parameters. The optimization technique produces results based on the constraints and stops calculations after 40 number of generations. All the characteristics of the GA are kept constant except the size of the initial population. In the present study, only cross-ply layup formations are considered because the analytical model is derived based on the cross-ply layups only. The optimum values of parameters for layup [0s/90t/0u] are presented in Table 6. The lowest buoyancy factor (BF) is calculated for Carbon/Epoxy and equals 0.24 with a decrease of 52% in buoyancy factor (DBF). The buckling load factor (λ) for Carbon/Epoxy equals 2.5 with an increase of 22% in buckling load factor (IBSF) For the same material, the factors of safety of Tsai-Wu failure criterion (FSTW), Tsai-Hill failure criterion (FSTH), and Maximum stress failure criterion (FSMS) are equal to 14.5, 10.6, and 11.2 respectively. For Glass/Epoxy, the buoyancy factor (BF) equals 0.44 with DBFof 12%. The buckling load factor (λ) for the same material equals 2.6 with IBSFof 26.8% and its values of Tsai-Wu failure criterion (FSTW), Tsai-Hill failure criterion (FSTH), and Maximum stress failure criterion (FSMS) are equal to 10, 8.9, and 9.2 respectively. Similarly, the buoyancy factor (BF) for Boron/Epoxy is also calculated as 0.24 with DBF=52% and its value of buckling load factor (λ) equals 3.3 with IBSFof 61%. The optimized values of FSTW, FSTH, and FSMSfor Boron/Epoxy are obtained as 19.4, 16.7, and 16.8 respectively, which are larger than the respective values of FSTW, FSTH, and FSMSfor Carbon/Epoxy and Glass/Epoxy. For almost similar buckling performance, the material failure performance of Boron/Epoxy is better than that of Carbon/Epoxy. The material failure performance of Carbon/Epoxy is comparatively better than that of Glass/Epoxy. However, the decrease in the buoyancy factor (DBF) of Glass/Epoxy is considerably small as compared to those of Carbon/Epoxy and Glass/Epoxy.

    Table 6 Layup [0s/90t/0u].

    4.2. Layup [0s/90t/0u]s

    The symmetrical layup formation [0s/90t/0u]s is selected to compare its performance with the unsymmetrical layup formation [0s/90t/0u]. Moreover, the symmetric layup formation provides better resistance against thermal twisting of the cylindrical shell after curing when it cools down. The optimum values for design parameters of symmetric layup [0s/90t/0u]s are given in Table 7. The lowest buoyancy factor (BF) is calculated for Carbon/Epoxy and equals 0.24 with a decrease of 52% in buoyancy factor (DBF). Similarly, the buckling load factor (λ) for Carbon/Epoxy equals 5.3 with an increase of 158.5% in buckling load factor (IBSF) and its values of FSTW, FSTH, FSMSare calculated as 5.5, 4.1, and 4.3 respectively. For Glass/Epoxy, the buoyancy factor (BF) equals 0.35 with DBFof 30%. The buckling load factor (λ) for the same material equals 2.6 with IBSFof 26.8% and its values of FSTW, FSTH, and FSMSare equal to 5.9, 5.3, and 5.5 respectively. The buoyancy factor (BF) for Boron/Epoxy is calculated as 0.25 with DBF= 50% and its value of buckling load factor (λ) equals 3.9 with IBSFequal to 90.2%. Similarly, for Boron/Epoxy, FSTW, FSTH, FSMSare obtained as 9.2, 7.7, and 7.8 respectively, which are larger than the respective values of FSTW, FSTH, and FSMSfor Carbon/Epoxy and Glass/Epoxy. Carbon/Epoxy, Glass/Epoxy and Boron/Epoxy demonstrate better buckling performance in symmetrical layup formation [0s/90t/0u]s than their performances in the unsymmetrical layup formation [0s/90t/0u]. However, the material failure performance of all the three composites is better in the unsymmetrical layup formation [0s/90t/0u] than their performance in the symmetrical layup formation [0s/90t/0u]s . The overall performance of all the three composites is higher in symmetric layup [0s/90t/0u]s as compared to layup [0s/90t/0u]. Therefore, it is recommended to use the symmetrical layup formation instead of unsymmetrical layup to obtain better performance and to avoid thermal twisting of cylindrical shell during curing as well.

    Table 7 Layup [0s/90t/0u]s.

    4.3. Layup [0s/90t]s

    Optimum layups and optimum values for design parameters of layup [0s/90t]s are listed in Table 8. The lowest buoyancy factor (BF) is calculated for Boron/Epoxy and equals 0.21 with decrease of 58% in buoyancy factor (DBF), while its buckling load factor (λ) is calculated to be 2.52 with an increase of 22.9% in buckling load factor (IBSF). For Boron/Epoxy, the values of FSTW, FSTH, FSMSare obtained as 10.3, 8.6, and 8.7 respectively. Similarly, for Glass/Epoxy, the buoyancy factor (B.F) equals 0.42 with DBF= 16%. The buckling load factor (λ) for the same material equals 4.1 with IBSF= 100%, and its values of FSTW, FSTH, and FSMSare calculated to be 8.8, 7.9, and 8.2 respectively. The buoyancy factor (B.F) for Carbon/Epoxy is calculated as 0.22 with DBF= 56 % and its value of buckling load factor (λ) equals 3.6 with IBSF= 75.6%. For Carbon/Epoxy, the values of FSTW, FSTH, and FSMSare computed as 6.0, 4.5, and 4.6 respectively. The above results show that the material failure performance and weight reduction in the case of Boron/Epoxy are larger than that of Glass/Epoxy and Carbon/Epoxy. However, the buckling stability performance of Boron/Epoxy is smaller than that of Carbon/Epoxy and Boron/Epoxy. Similarly, in the same layup formation, the material failure and buckling performance of Glass/Epoxy is higher than that of Carbon/Epoxy. But the weight of the composite pressure hull constructed with Glass/Epoxy will be substantially higher as it gives smaller decrease in buoyancy factor as compared to Carbon/Epoxy and Glass/Epoxy.

    4.4. Layup [90s/0t]s

    The optimum values of design parameters of layup [90s/0t]s are listed in Table 9. For this layup, the lowest buoyancy factor (BF) is computed for Carbon/Epoxy and equals 0.22 with a decrease of 56 % in buoyancy factor (DBF), while its value of buckling load factor (λ) is computed as 3.9 with increase of 90.2% in buckling load factor (IBSF). The values of FSTW, FSTH, and FSMSfor Carbon/Epoxy are computed as 9.3, 7.7, and 7.2 respectively. Similarly, for Glass/Epoxy, the buoyancy factor (BF) equals 0.42 with DBF= 16%, while its value of buckling load factor (λ) equals 4.3 with IBSF= 109.8%. The values of FSTW, FSTH, and FSMSfor Glass/Epoxy composites are calculated as 6.7, 6.5, and 6.4 respectively. The buoyancy factor (B.F) for Boron/Epoxy composite is calculated as 0.24 with decrease of 52% in buoyancy factor (DBF), while its value of buckling load factor (λ) is obtained as 3.6 with an increase of 75.6% in buckling load factor (IBSF). For Boron/Epoxy composite, the values of FSTW, FSTH, FSMSare obtained as 15.2, 13, and 13.2 respectively, which are larger than the respective values of FSTW, FSTH, and FSMSfor Carbon/Epoxy and Glass/Epoxy. For almost similar value of decrease in buoyancy factor (DBF), the buckling performance of all the three composite materials in layup [90s/0t]s is better than their performance in layup [0s/90t]s. For almost similar decrease in buoyancy factor (DBF) the overall performance of all the three composites is better in layup formation [90s/0t]s as compared to layup formation [0s/90t]s. These results show that, to obtain high factors of safety of material failure criteria and buckling load factor, the fibers in the bottom most layer of the composite layer should be oriented along the radial direction of cylindrical pressure hull.

    Table 8 Layup [0s/90t]s.

    Table 9 Layup [90s/0t]s.

    5. Conclusion

    This paper describes the optimization of the composite submersible cylindrical pressure hull considering a genetic algorithm (GA) and analytical formulation. The overall design optimization process is accomplished by coupling ISIGHT and MATLAB. The cross-ply layup formations [0s/90t/0u], [0s/90t/0u]s, [0s/90t]s and [90s/0t]s utilizing three uni-directional composites, Carbon/Epoxy, Glass/Epoxy, and Boron/Epoxy are employed for creation of composite layups of the pressure hull. Minimizing the buoyancy factor (B.F) and maximizing the buckling load factor (λ) are used as the objective functions of the optimization under 3 MPa hydrostatic pressure. The Tsai-Wu, Tsai-Hill and Maximum stress composite materials failure criteria and buckling load factor are utilized as constraints on the optimization. The optimum layups are compared based on the decrease in buoyancy factor (DBF) and increase in buckling load factor (IBSF) over equivalent steel hull. The main conclusions of the present study are:

    ? Analytical formulation based on Classical Lamination Theory (CLT) and Classical Shell Theory (CST), and design optimization tools can be efficiently exploited for design optimization of structures modeled with composites like the composite submersible pressure hull. The procedure adopted in this paper not only avoids cumbersome finite element analysis but also greatly reduces computational time.

    ? A maximum decrease of 58% in buoyancy factor (DBF) over the reference steel hull is calculated for Boron/Epoxy. Also, Boron/Epoxy gives higher factors of safety for material failure criteria in the case of all 4 layup formations.

    ? A maximum increase of 158.5% in buckling load factor (IBSF) over the reference steel hull is calculated for Carbon/Epoxy. Similarly, Carbon/Epoxy gives overall higher decrease in buoyancy factor (DBF) over the reference steel hull as compared to Boron/Epoxy and Glass/Epoxy.

    ? Furthermore, the overall performance of all the three composites in symmetric layup formations is better than their performance in asymmetric layup formations, since the former result in higher values of design parameters as compared to the latter.

    Declaration of interests

    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.

    Acknowledgement

    This work is supported by the National Natural Science Foundation of China research grant “Study on the characteristic motion and load of bubbles near a solid boundary in shear flows” (51679056) and Natural Science Foundation of Heilongjiang Province of China (E2016024).

    久久精品国产鲁丝片午夜精品| 亚洲av不卡在线观看| 亚洲美女视频黄频| 男人舔奶头视频| av天堂中文字幕网| 国产一区二区三区综合在线观看 | 美女被艹到高潮喷水动态| 高清在线视频一区二区三区| 免费看av在线观看网站| 99热全是精品| 三级毛片av免费| 国内精品宾馆在线| 亚洲欧美中文字幕日韩二区| 在现免费观看毛片| 欧美zozozo另类| 噜噜噜噜噜久久久久久91| 有码 亚洲区| 菩萨蛮人人尽说江南好唐韦庄| 3wmmmm亚洲av在线观看| 国产精品av视频在线免费观看| 国产三级在线视频| 91aial.com中文字幕在线观看| 国产真实伦视频高清在线观看| av在线观看视频网站免费| 久久久亚洲精品成人影院| 纵有疾风起免费观看全集完整版 | 国产 一区精品| 亚洲熟妇中文字幕五十中出| 国产精品不卡视频一区二区| 免费看不卡的av| 小蜜桃在线观看免费完整版高清| 青春草视频在线免费观看| 中文字幕久久专区| 99热这里只有是精品50| 久久久久久久亚洲中文字幕| 成人综合一区亚洲| 国产色婷婷99| av专区在线播放| 99久久精品一区二区三区| 国产激情偷乱视频一区二区| 免费人成在线观看视频色| 秋霞伦理黄片| 天堂av国产一区二区熟女人妻| 熟妇人妻不卡中文字幕| 欧美成人精品欧美一级黄| 久久久久久久大尺度免费视频| 国产精品一二三区在线看| 午夜精品在线福利| 亚洲国产色片| 在线免费观看不下载黄p国产| 午夜爱爱视频在线播放| 久久99热6这里只有精品| 水蜜桃什么品种好| 一级黄片播放器| 免费看日本二区| 又爽又黄无遮挡网站| 亚洲国产精品国产精品| 亚洲天堂国产精品一区在线| 国产精品麻豆人妻色哟哟久久 | 亚洲伊人久久精品综合| 亚州av有码| 亚洲精品乱久久久久久| 国产亚洲午夜精品一区二区久久 | 色吧在线观看| av又黄又爽大尺度在线免费看| 国产av不卡久久| 久久综合国产亚洲精品| 国产黄a三级三级三级人| 亚洲国产欧美在线一区| 久久人人爽人人爽人人片va| 久久久精品94久久精品| 一区二区三区四区激情视频| 男插女下体视频免费在线播放| 国产老妇伦熟女老妇高清| 亚洲一区高清亚洲精品| 亚洲久久久久久中文字幕| 久久精品国产亚洲网站| 成年av动漫网址| 久久久精品免费免费高清| 少妇熟女aⅴ在线视频| 亚洲成色77777| 日本-黄色视频高清免费观看| 日韩欧美 国产精品| 国产精品一区二区性色av| 日韩欧美精品v在线| 搡女人真爽免费视频火全软件| 日韩大片免费观看网站| xxx大片免费视频| 亚洲精品国产成人久久av| 搡老妇女老女人老熟妇| 毛片女人毛片| 成人性生交大片免费视频hd| 欧美成人精品欧美一级黄| 别揉我奶头 嗯啊视频| 夜夜爽夜夜爽视频| 国产精品不卡视频一区二区| 看免费成人av毛片| 老司机影院毛片| 在线观看av片永久免费下载| av国产免费在线观看| 国产精品久久久久久精品电影| 一边亲一边摸免费视频| 国产亚洲午夜精品一区二区久久 | 色哟哟·www| 久久综合国产亚洲精品| 国产精品久久久久久av不卡| 国产片特级美女逼逼视频| a级毛色黄片| 亚洲av在线观看美女高潮| 联通29元200g的流量卡| 在线播放无遮挡| 国产精品一区二区三区四区免费观看| 国产一级毛片七仙女欲春2| 免费看a级黄色片| 国内揄拍国产精品人妻在线| av免费在线看不卡| 男女边吃奶边做爰视频| 国产黄a三级三级三级人| 日韩成人av中文字幕在线观看| 欧美xxxx黑人xx丫x性爽| 日韩亚洲欧美综合| 免费看av在线观看网站| 又爽又黄a免费视频| 少妇的逼水好多| .国产精品久久| 日韩av不卡免费在线播放| 少妇被粗大猛烈的视频| 激情五月婷婷亚洲| 国产在视频线精品| 亚洲精品456在线播放app| 日日啪夜夜爽| 久久久精品免费免费高清| av网站免费在线观看视频 | 精品一区在线观看国产| 99久久精品国产国产毛片| 啦啦啦啦在线视频资源| 老师上课跳d突然被开到最大视频| 中文资源天堂在线| 国产一级毛片七仙女欲春2| 日产精品乱码卡一卡2卡三| 国产伦精品一区二区三区视频9| 亚洲18禁久久av| 免费观看性生交大片5| 男人和女人高潮做爰伦理| 蜜桃亚洲精品一区二区三区| 99热这里只有是精品50| 伊人久久精品亚洲午夜| 国产女主播在线喷水免费视频网站 | av专区在线播放| 久久久久久久久久久免费av| 熟女电影av网| 国产av在哪里看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 男女啪啪激烈高潮av片| 26uuu在线亚洲综合色| 男女边吃奶边做爰视频| 日韩av在线大香蕉| 日韩 亚洲 欧美在线| 别揉我奶头 嗯啊视频| 亚洲精华国产精华液的使用体验| 久久久久久久久久成人| 国产高清不卡午夜福利| 婷婷色麻豆天堂久久| 亚洲综合色惰| 一级毛片电影观看| 婷婷色综合www| 如何舔出高潮| 国产成人午夜福利电影在线观看| 免费无遮挡裸体视频| freevideosex欧美| 国产一区亚洲一区在线观看| 超碰av人人做人人爽久久| 亚洲精品第二区| 高清毛片免费看| 日日摸夜夜添夜夜添av毛片| 国产午夜福利久久久久久| 又大又黄又爽视频免费| 国产亚洲最大av| 在线观看免费高清a一片| 免费在线观看成人毛片| 亚洲天堂国产精品一区在线| 成人午夜精彩视频在线观看| 亚洲成人av在线免费| 日本色播在线视频| 一级毛片 在线播放| 69人妻影院| 国产精品国产三级国产专区5o| 精品久久久噜噜| h日本视频在线播放| 国产一区二区三区av在线| 午夜激情久久久久久久| 好男人视频免费观看在线| 亚洲精品中文字幕在线视频 | 国产精品久久久久久久久免| 亚洲丝袜综合中文字幕| 亚洲成色77777| 日韩伦理黄色片| av在线天堂中文字幕| 国产淫片久久久久久久久| 国产乱人偷精品视频| 大香蕉久久网| 久久久亚洲精品成人影院| 五月玫瑰六月丁香| 亚洲精品影视一区二区三区av| 免费看av在线观看网站| 免费av观看视频| 男女边摸边吃奶| 国产 亚洲一区二区三区 | 99九九线精品视频在线观看视频| 婷婷色综合www| 亚洲av不卡在线观看| 亚洲av中文字字幕乱码综合| 午夜福利在线观看免费完整高清在| 亚洲精品成人久久久久久| 97在线视频观看| 欧美xxxx黑人xx丫x性爽| 晚上一个人看的免费电影| 国产伦一二天堂av在线观看| 在线观看av片永久免费下载| 亚洲国产欧美在线一区| av国产久精品久网站免费入址| 一级av片app| 免费观看精品视频网站| 亚洲最大成人中文| 婷婷色麻豆天堂久久| 爱豆传媒免费全集在线观看| 嫩草影院入口| 免费看日本二区| 免费av观看视频| 只有这里有精品99| 亚洲内射少妇av| 午夜福利视频1000在线观看| 精品99又大又爽又粗少妇毛片| 亚洲,欧美,日韩| 婷婷色麻豆天堂久久| 久久久久久久午夜电影| 2021少妇久久久久久久久久久| 国产精品蜜桃在线观看| 国产午夜精品论理片| 午夜日本视频在线| 久久久国产一区二区| 男人和女人高潮做爰伦理| 国产中年淑女户外野战色| 婷婷色av中文字幕| 色5月婷婷丁香| 男女边摸边吃奶| 色哟哟·www| 国产美女午夜福利| 日本一本二区三区精品| 亚洲成色77777| 在线播放无遮挡| 亚洲第一区二区三区不卡| 亚洲欧洲国产日韩| 久久6这里有精品| 韩国高清视频一区二区三区| 国产伦精品一区二区三区四那| 在线观看av片永久免费下载| 一级毛片电影观看| 精品人妻熟女av久视频| 亚洲天堂国产精品一区在线| 国产精品一区二区三区四区免费观看| 亚洲精品久久久久久婷婷小说| 中国国产av一级| 久久精品久久久久久久性| 日韩成人伦理影院| freevideosex欧美| 国产高清国产精品国产三级 | 日本黄大片高清| 最近最新中文字幕免费大全7| 一本久久精品| 男女那种视频在线观看| 老司机影院毛片| 美女国产视频在线观看| av国产久精品久网站免费入址| 九九爱精品视频在线观看| 国产成人freesex在线| 国产极品天堂在线| 成人特级av手机在线观看| 国产乱来视频区| 人妻一区二区av| 免费观看a级毛片全部| 校园人妻丝袜中文字幕| 国产精品一区www在线观看| 联通29元200g的流量卡| 午夜日本视频在线| 精品久久国产蜜桃| 日韩一区二区三区影片| 国产伦一二天堂av在线观看| 国产午夜精品久久久久久一区二区三区| 淫秽高清视频在线观看| 麻豆乱淫一区二区| 色5月婷婷丁香| 日日啪夜夜撸| 嫩草影院精品99| 亚洲精品成人av观看孕妇| 久久久欧美国产精品| av在线观看视频网站免费| 啦啦啦中文免费视频观看日本| 精品99又大又爽又粗少妇毛片| 波野结衣二区三区在线| 国内揄拍国产精品人妻在线| 久久韩国三级中文字幕| 亚洲国产成人一精品久久久| 看十八女毛片水多多多| 菩萨蛮人人尽说江南好唐韦庄| 最后的刺客免费高清国语| 日本一本二区三区精品| ponron亚洲| 狂野欧美激情性xxxx在线观看| 在线免费十八禁| 免费人成在线观看视频色| 国产精品国产三级专区第一集| 嫩草影院精品99| 婷婷色av中文字幕| 国产精品一区二区在线观看99 | 最近中文字幕高清免费大全6| 99热这里只有精品一区| 人人妻人人澡人人爽人人夜夜 | 日韩av在线大香蕉| 天堂影院成人在线观看| 一级毛片黄色毛片免费观看视频| 久久久久久久久久黄片| 在现免费观看毛片| 波多野结衣巨乳人妻| 国产亚洲5aaaaa淫片| 97在线视频观看| 久久草成人影院| 亚洲成色77777| 午夜激情久久久久久久| 99热全是精品| 黄色一级大片看看| 国产精品一及| 国产美女午夜福利| 人妻一区二区av| 99久久精品热视频| 国产在视频线精品| 又黄又爽又刺激的免费视频.| 精品人妻偷拍中文字幕| 成人特级av手机在线观看| 国产一级毛片在线| 一个人看的www免费观看视频| 久久久国产一区二区| 国产乱来视频区| 一个人看的www免费观看视频| 丰满人妻一区二区三区视频av| 午夜福利视频1000在线观看| 夜夜爽夜夜爽视频| 亚洲真实伦在线观看| 国产黄色小视频在线观看| 天堂影院成人在线观看| 午夜福利高清视频| 亚洲av一区综合| 色尼玛亚洲综合影院| 亚洲第一区二区三区不卡| av一本久久久久| 女的被弄到高潮叫床怎么办| 女人久久www免费人成看片| 免费av毛片视频| 久久久久精品性色| 高清视频免费观看一区二区 | 久久精品夜色国产| 少妇高潮的动态图| 一级片'在线观看视频| 韩国av在线不卡| 国产精品.久久久| 美女黄网站色视频| 午夜激情福利司机影院| 国内精品一区二区在线观看| 亚洲欧洲国产日韩| 直男gayav资源| 欧美日韩综合久久久久久| 男女边吃奶边做爰视频| 亚洲国产成人一精品久久久| 免费观看在线日韩| 色综合站精品国产| 蜜桃久久精品国产亚洲av| av.在线天堂| 99九九线精品视频在线观看视频| 极品少妇高潮喷水抽搐| 身体一侧抽搐| 视频中文字幕在线观看| 狂野欧美白嫩少妇大欣赏| 久久久久久久久久久免费av| 国产欧美另类精品又又久久亚洲欧美| 少妇人妻精品综合一区二区| 少妇高潮的动态图| 爱豆传媒免费全集在线观看| 天堂俺去俺来也www色官网 | 免费播放大片免费观看视频在线观看| 亚洲四区av| 在线播放无遮挡| av免费在线看不卡| 国产69精品久久久久777片| 欧美不卡视频在线免费观看| 天堂√8在线中文| 女的被弄到高潮叫床怎么办| 天堂中文最新版在线下载 | 精品久久久久久电影网| 国产麻豆成人av免费视频| 熟妇人妻久久中文字幕3abv| 在线免费观看的www视频| 青春草国产在线视频| 99视频精品全部免费 在线| 国产精品爽爽va在线观看网站| 天堂影院成人在线观看| 亚洲真实伦在线观看| 欧美成人一区二区免费高清观看| 国产精品麻豆人妻色哟哟久久 | 秋霞在线观看毛片| 国产白丝娇喘喷水9色精品| 色哟哟·www| 欧美高清性xxxxhd video| 日韩一本色道免费dvd| 最近中文字幕高清免费大全6| 在线观看人妻少妇| 蜜臀久久99精品久久宅男| 国语对白做爰xxxⅹ性视频网站| 国产精品不卡视频一区二区| 99久久中文字幕三级久久日本| 婷婷色av中文字幕| 国内精品美女久久久久久| 免费观看性生交大片5| 可以在线观看毛片的网站| 久久精品夜夜夜夜夜久久蜜豆| 成人性生交大片免费视频hd| 亚洲色图av天堂| 丰满乱子伦码专区| av线在线观看网站| av国产免费在线观看| 亚洲欧美日韩卡通动漫| 成人午夜精彩视频在线观看| 草草在线视频免费看| 国产亚洲一区二区精品| 身体一侧抽搐| 卡戴珊不雅视频在线播放| 国产在线男女| 亚洲国产精品专区欧美| 国产亚洲一区二区精品| av在线蜜桃| 性插视频无遮挡在线免费观看| 尾随美女入室| 久久久成人免费电影| 青春草国产在线视频| 1000部很黄的大片| 两个人视频免费观看高清| 80岁老熟妇乱子伦牲交| 久久精品久久久久久噜噜老黄| 97精品久久久久久久久久精品| 精品不卡国产一区二区三区| 一级毛片电影观看| 真实男女啪啪啪动态图| av免费观看日本| 最近2019中文字幕mv第一页| 欧美日韩综合久久久久久| 国产在线男女| 国产午夜精品论理片| 国产一区亚洲一区在线观看| 婷婷色av中文字幕| 成人毛片60女人毛片免费| 欧美成人a在线观看| 在线播放无遮挡| 国产亚洲91精品色在线| 成人亚洲精品av一区二区| 婷婷色麻豆天堂久久| 性色avwww在线观看| 亚洲在久久综合| 国产成人精品婷婷| 亚洲精品国产av蜜桃| freevideosex欧美| 久久久色成人| 黄色配什么色好看| 国产国拍精品亚洲av在线观看| av黄色大香蕉| 少妇丰满av| 三级国产精品片| 大话2 男鬼变身卡| 国产精品久久久久久av不卡| 尾随美女入室| 九色成人免费人妻av| 国产综合懂色| 男人舔女人下体高潮全视频| 色视频www国产| 国产免费一级a男人的天堂| 卡戴珊不雅视频在线播放| 尾随美女入室| 亚洲国产欧美人成| 亚洲在线自拍视频| 国产不卡一卡二| 久久精品久久久久久久性| 日韩人妻高清精品专区| 午夜免费激情av| 亚洲经典国产精华液单| 免费大片18禁| 亚洲av中文字字幕乱码综合| 在线免费观看的www视频| 熟妇人妻不卡中文字幕| 国产 一区精品| 久久久久性生活片| 三级国产精品欧美在线观看| 欧美成人一区二区免费高清观看| 成人高潮视频无遮挡免费网站| 激情 狠狠 欧美| 免费人成在线观看视频色| 国产精品一区二区三区四区久久| 又爽又黄a免费视频| 最新中文字幕久久久久| 如何舔出高潮| 九色成人免费人妻av| 亚洲av不卡在线观看| 亚洲精品日韩av片在线观看| 免费在线观看成人毛片| 男女啪啪激烈高潮av片| 少妇裸体淫交视频免费看高清| 国产久久久一区二区三区| 国产午夜精品久久久久久一区二区三区| 免费观看无遮挡的男女| 联通29元200g的流量卡| 日日撸夜夜添| 亚洲av在线观看美女高潮| 国内揄拍国产精品人妻在线| 女人被狂操c到高潮| 日产精品乱码卡一卡2卡三| 午夜激情福利司机影院| 久久精品人妻少妇| 熟女人妻精品中文字幕| 日本黄大片高清| 久久久久网色| 1000部很黄的大片| 国产精品嫩草影院av在线观看| 亚洲av二区三区四区| 少妇猛男粗大的猛烈进出视频 | 国产激情偷乱视频一区二区| 97超视频在线观看视频| 国产在视频线在精品| 日本wwww免费看| 极品教师在线视频| 哪个播放器可以免费观看大片| 99热这里只有是精品在线观看| 亚洲精品亚洲一区二区| 五月天丁香电影| 午夜福利网站1000一区二区三区| 日本一二三区视频观看| 国产精品1区2区在线观看.| 久久6这里有精品| 亚洲最大成人中文| 天美传媒精品一区二区| 尾随美女入室| 乱人视频在线观看| 在线播放无遮挡| 综合色av麻豆| kizo精华| 超碰av人人做人人爽久久| 国产午夜精品论理片| 乱人视频在线观看| 亚洲无线观看免费| 青春草亚洲视频在线观看| 小蜜桃在线观看免费完整版高清| 狂野欧美激情性xxxx在线观看| 亚洲av中文字字幕乱码综合| 国产一区二区亚洲精品在线观看| 成人一区二区视频在线观看| 床上黄色一级片| 啦啦啦韩国在线观看视频| 亚洲国产欧美人成| 欧美潮喷喷水| 国产高清国产精品国产三级 | 乱人视频在线观看| 午夜视频国产福利| 搡老乐熟女国产| 尾随美女入室| 成人漫画全彩无遮挡| 麻豆成人av视频| 天美传媒精品一区二区| 人妻制服诱惑在线中文字幕| 毛片一级片免费看久久久久| 女人久久www免费人成看片| 男女视频在线观看网站免费| 国产黄片美女视频| 日韩欧美精品免费久久| 亚洲欧美一区二区三区黑人 | 国产 亚洲一区二区三区 | 尤物成人国产欧美一区二区三区| 亚洲av.av天堂| 国产免费又黄又爽又色| 国产亚洲午夜精品一区二区久久 | 特级一级黄色大片| 卡戴珊不雅视频在线播放| 亚洲色图av天堂| 国产亚洲精品久久久com| 91午夜精品亚洲一区二区三区| 女人被狂操c到高潮| 亚洲怡红院男人天堂| 黄色欧美视频在线观看| 中文精品一卡2卡3卡4更新| 久久这里有精品视频免费| 在线观看美女被高潮喷水网站| 免费无遮挡裸体视频| 2021少妇久久久久久久久久久| 麻豆久久精品国产亚洲av| 日韩 亚洲 欧美在线| 欧美日韩在线观看h| 极品教师在线视频| 久久久成人免费电影| 99热这里只有是精品50| 丰满乱子伦码专区| 亚洲欧美日韩东京热| 国精品久久久久久国模美| 免费av不卡在线播放| 久久国内精品自在自线图片| 国产高潮美女av| 国产69精品久久久久777片| 亚洲欧美成人综合另类久久久| 中文字幕人妻熟人妻熟丝袜美| 97超视频在线观看视频| 联通29元200g的流量卡|