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

    Design optimization of composite egg-shaped submersible pressure hull for minimum buoyancy factor

    2022-01-05 09:41:04MuhmmdImrnDongynShiLiliTongHfizMuhmmdWqsMuqeemUddin
    Defence Technology 2021年6期

    Muhmmd Imrn , Dongyn Shi ,*, Lili Tong , Hfiz Muhmmd Wqs ,Muqeem Uddin

    a College of Mechanical and Electrical Engineering Harbin Engineering University, Harbin,150001, China

    b College of Aerospace and Civil Engineering Harbin Engineering University, Harbin,150001, China

    c Pakhtunkhwa Energy Development Organization (PEDO), Khyber Pakhtunkhwa, Pakistan

    Keywords:Composite egg-shaped pressure hull Design optimization Buoyancy factor Material failure Buckling instability

    ABSTRACT This paper describes a design optimization study of the composite egg-shaped submersible pressure hull employing optimization and finite element analysis(FEA)tools as a first attempt to provide an optimized design of the composite egg-shaped pressure hull for manufacturing or further investigations. A total of 15 optimal designs for the composite egg-shaped pressure hull under hydrostatic pressure are obtained in terms of fibers’ angles and the number of layers for 5 lay-up arrangements and 3 unidirectional (UD)composite materials. The optimization process is performed utilizing a genetic algorithm and FEA in ANSYS. The minimization of the buoyancy factor (B?F) is selected as the objective for the optimization under constraints on both material failure and buckling strength. Nonlinear buckling analysis is conducted for one optimal design considering both geometric nonlinearity and imperfections. A sensitivity study is also conducted to further investigate the influence of the design variables on the optimal design of the egg-shaped pressure hull.

    1. Introduction

    Owing to their excellent characteristics such as increased strength-to-weight ratio, good corrosion resistance, and low magnetic and acoustic signatures, fiber reinforced plastics (FRP) composite materials are increasingly used in the production of highspeed crafts, mine countermeasure platforms and numerous subparts of big naval ships. FRP composites are also utilized in the manufacturing of submersible pressure hulls and the pressure hulls manufactured with composite materials have larger collapse depth than the pressure hulls manufactured with steels and other alloys[1,2]. Composite pressure hulls are constructed with composite laminates using different lay-up sequences and material systems.The design of composite pressure hulls greatly depends on the number of ply’s layers, fiber’s angles and material systems. It is of very much interest to researchers and design engineers to optimize the design of composite pressure hull before going into the manufacturing process. An optimized design will further reduce the overall weight of the composite pressure hull and hence the cost will also decrease. In the field of design optimization, a lot of research has been performed concerning the composite constructions. A detailed review of the optimization of various composite structures has been reported in Ref.[3].Several studies have looked into the optimization of the plates modeled with composites for different objectives by means of a genetic algorithm (GA) and numerical study[4-12].A considerable amount of research has been performed to obtain the optimum designs of the pressure hulls constructed with composites by means of optimization tools combined with numerical analysis subjected to instability constraint or both instability and material failure constraints. The research concerning the optimization of the composite submersible pressure hulls has been performed for diverse objectives such as weight reduction, minimization of the buoyancy factor and maximization of the buckling load multiplier. Different types of hull configurations have been utilized in the studies on the optimization of the pressure hull such as single skin cylindrical shell configuration [13-22], elliptical hull configuration [23-25] and sandwich core configuration[26-28].Due to the useful properties such as reduced buoyancy factor and effective distribution of stress and displacement, the spherical hulls manufactured with steels or other alloys have been used for deep-sea manned submersible platforms[29].The research conducted about the spherical hulls is mostly related to linear and non-linear buckling analyses. Both geometric and material non-linearity have been considered in the buckling analysis of spherical hulls [30-38]. Although spherical hulls have some advantages, there are still important issues of concern to designers such as highly geometric imperfections, material plasticity sensitivity, and less efficient space utilization. Due to the aforementioned issues, researchers have conducted studies to develop and analyze non-spherical hulls such as egg-shaped pressure hulls. Eggshell is one of the most effective structures among the naturally occurring shells. The most important properties of the eggshell are high strength-to-weight ratio,high span-tothickness ratio, good stability, and good material distribution[29,39]. According to the available literature, there has been no research concerning the optimization of the composite egg-shaped pressure hull. Most of the studies have been conducted on the linear and non-linear buckling analyses of the egg-shaped pressure hulls constructed with steels or other alloys. The designs of those egg-shaped pressure hulls under external pressure have been compared with spherical shape hulls constructed with the same materials [29,40-43].

    The most commonly used composite materials in the modeling and manufacturing of pressure hulls in previous studies are Epoxy/Carbon [14,16,17,20,22,24,25,28,44,45], Epoxy/Glass [14,26,28,45],Epoxy/Boron [24,26,45] and Epoxy/Graphite [24,26]. Owing to the anisotropy of the composites, their properties can be adjusted by varying the angles of fibers in every laminate’s ply [46]. Studies conducted on composite pressure hull have used several angle configurations for composite lay-up such as [±55] [14,28], [0/90]and [90/0] [24,28], [0/±45/90] [23,28] and [0/90/0],[0/-10/90/-10/90] [28] and [φ1/φ2/φ3/φ4/φ5] where φibelonged to 30°, 45°, 60°,75°and 90°[14]. The angle configurations mean the orientation angles of different layers of the laminate. The orientation angle is measured with respect to a reference direction or axis of the laminate. The present paper uses 5 different lay-up arrangements comprising [0a/90b/0c], [10a/-10b/90c/-10d/10e], [β1a/β2b], [β1a/β2b/β3c] and [β1a/β2b/β3c/β4d/β5e] for creating the composite laminates.Here a, b and c are the number of layers of different plies in the laminate.The same lay-up arrangements have also been utilized in the research reported in the present paper. The lay-up arrangements [0a/90b/0c] and [10a/-10b/90c/-10d/10e] are chosen since these can be generated by means of the hand lay-up procedure employing pre-preg composites [47], utilized in the present research.The third lay-up arrangement[β1a/β2b]is chosen to know whether a lay-up with only two fiber’s angles can be optimized or not.The lay-up arrangements[β1a/β2b/β3c]and[β1a/β2b/β3c/β4d/β5e]are chosen for comparison with pre-established lay-up arrangements [0a/90b/0c] and [10a/-10b/90c/-10d/10e] and to determine whether any further optimization can be obtained for preestablished fiber’s angles.

    Recently egg-shaped pressure hulls have been demonstrated to have advantages over other shape hulls.However,these hulls have been constructed with steels or other alloys and the egg-shaped pressure hulls constructed with composite materials have not been fully investigated yet. Moreover, due to the big change in its circumferential and meridional mechanical performance, eggshaped pressure hull constructed with composite materials may give better results as compared to the egg-shaped hull constructed with steel and other alloys [48]. Therefore, it is very important to conduct an optimization study on the composite egg-shaped pressure hull. The present paper describes the design optimization study of the submerged composite egg-shaped pressure hull to obtain the minimum buoyancy factor (B?F) by means of a GA and FEA as a first attempt to provide a base for future experimental study. The design optimization study is performed for 5 lay-ups arrangements and 3 UD composites, Epoxy/Carbon, Epoxy/Glass,and Epoxy/Boron. The safety factors of Tsai-Wu and Tsai-Hill composite failure criteria and the buckling load multiplier are employed as the material failure and instability constraints respectively on the optimization process. The optimization simulations are performed based on linear elastic material models and perfect structures. A minimum safety factor of 1.5 is set for both composite failure criteria and buckling load multiplier.A sensitivity investigation is also carried out to further explore the influence of the design variables on the optimum design. To investigate the influence of initial geometric imperfection on the optimal design,non-linear buckling analysis is also conducted for one optimum design.

    2. Geometric features of the egg-shaped pressure hull

    The contour of the egg-shaped pressure hull is illustrated in Fig. 1, which is based on the goose-egg shape. The radius of circumference R(x), of the shell of revolution, is calculated by Eq.(1). This equation was originally developed by Narushin [49] and later adopted by researchers in Refs. [29,40-43].

    Fig.1. Geometry of the egg-shaped pressure hull.

    The surface area of the egg-shaped pressure hull is determined by Eq. (4) [29,40].

    In this paper, the volume of the egg-shaped pressure hull is considered to be equal to that of the composite cylindrical pressure hull having hemispherical ends used in Refs.[28,51].The previously studied cylindrical pressure hull was 11.04 m long and radii of the hemispherical ends were 1.370 m. The length of the major and minor axis is determined from Eq.(3)and Eq.(2),using the volume of the aforementioned cylindrical pressure hull and previously used value of the shape index (SI).The value of previously used shape index (SI) is 0.69 [29,40,43]. According to Eq. (3) and Eq. (2), the length of the major axis(L)is 6.199 m and the length of the minor axis (B) is 4.2773 m.

    3. Finite element model of the egg-shaped composite pressure hull

    The geometric model of the egg-shaped pressure hull is created with the help of the revolve feature in ANSYS DesignModeler.First,a curve is sketched with values obtained from Eq.(1).The curve is then revolved through 360°to get the geometric model of the pressure hull.The Shell elements of ANSYS having an approximate size of 100 mm are used to create the mesh model for further analysis. Fig. 2(a) represents the mesh model of the pressure hull.The study performed in the present paper is a part of our ongoing design optimization study of different types of composite pressure hulls of the same volumes including cylindrical, egg-shaped and spherical pressure hulls.The procedure adopted for the creation of the finite element model and carrying out subsequent analysis is essentially the same as described in our previous study on the cylindrical composite pressure hull [51].

    3.1. Composite lay-up creation

    The ANSYS Composite PrePost (ACP) is used for the creation of composite lay-up of the pressure hull.The orthotropic properties of the 3 UD composites,Epoxy/Carbon,Epoxy/Glass,and Epoxy/Boron are listed in Table 1. To decrease the computation duration, the thickness of the ply is considered to be equal to 1 mm. In earlier studies, an equivalent size of the ply thickness has been used[28,51].In the ACP,first,a Fabric is defined and then the Orientation and Reference Direction of fiber are generated with the help of their respective features in ACP. The plies for each composite’s layer are defined in the Modeling Groups of the ACP.

    3.2. Boundary and loading conditions

    In order to avoid rigid body motion, boundary conditions are applied to the forward and aft ends of the egg-shaped pressure hull.Fig. 3 shows the boundary and loading conditions of the eggshaped pressure hull. The outside surface of the pressure hull is applied with 3 MPa uniformly distributed pressure (A). At the forward end, the translation of one node in the Y, and Z directions is restricted(B).At the aft end,the translation of one node is restricted in the X,Y, and Z directions (C). Additionally, at the aft end, the rotation of the same node about the X axis is restricted(D).Similar boundary conditions have also been used in previous studies[28,51].

    Fig. 2. (a) Mesh model (b) Sectional view of the mesh model.

    The 3 MPa pressure has been arbitrarily chosen. In the linear analysis, a minimum factor of safety of 1.5 is selected for both material failure and buckling instability.This means that the design pressure of the pressure hull is 3 MPa(300 m of water depth).The FEA requires pressure or other loading conditions to complete and even greater value of outside pressure can be applied. However, it will increase the required thickness of the pressure hull and hence the computational time and cost will increase. Therefore, 3 MPa external pressure has been used in all the optimization simulations presented in this study.

    3.3. Procedure for finite element model solution

    The ANSYS Composite PrePost (ACP) is utilized to create the composite model of the egg-shaped pressure hull. In the ACP, parameters are defined for the fibers’ angles and the corresponding number of layers.The data containing the composite shell model is then linked to the Static Structural module of ANSYS to conduct the Static Structural analysis. In the Static Structural module, the FE model is solved and parameters are defined for the safety factors of the composite failure criteria. A number of earlier research works that focused on the same research field have adopted Tsai-Wu[14,17,19,20,23-26,28] and Tsai-Hill [16,28] failure criteria of composites.To conduct the linear Buckling analysis the FE model is linked to the Eigenvalue Buckling module.In the same module,the model is solved and the parameter for the linear buckling load multiplier is defined. All the three modules are then connected to the Direct Optimization module of the ANSYS’s Design Explorer.In the same module, values for all pre-established parameters and constraints are entered and the objective is defined. Finally, the project is run and the optimum design candidates are calculated.A parametric design analysis of a reference egg-shaped HY100 steelhull meshed with shell elements is conducted in ANSYS for calculation of the percent decrease in the buoyancy factor (B?F). The density and the ultimate yield strength of the HY100 steel are 7828 kg/m3and 790 MPa respectively[52,53].

    Table 1 Orthotropic material properties of Epoxy/Carbon, Epoxy/Glass and Epoxy/Boron[24,54].

    Fig. 3. Loading and boundary conditions of the egg-shaped composite submersible pressure hull,A:3 MPa external pressure,B:At forward end,Y and Z displacements are constrained, C: At aft end, X, Y and Z displacements are constrained, D: At aft end rotation about X axis is constrained.

    4. Procedure for design optimization

    The optimization process is conducted in the Direct Optimization module of ANSYS. The Controlled Elitist GA employed in this study, is an alternate form of the Non-dominated Sorting Genetic Algorithm-II(NSGA-II).It utilizes a number of features to determine the Global Optimum. In order to use the optimization module,atleast one objective is required.However, multiple objectives can also be defined.The optimization procedure flow chart is shown in Fig. 4. The main features of the GA are given below[55,56].

    Optimal Space-Filling:The Optimal Space-Filling(OSF)is a Latin Hypercube Sampling (LHS) with post-processing ability. More details about the Space-Filling property can be found in Ref. [57,58].The OSF achieves a more uniform distribution of points. In this study,the OSF is used for the generation of initial sampling.The size of the initial population and the number of samples per iteration are selected based on the number of input parameters.In this study,the size of the initial population is equal to ten times the number of input parameters and the number of samples per iteration is selected to be greater than half of the size of the initial population in the case of every lay-up arrangement studied in the present work.

    Fig. 4. Procedure flow chart for design optimization.

    Crossover: The Crossover generates new chromosomes (children)by uniting two old chromosomes(parents).According to the principle of evolution, when the finest characteristics of each parent are chosen then the children can be more enhanced as compared to their parents. Three different types of Crossover approaches including One-point, Two-points, and Uniform are available. In the present work, the One-point Crossover approach with 0.98 crossover probability is used.A One-point Crossover randomly chooses a Crossover point and interchanges the two parent chromosomes to generate two new offspring. A two-points crossover randomly chooses two crossover points and interchanges the two parent chromosomes to generate two new offspring. A uniform crossover mixes the parent chromosomes at the gene level rather than the segment level (as in the One and Two-point Crossover).

    Mutation: The Mutation is a genetic phenomenon that alters the values of genes within a genome (chromosomes) from their original values. Due to the incorporation of different gene values,the GA may arrive at improved outcomes than earlier produced.The mutation is helpful in preventing the outcomes from becoming stagnant at the local optimum.In this study,a mutation probability of 0.01 is used.

    Maximum Allowable Pareto Percentage Criterion: This is a convergence criterion that checks the percentage of a predetermined ratio of Pareto points-to-number of samples in one iteration. The optimization converges when the pre-determined percentage is reached. In this study, value of 70% is used for the Maximum Allowable Pareto Percentage.

    Convergence Stability Percentage Criterion: This convergence criterion checks the stability of the population by using the mean and standard deviations of the output variables. Starting from the second generation,the mean and standard deviations are checked.The optimization converges, when these both deviations are smaller than the pre-determined value.In this study,value of 2%is used for Convergence Stability Percentage.

    4.1. Objective function of optimization

    The design optimization study of the composite egg-shaped pressure hull is performed based on minimizing the buoyancy factor (B?F) as objective function subjected to material failure and instability constraints.The same objective has also been considered in earlier studies conducted on optimization of the composite pressure hulls [23-25,51,65]. The Buoyancy factor (B?F) is calculated by dividing the weight of the hull by the weight of the water displaced by the hull (Eq. (5)).

    4.2. Design variables (parameters) of optimization

    Design variables are the objects that can change the properties of the model within a pre-defined range during the design optimization process. In the present paper, the number of plies’ layers and the fiber’s angle in each layer are used as design variables for the optimization study.

    4.3. Constraints on optimization

    Constraints are the limits imposed on various design parameters during a constrained optimization problem. Three types of constraints,material failure,structural instability,and side constraints are used in the present optimization study. These constraints are explained as follows.

    Constraints on material failure: Safety factors of the Tsai-Wu(FSTW) and Tsai-Hill(FSTH) composite materials failure criteria are considered as constraints on material failure. To keep the material of the pressure hull safe,the values of both FSTWand FSTHmust be higher than or equal to 1. Similarly, for the steel hull, Von Mises stress(σv)must be smaller than or equal to the yield strength(σy)of the steel.Eqs. (6)-(9) represent these constraints.

    Constraint on structural Instability: The buckling load multiplier (λ) is considered as a constraint on the structural instability and must be higher than or equal to1 in order to keep the composite pressure hull stable. It is calculated by dividing the critical buckling pressure(P) by the applied hydrostatic pressure (Po).

    The instability constraint to be equal to or greater than 1, is selected based on the assumption that the material of the composite pressure hull is linear elastic and free from any imperfection.This is the minimum requirement when the material is linear elastic and the pressure hull is defect free.But in the real situation,the structure may not remain linear elastic and free of imperfection due to a number of reasons. Therefore, in the simulations conducted in this paper, a minimum safety factor of 1.5 is set for both material and instability constraints.

    Side constraints: The upper and lower bounds of different parameters are known as side constraints. In the present design optimization study, the angles of orientation of fibers and the number of plies (layers) in the laminate are the side constraints.These side constraints are represented by Eqs.10 and 11.

    In case of fiber’s orientation angles,

    5. Results and discussions

    First, a parametric study of the reference egg-shaped HY100 steel hull is carried out under 3 MPa hydrostatic pressure, for comparison of the percent decrease (DP) in the buoyancy factor.The parametric study aims at minimization of the buoyancy factor(B?F), while keeping Von Mises stress (σv) and buckling load multiplier(λ)as the material failure and instability constraints.The results of the parametric study for steel hull are given in Table 2.As can be seen from this Table,with the optimum thickness of 18 mm,the buoyancy factor(B?F)equals 0.17306.The Von Mises stress(σv)is 248.3 MPa and the buckling load multiplier (λ) is 2.1872. The results of the design optimization of the composite egg-shaped submerged pressure hull for all considered lay-ups and material systems are listed in Tables 3-7.

    In the case of all lay-up arrangements studied in the present work, 20 maximum number of iterations is used. However, in all cases the optimum designs were obtained based on the Maximum Allowable Pareto Percentage and Convergence Stability Percentage criteria before reaching the maximum number of iterations. The optimization simulations are conducted on a workstation with specifications as follows: Intel (R) Xeon (R) CPU E5-2620 V2 @2.10 GHz(24 CPU),32 GB RAM and 2 TB HD.The minimum number of simulations required for convergence of one optimization case is approximately 250-600 depending on the number of input parameters. The average computational time for one complete optimization case is approximately 10 h. The convergence of the objective function for Epoxy/Glass lay-up [β1a/β2b/β3c] is shown in Fig. 5. In the case of this lay-up, both 70% Maximum Allowable Pareto Percentage and 2%Convergence Stability Percentage criteria are fulfilled.

    5.1. Design optimization of [0a/90b/0c]

    The results of the design optimization of lay-up [0a/90b/0c] are listed in Table 3.In the same lay-up arrangement,a,b and c denote the number of respective plies’layers and each one can take value from 1 to 75. The results of the design optimization of lay-up [0a/90b/0c] shows that the minimum buoyancy factor (B?F)of 0.079 is computed for Epoxy/Carbon and hence demonstrates a maximum decrease(DP)of 54.2%among all the three considered composites.For the same material, the safety factors (FSTW) and (FSTH) and buckling load multiplier (λ) are computed as 3.02, 2.17 and 1.81 respectively.For the same lay-up arrangement,the buoyancy factor(B?F) of 0.13 is computed for Epoxy/Glass and demonstrates a decrease(DP)of 24.78%.For the same composite,FSTW,FSTH,and λ are computed as 2.6, 2.04 and 1.88 respectively. Similarly, the buoyancy factor (B?F) of 0.093 is computed for Epoxy/Boron and results in a decrease (DP) of 46.07%. For the same composite, the values of FSTW, FSTH, and λ are computed as 4.93, 4.16 and 1.81 respectively. The values of FSTWand FSTHfor Epoxy/Boron are higher than their respective values for both Epoxy/Carbon and Epoxy/Glass.

    5.2. Design optimization of [10a/-10b/90c/-10d/10e]

    The results of design optimization for lay-up[10a/-10b/90c/-10d/10e]are given in Table 4.In the same lay-up arrangement,a,b,c,d,and e denote the number of respective plies’ layers and each one can take value from 1 to 75.The minimum buoyancy factor(B?F)of 0.085 is computed for Epoxy/Carbon and hence gives the maximum decrease(DP)of 51.01%.For the same composite,FSTWand FSTHand λ are computed as 2.4334, 1.7532 and 1.7824 respectively. The buoyancy factor (B?F) for Epoxy/Glass is computed as 0.125 with DP=27.62%.The values of FSTW,FSTHand λ for Epoxy/Glass are 2.26,1.98 and 1.73 respectively. Similarly, the buoyancy factor (B?F) of 0.106 is computed for Epoxy/Boron and demonstrates a decrease(DP)of 38.97%.For the same composite,the values of FSTWand FSTHare 3.36 and 3.33 respectively.The optimized value of λ for Epoxy/Boron is computed as 2.61,which is larger than its respective values for both Epoxy/Glass and Carbon/Epoxy. Moreover, the decrease(DP)of Epoxy/Boron is larger than that of Epoxy/Glass and smaller than that of Epoxy/Carbon.

    5.3. Design optimization of [β1a/β2b]

    For lay-up[β1a/β2b]both the number of layers and fibers’angles are optimized. In lay-up [β1a/β2b], β1and β2represent the fibers’angles and “a” and “b” represent the number of respective ply’s layers. The angles β1and β2can take values from 0°to 90°, while the numbers of ply’s layers a and b can take value from 1 to 75.The results of design optimization for lay-up [β1a/β2b] are listed in Table 5. The minimum buoyancy factor (B?F)of 0.081 is computed for Epoxy/Boron and therefore results in the highest decrease(DP)of 53.17%.The optimized values of FSTW,FSTHand λ for Epoxy/Boron are 2.04, 1.87 and 1.84 respectively. For Epoxy/Carbon, the buoyancy factor(B?F)of 0.099 is computed with a decrease DP=42.49%.The optimized values of FSTW,FSTHand λ for Epoxy/Carbon are 1.82,1.65 and 2.92 respectively. For the same lay-up arrangement, the buoyancy factor (B?F) of 0.130 is computed for Epoxy/Glass which gives a decrease (DP) of 24.78%. The values of FSTW, FSTHand λ for Epoxy/Glass are 1.80, 1.6 and1.78 respectively. The value of λ for Epoxy/Glass is larger than its value for Epoxy/Carbon and lesser than its value for Epoxy/Boron.

    Table 2 Parametric design analysis of egg-shaped submerged steel hull.

    Table 3 Design optimization of lay-up [0a/90b/0c].

    5.4. Design optimization of [β1a/β2b/β3c]

    For lay-up [β1a/β2b/β3c] design optimization is carried out both in terms of fibers’ angles and the number of layers. In lay-up [β1a/β2b/β3c], β1, β2, and β3represent the fiber’s angles and a, b, and c represent the number of respective ply’s layers. The angles β1, β2and β3can take values from 0°to 90°, while the number of ply’s layers a, b and c can take value from 1 to 75. The results of the optimization for this lay-up are listed in Table 6.In the case of this lay-up arrangement,the minimum buoyancy factor(B?F)of 0.081 is computed for Epoxy/Boron and hence gives the highest decrease(DP)of 53.17%.The values of FSTW,FSTH,and λ for Epoxy/Boron are 1.85 and 1.58 and 1.78 respectively. The buoyancy factor (B?F)of 0.085 is computed for Epoxy/Carbon with the decrease DP=51.01%.The values of FSTW,FSTH, and λ for Epoxy/Carbon are 2.04 and 1.76 and 2.49 respectively.The value of λ for Epoxy/Carbon is larger than its values for both Epoxy/Boron and Epoxy/Glass. Similarly, for Epoxy/Glass the buoyancy factor (B?F) of 0.12 is computed with DP=29.04%.The values of FSTW,FSTH,and λ for Epoxy/Glass are 1.68 and 1.58 and 1.67 respectively.The buoyancy factor(B?F)for lay-up[β1a/β2b/β3c] is 6.97% larger than its respective value in lay-up [0a/90b/0c] for Carbon/Epoxy. The FSTWand FSTHfor Epoxy/Carbon in lay-up [0a/90b/0c] are 47.95% and 22.95% larger than their respective values in lay-up [β1a/β2b/β3c], while its value of λ is 37.86%larger in lay-up[β1a/β2b/β3c]than its value in lay-up[0a/90b/0c].The overall efficiency of Epoxy/Carbon is better in lay-up [0a/90b/0c]than in lay-up [β1a/β2b/β3c]. The buoyancy factor (B?F) for Epoxy/Glass in lay-up [0a/90b/0c] is 6% larger than its respective value in lay-up[β1a/β2b/β3c].The FSTW,FSTHand λ for Epoxy/Glass in lay-up[0a/90b/0c]are 55.04%,28.86%and 12.1%larger than their respective values in lay-up[β1a/β2b/β3c].This shows that with a slight increase in weight,Epoxy/Glass performed well in case of lay-up[0a/90b/0c].The buoyancy factor (B?F) of Epoxy/Boron in lay-up [0a/90b/0c] is 15.15% larger than its respective value in lay-up [β1a/β2b/β3c]. TheFSTW,FSTHand λ for Epoxy/Boron in lay-up[0a/90b/0c]are 167.08%,163.67%and 1.95%larger than their respective values in lay-up[β1a/β2b/β3c].The material failure performance of Epoxy/Boron in lay-up[0a/90b/0c] is better than its performance in lay-up [β1a/β2b/β3c].However,Epoxy/Boron has considerably more weight in lay-up[0a/90b/0c] as compared to its weight in lay-up[β1a/β2b/β3c].

    Table 4 Design optimization of lay-up [10a/-10b/90c/-10d/10e].

    Table 5 Design optimization of lay-up [β1a/β2b].

    Table 6 Design optimization of lay-up [β1a/β2b/β3c].

    5.5. Design optimization of [β1a/β2b/β3c/β4d/β5e]

    Fig. 5. Convergence of the optimization for Epoxy/Glass lay-up [β1a/β2b/β3c].

    For lay-up [β1a/β2b/β3c/β4d/β5e], both the number of layers and fibers’ angles are optimized. In lay-up [β1a/β2b/β3c/β4d/β5e], β1, β2,β3, β4, and β5represent the fiber’s angles and a, b, c, d, and e represent the number of respective ply’s layers.The angles β1,β2,β3,β4,and β5can take values from 0°to 90°,while the number of ply’s layers a, b, c, d, and e can take value from 1 to 75. The results of design optimization for lay-up [β1a/β2b/β3c/β4d/β5e] are listed in Table 7. The minimum buoyancy factor (B?F) of 0.098 is computed for Epoxy/Boron and therefore results in a highest decrease(DP)of 41.8%. The values of FSTWand FSTHfor Epoxy/Boron are 2.44 and 1.99 respectively and these values are larger than the respective values of FSTWand FSTHfor Epoxy/Glass(1.88 and 1.97)and smaller than those of Epoxy/Carbon(2.58 and 2.31).The λ of Epoxy/Boron is 3.035 and it is larger than its respective value of λ for Epoxy/Glass(1.6) and smaller than that of Epoxy/Carbon (3.85). The buoyancy factor (B?F) in lay-up [10a/-10b/90c/-10d/10e] for Epoxy/Boron is 7.5% larger than its respective value in lay-up [β1a/β2b/β3c/β4d/β5e].The FSTWand FSTHin lay-up [10a/-10b/90c/-10d/10e] for Epoxy/Boron are 37.55%and 67.95%larger than their respective values in layup [β1a/β2b/β3c/β4d/β5e]. The λ for Epoxy/Boron in lay-up [β1a/β2b/β3c/β4d/β5e]is 16.07%larger than its respective value in lay-up[10a/-10b/90c/-10d/10e]. The buoyancy factor (B?F) in lay-up [10a/-10b/90c/-10d/10e] for Epoxy/Glass is 15.7% larger than its respectivevalue in lay-up [β1a/β2b/β3c/β4d/β5e], which equals 0.14 with the decrease DP=16.27%.The values of FSTWand FSTHin lay-up[10a/-10b/90c/-10d/10e] for Epoxy/Glass are 15.5% and 0.66% larger than their respective values in lay-up [β1a/β2b/β3c/β4d/β5e]. The λ for Epoxy/Glass in lay-up [10a/-10b/90c/-10d/10e] is 8.25% larger than its respective value in lay-up [β1a/β2b/β3c/β4d/β5e]. The buoyancy factor(B?F)in lay-up[β1a/β2b/β3c/β4d/β5e]for Epoxy/Carbon is equal to 0.12 with decrease DP= 30.77%, which is 41.31% larger than its respective value in lay-up [10a/-10b/90c/-10d/10e]. The FSTW, FSTHand λ in lay-up [β1a/β2b/β3c/β4d/β5e] for Epoxy/Carbon are 6.2%,31.98% and 115.8% larger than their respective values in lay-up[10a/-10b/90c/-10d/10e].

    Table 7 Design optimization of lay-up [β1a/β2b/β3c/β4d/β5e].

    5.6. Comparison of different types of composite pressure hulls

    The results of the optimization of the present study are compared with our previous studies on the design optimization of the composite cylindrical and spherical pressure hulls [51,59]. To have equal basis for comparisons, the volumes, lay-ups, materials,mesh sizes, boundary and loading conditions, and optimization procedures are kept the same for all composite pressure hulls.The results of comparison are listed in Tables 8-10.It can be noted from these tables that the composite cylindrical pressure hull has higher buoyancy factors as compared to the composite spherical and eggshaped pressure hulls for all lay-ups and materials considered.This means that composite cylindrical pressure hull has more weight as compared to the composite spherical and composite egg-shaped pressure hulls for the same design pressure. Similarly, it can also be noted from Tables 8-10 that the buoyancy factors of the composite spherical and egg-shaped pressure hulls are almost the same in case of all lay-ups and materials considered.This means that the composite spherical and egg-shaped pressure hulls have almost the same weights for the same design pressure. The safety factors of failure criteria for composite cylindrical pressure hull are higher than the corresponding safety factors of failure criteria for both composite spherical and egg-shaped pressure hulls, whereas the buckling strength factor of composite cylindrical pressure hull is lower than its corresponding values for the composite spherical and egg-shaped pressure hulls. From this comparison it can be concluded that if composite egg-shaped pressure hull or composite spherical pressure hull are manufactured with accuracy and the imperfections are avoided,then these pressure hulls will reduce the manufacturing as well as operating costs to a greater extent.

    6. Non-linear buckling analysis

    All the simulations performed in this study for the optimization of composite egg-shaped pressure hull use linear buckling analysis,which can be employed to assess the initial design of the pressure hull.Linear Static Structural and Eigenvalue Buckling analyses have also been employed in a number of earlier studies on the design optimization of pressure hull modeled with composites[14,16,17,19,20,23-26,28].In linear buckling analysis,the structures are assumed to be elastic and free of any imperfection.However,in reality, the structures are not perfectly elastic and also have imperfections. Moreover, at increased depth geometric as well material non-linearity may occur.Therefore,it is necessary to carry out non-linear buckling analysis in addition to linear buckling analysis to accurately model the behavior of structures undergoing buckling deformation. For steel hulls, non-linear buckling analysis is performed using elastic-plastic material properties of steel and geometric non-linearity with and without imperfection. In this process, the steel hull is loaded with external pressure until the material yields, and the hull buckles and collapses completely[28,60]. In the present study, for carrying out non-linear buckling analysis of the egg-shaped pressure hull, optimized Epoxy/Carbon lay-up [02/9024/017] is used. Similarly, for comparison, non-linear buckling analysis of the previously optimized composite spherical pressure hull having Epoxy/Carbon lay-up [012/9017/05] is also conducted[59].The first 12 linear buckling modes of the composite egg-shaped pressure hull obtained from ABAQUS are shown in Fig. 6. It can be observed that composite egg-shaped pressure hull has 2-4 horizontal half waves and several circumferential waves.The linear critical buckling pressure for the lowest mode is 5.44 MPa.

    In order to determine the accurate buckling behavior, imperfection is incorporated in the non-linear analysis of the composite egg-shaped pressure hull. Different methods of incorporating the imperfection into the analysis are used such as linear buckling mode-shaped imperfection(LBMI),the worst multiple perturbation load imperfection (WMPI), single boundary perturbation imperfection (SBPI), single load perturbation imperfection (SLPI),geometric dimple-shaped imperfection (GDI), mid-surface imperfection (MSI) and axisymmetric imperfection (ASI) [61,62]. In this paper, the linear buckling mode-shaped imperfection (LBMI)method is used for the investigation of the imperfection sensitivity of the composite egg-shaped pressure hull. This is a commonprocedure and has been used in previous studies for the determination of the imperfection sensitivity of the composite shells[28,61,63,64]. In the linear buckling mode-shaped imperfection(LBMI)method,first,a linear buckling analysis is performed for the perfect shell and the coordinates of Eigenmodes are obtained. In the second step,the worst buckling mode is applied with a scaling factor to the imperfect shell in a geometrically non-linear analysis by editing the keyword file. Generally, the first buckling mode is assumed to be the worst buckling mode. However, when the difference between the Eigenvalues is small,higher modes may result in smaller buckling pressure. Therefore, many buckling modes should be examined to find out the worst buckling mode [60]. In this paper,a total of 15 linear buckling modes are examined and the worst buckling mode is found out and employed in the non-linear buckling analysis of the composite pressure hulls.Similarly,a total of 6 different imperfections values including 1, 3, 4, 5 and 10 mm are added to the initial perfect geometric shapes of the pressure hulls.In ABAQUS,this can be accomplished by using the static RIKS analysis with NLGEOM feature. The pressure hull are loaded with 6 MPa external pressure. The external pressure is applied in 100 automatically controlled increments. For the lowest and worst mode imperfections, the critical buckling modes and the postbuckling modes obtained from the nonlinear buckling analysis related to all 6 imperfections are shown in Figs. 7 and 8 respectively. The post-buckling modes appear to be in the form of local dents, which are similar to the previously experimentally determined buckling modes of Epoxy/Carbon composite cylindrical pressure hulls in(Fig.9)[44,64].Similarly,for the lowest and worst mode imperfections, the load-displacement relationships for all 6 imperfections are plotted in Figs.10 and 11 respectively. In these figures,the points(a-f)and(g-l)correspond to the critical buckling modes and the points (A-F) and (G-L) correspond to the postbuckling modes of the pressure hull shown in Figs. 7 and 8. The load(PNL)in Figs.10 and 11,is the buckling pressure obtained from the non-linear buckling analysis and the displacement represents the inward displacement of a reference node located at the center of the longitudinal half wave of the buckling mode of the pressure hull in the non-linear buckling analysis. The reference node gives the maximum displacement. To compare the imperfection sensitivity of the composite egg-shaped and spherical hulls, the results are compared in terms of the ratio of buckling pressure obtained from non-linear analysis -to- buckling pressure obtained from linear analysis().These results are listed in Tables 11 and 12 for lowest mode imperfections and worst mode imperfections respectively. The imperfection size has a great effect on the buckling pressure of both pressure hulls. In the case of lowest mode imperfections,as can be observed from Table 11,theratio for the egg-shaped pressure hull is greater than that of the spherical hull.Similarly, in the case of worst mode imperfections, as can be seen from Table 12,theratio for the egg-shaped pressure hull is also greater than that of the spherical hull.This means that for the same imperfection size,the composite egg-shaped pressure hull exhibits lower decrease in the buckling pressure as compared to the composite spherical pressure hull. In other words, the imperfection sensitivity of the composite egg-shaped pressure hull is lower than that of the composite spherical pressure hull.

    Table 8 Comparison of Epoxy/Carbon composite pressure hulls.

    Table 9 Comparison of Epoxy/Glass composite pressure hulls.

    Table 10Comparison of Epoxy/Boron composite pressure hulls.

    Table 11 Effect of lowest mode imperfection on the buckling pressures ratio.

    Fig. 6. First 12 linear buckling modes of the composite egg-shaped pressure hull.

    Fig. 7. Critical buckling modes (a-f) and their corresponding post buckling modes (A-F) for 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and 10 mm lowest mode imperfections.

    7. Sensitivity study

    A sensitivity investigation is conducted for optimized Epoxy/Carbon lay-up [β1a/β2b], to further investigate the effect of fibers’angles and applied hydrostatic pressure (depth) on the material and stability performance of the optimal design of the composite egg-shaped pressure hull. The influence of fibers’ angles on the composites failure criteria and buckling load multiplier is demonstrated in Figs. 12 and 13 and discussed in the following subsections.

    7.1. Influence of fibers’ angles on the composite’s failure criteria

    The influence of fibers’angles β1and β2on the safety factors of composites failure criteria is plotted in Fig.12.Fig.12(a)illustrates that both safety factors for Epoxy/Carbon show similar trends with respect to fiber’s angle (β1). Both safety factors FSTWand FSTHget larger with the increase of fiber’s angle (β1) and highest values of both safety factors occur at 90°value of fiber’s angle(β1).Likewise,it can be noted from Fig. 12 (b) that both FSTWand FSTHfirst decrease and then get larger with the rise of fiber’s angle (β2) and the FSTWand FSTHsafety factors approach their highest values at β2= 24°and β2= 22°respectively. Afterward, both safety factors slightly decrease and then become almost constant with respect to fiber’s angle (β2) until β2= 90°. The influence of fibers’ angles β1and β2on both safety factors for Epoxy/Glass is plotted in Fig.12(c)and(d).It can be observed from Fig.12(c)that FSTWand FSTHfollow almost identical tendencies.Both safety factors gradually get larger with the rise of fiber’s angle (β1) and attain the highest values at β1=90°.Likewise,it is clear from Fig.12(d)that both safety factors also follow nearly identical tendencies and at the start,both safety factors decrease and then get larger with the rise of fiber’s angle(β2) and FSTWand FSTHapproach their highest values at β2= 29°and β2=26°respectively.Afterward,both safety factors marginally decrease and then remain constant with respect to fiber’s angle β2.The influence of the fibers’ β1and β2on both safety factors for Epoxy/Boron is plotted in Fig.12(e)and(f).It can be observed from Fig.12 (e) that both FSTWand FSTHdecline with the rise of fiber’s angle (β1) until β1=19°. Afterward, both safety factors get larger with the rise of fiber’s angle(β1)and approach their highest values at β1=90°. Likewise, both FSTWand FSTHdecline with the rise of fiber’s angle(β2)and approach their smallest values at β2=90°as demonstrated in Fig.12 (f).

    Fig. 8. Critical buckling modes (g-l) and their corresponding post buckling modes (G-L) for 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and 10 mm worst mode imperfections.

    7.2. Influence of fibers’ angles on the structural stability

    The influence of fibers’ angles β1and β2on the buckling load multiplier is demonstrated in Fig.13.It can be observed from Fig.13(a) that buckling load multiplier (λβ1) for Epoxy/Carbon remains almost constant until β1=35°, it then suddenly drops at β1=36°,after which it begins to increase again and suddenly drops again at β1=50°.After this point,it begins to increase again and reaches the highest value at β1=67°. The buckling load multiplier (λβ1) then linearly gets smaller with the rise of fiber’s angle (β1) and approaches its lowest value at β1=90°. Likewise, the buckling load multiplier(λβ2)shows a decreasing tendency with the rise of fiber’s angle (β2) until β2= 46°. Afterward, it begins to get large and approaches its highest value at β2=74°.The buckling load multiplier(λβ2) then quickly gets smaller and approaches its lowest value at β2= 90°. The influence of fibers’ angles β1and β2on the buckling stability of Epoxy/Glass is plotted in Fig. 13 (b). This figure demonstrates that buckling load multiplier(λβ1)gets large linearly with the rise of fiber’s angles(β1)until β1=17°,after which it begins to get smaller until β1=30°and it begins to rise again and approaches the highest value at β1=78°. Afterward, a declining tendency can be noted until β1= 90°. The buckling load multiplier (λβ2) for Epoxy/Glass gets smaller with the rise of the fiber’s angles(β2)and approaches its lowest values at β2= 45°. It then slowly increases with the rise of the fiber’s angle (β2). The influence of the fibers’angles β1and β2on the buckling stability of Epoxy/Boron is plotted in Fig.13(c).It can be noted from this figure that the buckling load multiplier (λβ1) rapidly drops close to the beginning point, after which a slow declining tendency can be noted in its value until β1= 31°. At this position, the buckling load multiplier (λβ1) approaches its lowest value.It then begins to get large again with the rise of the fiber’s angle (β1) and approaches its highest value at β1=86°. Likewise, the value of buckling load multiplier (λβ2) for Epoxy/Boron composite follows a slowly declining tendency until β2= 49°. After this position, it begins to decline quickly with the rise of the fiber’s angle (β2) until β2= 62°, after which its value remains almost constant.

    Fig. 9. Failure modes of Epoxy/Carbon composite pressure hulls after testing [44,64].

    Fig.10. Load displacement graph for 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and 10 mm lowest mode imperfections.

    Fig.11. Load displacement graph for 1 mm, 2 mm, 3 mm, 4 mm, 5 mm and 10 mm worst mode imperfections.

    8. Conclusion

    This paper focuses on the design optimization of the composite egg-shaped submersible pressure hull by means of a geneticalgorithm (GA) and FEA in ANSYS as a first attempt to obtain optimal designs of the composite egg-shaped pressure hull for manufacturing or further investigations. Non-linear buckling analysis of the composite egg-shaped pressure hull is also carried out and its results are compared with that of the composite spherical pressure hull. The following points may be concluded from the study conducted in this paper.

    Table 12 Effect of worst mode imperfection on the buckling pressures ratio.

    a) Genetic algorithms and FEA can be efficiently used to obtain optimal designs of the composite egg-shaped pressure hull.Instead of carrying out parametric design study, the design engineers and researchers,by using the approach adopted in this study, can efficiently obtain optimal designs of the composite pressure hulls for construction or further investigation such as non-linear buckling analysis,impact analysis and underwater explosion analysis.

    b) Epoxy/Carbon composite gives the highest decrease of 54.2%in the buoyancy factor as compared to HY100 steel. Epoxy/Boron demonstrates the highest decrease in the buoyancy factor for 3 out of 5 lay-up arrangements. Likewise, Epoxy/Carbon demonstrates the highest decrease in the buoyancy factor for 2 out of 5 lay-up arrangements.In addition,Epoxy/Glass demonstrates little decrease in the buoyancy factor than Epoxy/Carbon and Epoxy/Boron.

    c) The composite cylindrical pressure hull has higher buoyancy factors as compared to composite spherical and egg-shaped pressure hulls for all lay-ups and materials considered. This means that composite cylindrical pressure hull has more weight as compared to composite spherical and egg-shaped pressure hulls for the same design pressure. The overall performance of the composite egg-shaped pressure hull and composite spherical pressure hull in terms of reduction in weight,material failure and structural stability is almost the same.

    Fig.12. Influence of fibers’ angles β1 and β2 on the safety factors FSTW and FSTH.

    d) The size of imperfection has a very profound effect on the buckling pressure of the composite egg-shaped and spherical submersible pressure hulls. In the case of both lowest and worst modes imperfections, the composite egg-shaped pressure hull demonstrates smaller decrease in buckling pressure as compared to the composite spherical pressure hull. This means that the imperfection sensitivity of the composite egg-shaped pressure hull is lower than that of the composite spherical hull.

    Fig.13. Influence of fibers’ angles β1 and β2 on the buckling load multiplier.

    Declaration of competing interest

    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.

    Acknowledgment

    This work is supported by the National Natural Science Foundation of China research grant #51679056 and Natural Science Foundation of Heilongjiang Province of China grant#E2016024.

    嫁个100分男人电影在线观看| 亚洲伊人色综图| 精品久久蜜臀av无| 91av网站免费观看| 婷婷精品国产亚洲av在线 | 热99国产精品久久久久久7| 嫁个100分男人电影在线观看| tube8黄色片| 99国产精品一区二区三区| 成熟少妇高潮喷水视频| 久久精品91无色码中文字幕| 男女午夜视频在线观看| 激情在线观看视频在线高清 | 色播在线永久视频| 男人操女人黄网站| 韩国精品一区二区三区| 美女 人体艺术 gogo| 俄罗斯特黄特色一大片| 国产在线观看jvid| 又黄又爽又免费观看的视频| 国产成人啪精品午夜网站| 夜夜爽天天搞| 国产蜜桃级精品一区二区三区 | 亚洲成av片中文字幕在线观看| 建设人人有责人人尽责人人享有的| 国产aⅴ精品一区二区三区波| 最近最新免费中文字幕在线| 午夜福利,免费看| 好看av亚洲va欧美ⅴa在| 这个男人来自地球电影免费观看| 大型黄色视频在线免费观看| 日韩欧美在线二视频 | 久久久久视频综合| 亚洲av美国av| 久久国产亚洲av麻豆专区| cao死你这个sao货| 久久人人97超碰香蕉20202| 午夜老司机福利片| 高清视频免费观看一区二区| 午夜日韩欧美国产| 国产成人精品无人区| 天天躁日日躁夜夜躁夜夜| 天堂√8在线中文| 成年版毛片免费区| 欧美日韩福利视频一区二区| 亚洲精品在线观看二区| www.熟女人妻精品国产| 深夜精品福利| 丝袜美足系列| 日本黄色日本黄色录像| 国产成人免费无遮挡视频| 亚洲七黄色美女视频| 天天躁夜夜躁狠狠躁躁| 久久草成人影院| 国产一区二区激情短视频| 亚洲精品一卡2卡三卡4卡5卡| 99久久综合精品五月天人人| 日韩熟女老妇一区二区性免费视频| 成年人午夜在线观看视频| 国产精华一区二区三区| 一进一出好大好爽视频| 咕卡用的链子| 久99久视频精品免费| 啦啦啦 在线观看视频| 脱女人内裤的视频| a级毛片在线看网站| 日韩成人在线观看一区二区三区| 日日爽夜夜爽网站| 色精品久久人妻99蜜桃| 99re在线观看精品视频| 热99国产精品久久久久久7| 中文字幕色久视频| 日韩熟女老妇一区二区性免费视频| 国产精品免费视频内射| 狠狠婷婷综合久久久久久88av| 欧洲精品卡2卡3卡4卡5卡区| 精品一区二区三区视频在线观看免费 | 中文字幕人妻丝袜一区二区| 午夜亚洲福利在线播放| 99热网站在线观看| 黑人巨大精品欧美一区二区蜜桃| 无遮挡黄片免费观看| 午夜精品国产一区二区电影| 可以免费在线观看a视频的电影网站| 亚洲精品在线观看二区| 亚洲欧美日韩高清在线视频| av超薄肉色丝袜交足视频| 午夜影院日韩av| www.熟女人妻精品国产| 日韩熟女老妇一区二区性免费视频| 黑人巨大精品欧美一区二区mp4| 亚洲综合色网址| 日韩视频一区二区在线观看| 一级a爱视频在线免费观看| 国产亚洲欧美98| 亚洲第一青青草原| 欧美日韩国产mv在线观看视频| 国产极品粉嫩免费观看在线| 91精品三级在线观看| 热re99久久精品国产66热6| 黄色丝袜av网址大全| 久久中文字幕一级| 新久久久久国产一级毛片| 人人澡人人妻人| 欧美日韩视频精品一区| 国产亚洲欧美在线一区二区| 亚洲 国产 在线| 国产成人精品久久二区二区免费| 一级作爱视频免费观看| 国产欧美日韩精品亚洲av| 建设人人有责人人尽责人人享有的| 9色porny在线观看| 免费观看a级毛片全部| 成年人免费黄色播放视频| 国产精品98久久久久久宅男小说| www.精华液| 亚洲av美国av| 黄色 视频免费看| 色尼玛亚洲综合影院| 俄罗斯特黄特色一大片| 亚洲成人免费电影在线观看| 日本wwww免费看| 亚洲熟妇中文字幕五十中出 | 亚洲欧美一区二区三区黑人| 九色亚洲精品在线播放| 欧美日韩国产mv在线观看视频| 日本a在线网址| 美女高潮到喷水免费观看| tube8黄色片| 激情在线观看视频在线高清 | 亚洲自偷自拍图片 自拍| 国产精品免费大片| 欧美乱妇无乱码| 精品一区二区三卡| 91九色精品人成在线观看| 亚洲人成电影免费在线| 90打野战视频偷拍视频| 日韩成人在线观看一区二区三区| 色老头精品视频在线观看| 色尼玛亚洲综合影院| 制服诱惑二区| 99re6热这里在线精品视频| 亚洲色图综合在线观看| 成人免费观看视频高清| 欧美日韩瑟瑟在线播放| 欧美日韩视频精品一区| 欧美日韩精品网址| 亚洲一码二码三码区别大吗| 天天躁日日躁夜夜躁夜夜| 成在线人永久免费视频| 91九色精品人成在线观看| 亚洲精品中文字幕一二三四区| www.自偷自拍.com| 中文字幕色久视频| 国产精品成人在线| 久久午夜综合久久蜜桃| 黄色视频,在线免费观看| 99国产极品粉嫩在线观看| av不卡在线播放| 777米奇影视久久| 久久久久久久久久久久大奶| 国产免费av片在线观看野外av| 亚洲精品一卡2卡三卡4卡5卡| 黄色怎么调成土黄色| 久久久精品国产亚洲av高清涩受| av网站免费在线观看视频| 精品人妻熟女毛片av久久网站| 12—13女人毛片做爰片一| 欧美乱色亚洲激情| 美女 人体艺术 gogo| 亚洲七黄色美女视频| 免费在线观看完整版高清| 一进一出抽搐gif免费好疼 | 国产xxxxx性猛交| 国产精品一区二区在线不卡| 中文欧美无线码| 大型av网站在线播放| 亚洲精品粉嫩美女一区| 一区二区日韩欧美中文字幕| 亚洲少妇的诱惑av| 亚洲精品一卡2卡三卡4卡5卡| 亚洲情色 制服丝袜| 日本一区二区免费在线视频| 久久天躁狠狠躁夜夜2o2o| 色94色欧美一区二区| 日韩欧美一区视频在线观看| 老司机靠b影院| 久热爱精品视频在线9| 日本wwww免费看| 18在线观看网站| 美女高潮喷水抽搐中文字幕| 男人舔女人的私密视频| 侵犯人妻中文字幕一二三四区| 黄色怎么调成土黄色| 亚洲一区二区三区欧美精品| 亚洲成av片中文字幕在线观看| 乱人伦中国视频| 亚洲欧美激情在线| 啦啦啦在线免费观看视频4| aaaaa片日本免费| 国产三级黄色录像| 国产高清激情床上av| 亚洲精品在线美女| 99国产精品一区二区三区| 啦啦啦在线免费观看视频4| 欧美精品av麻豆av| 夜夜夜夜夜久久久久| a在线观看视频网站| 大陆偷拍与自拍| 欧美午夜高清在线| 午夜影院日韩av| 亚洲精品国产区一区二| 日韩免费高清中文字幕av| 91九色精品人成在线观看| 人妻久久中文字幕网| 久久久久国内视频| 欧美日韩乱码在线| 91精品三级在线观看| 国产xxxxx性猛交| 国产在线一区二区三区精| 亚洲精品中文字幕在线视频| 欧美性长视频在线观看| 欧美日韩国产mv在线观看视频| 精品卡一卡二卡四卡免费| 91老司机精品| 最新美女视频免费是黄的| 久久久久国产一级毛片高清牌| 99热国产这里只有精品6| 母亲3免费完整高清在线观看| 久久人妻福利社区极品人妻图片| 十八禁网站免费在线| 亚洲视频免费观看视频| 久久久精品免费免费高清| 日韩制服丝袜自拍偷拍| 在线观看66精品国产| www.熟女人妻精品国产| 波多野结衣一区麻豆| 狠狠婷婷综合久久久久久88av| 亚洲欧美一区二区三区黑人| 免费人成视频x8x8入口观看| 人妻一区二区av| 少妇被粗大的猛进出69影院| 国产精品久久久久成人av| 亚洲人成电影免费在线| 亚洲色图 男人天堂 中文字幕| 亚洲精品自拍成人| 亚洲av美国av| 老司机影院毛片| 一区二区日韩欧美中文字幕| 精品国产一区二区三区四区第35| 美女午夜性视频免费| 精品久久久久久,| 久久人人97超碰香蕉20202| 亚洲中文字幕日韩| 午夜福利影视在线免费观看| 欧美乱码精品一区二区三区| 50天的宝宝边吃奶边哭怎么回事| 91大片在线观看| 欧美大码av| 久久久国产一区二区| 777米奇影视久久| 大片电影免费在线观看免费| 色婷婷av一区二区三区视频| av片东京热男人的天堂| 国产欧美日韩一区二区精品| 69精品国产乱码久久久| 一级片免费观看大全| 亚洲一码二码三码区别大吗| 久久久久久久久久久久大奶| 久久精品国产亚洲av高清一级| 国产视频一区二区在线看| 又黄又爽又免费观看的视频| 好看av亚洲va欧美ⅴa在| 日韩精品免费视频一区二区三区| 一级黄色大片毛片| 欧美国产精品va在线观看不卡| 色老头精品视频在线观看| 欧美不卡视频在线免费观看 | 国产精品 国内视频| videosex国产| 国产av精品麻豆| 久久人妻熟女aⅴ| 大型av网站在线播放| 久久中文字幕一级| 欧美在线一区亚洲| 母亲3免费完整高清在线观看| 大型av网站在线播放| 在线看a的网站| 俄罗斯特黄特色一大片| 一夜夜www| 国产一区二区激情短视频| 无遮挡黄片免费观看| 国产精品香港三级国产av潘金莲| 黄色女人牲交| www.自偷自拍.com| 十八禁高潮呻吟视频| 国产一区有黄有色的免费视频| 国产精品久久久av美女十八| 久久久精品免费免费高清| 日韩欧美三级三区| av有码第一页| 亚洲九九香蕉| 狠狠婷婷综合久久久久久88av| 国产黄色免费在线视频| 国产精品美女特级片免费视频播放器 | 亚洲一区高清亚洲精品| 国产高清视频在线播放一区| 午夜免费成人在线视频| 亚洲国产看品久久| 大型黄色视频在线免费观看| av有码第一页| 国产欧美日韩一区二区精品| 黑人欧美特级aaaaaa片| 精品亚洲成国产av| 老熟妇仑乱视频hdxx| 国产精品香港三级国产av潘金莲| 国产欧美亚洲国产| 亚洲熟女精品中文字幕| svipshipincom国产片| 一区二区三区国产精品乱码| 十八禁高潮呻吟视频| 大香蕉久久网| a在线观看视频网站| 欧美另类亚洲清纯唯美| 搡老熟女国产l中国老女人| 纯流量卡能插随身wifi吗| 亚洲av美国av| 日韩人妻精品一区2区三区| 日韩欧美在线二视频 | 日韩精品免费视频一区二区三区| 天天添夜夜摸| 亚洲成国产人片在线观看| 久久午夜综合久久蜜桃| 久热爱精品视频在线9| av天堂在线播放| 国产精品一区二区在线观看99| 久久人妻熟女aⅴ| 亚洲综合色网址| 九色亚洲精品在线播放| 午夜福利欧美成人| 国产成人精品久久二区二区91| 精品一区二区三区四区五区乱码| www日本在线高清视频| 一区二区日韩欧美中文字幕| 亚洲av欧美aⅴ国产| 国产人伦9x9x在线观看| 欧美日韩av久久| 免费黄频网站在线观看国产| 在线看a的网站| 多毛熟女@视频| 欧美成人午夜精品| 久久中文看片网| 少妇裸体淫交视频免费看高清 | 1024视频免费在线观看| 中文字幕色久视频| 女人爽到高潮嗷嗷叫在线视频| 男人舔女人的私密视频| 满18在线观看网站| 国产精品 欧美亚洲| 欧美日韩精品网址| 他把我摸到了高潮在线观看| 久久亚洲真实| 国产精品一区二区在线不卡| 男人舔女人的私密视频| 一个人免费在线观看的高清视频| 99久久精品国产亚洲精品| 丝袜美足系列| 亚洲伊人色综图| 老鸭窝网址在线观看| 久久午夜综合久久蜜桃| 亚洲三区欧美一区| 精品一品国产午夜福利视频| 国产成人精品无人区| 亚洲在线自拍视频| 黄网站色视频无遮挡免费观看| 欧美黑人欧美精品刺激| 成人黄色视频免费在线看| 国产成人精品久久二区二区免费| 老熟妇仑乱视频hdxx| а√天堂www在线а√下载 | av线在线观看网站| 人人妻人人澡人人爽人人夜夜| e午夜精品久久久久久久| 欧美 亚洲 国产 日韩一| 色在线成人网| 精品人妻熟女毛片av久久网站| 另类亚洲欧美激情| 亚洲精品国产精品久久久不卡| 嫁个100分男人电影在线观看| svipshipincom国产片| 一边摸一边做爽爽视频免费| 国产精品国产av在线观看| 国产av又大| 免费在线观看亚洲国产| 最新美女视频免费是黄的| 国产成人欧美| 亚洲国产精品sss在线观看 | 美女午夜性视频免费| 久久精品熟女亚洲av麻豆精品| 久久久久久久精品吃奶| 少妇猛男粗大的猛烈进出视频| 久久精品aⅴ一区二区三区四区| 色综合欧美亚洲国产小说| 久久久国产欧美日韩av| 一级作爱视频免费观看| 亚洲国产欧美日韩在线播放| 丝瓜视频免费看黄片| 99国产精品免费福利视频| 亚洲精品久久成人aⅴ小说| 不卡av一区二区三区| 黄片小视频在线播放| 国产精华一区二区三区| 国产成人欧美在线观看 | 亚洲 国产 在线| 欧美黑人精品巨大| 女人久久www免费人成看片| 精品熟女少妇八av免费久了| 亚洲av片天天在线观看| 欧美av亚洲av综合av国产av| 国产成人影院久久av| 欧美久久黑人一区二区| 久久婷婷成人综合色麻豆| 久久精品国产亚洲av高清一级| 99久久综合精品五月天人人| 欧美乱色亚洲激情| 精品亚洲成a人片在线观看| 亚洲国产精品一区二区三区在线| 成人三级做爰电影| 一二三四在线观看免费中文在| 无遮挡黄片免费观看| 最近最新中文字幕大全电影3 | 欧美精品高潮呻吟av久久| 亚洲久久久国产精品| www.999成人在线观看| 18禁裸乳无遮挡免费网站照片 | 成年版毛片免费区| 叶爱在线成人免费视频播放| 亚洲va日本ⅴa欧美va伊人久久| 国产三级黄色录像| av中文乱码字幕在线| 久久人妻熟女aⅴ| 精品一区二区三区av网在线观看| 极品少妇高潮喷水抽搐| 91国产中文字幕| 亚洲va日本ⅴa欧美va伊人久久| 国产91精品成人一区二区三区| 男人操女人黄网站| 色综合婷婷激情| 国产淫语在线视频| 在线av久久热| 少妇粗大呻吟视频| 国产精品一区二区在线观看99| 怎么达到女性高潮| 国产一区二区三区视频了| 中文字幕人妻丝袜制服| xxxhd国产人妻xxx| av视频免费观看在线观看| 18禁裸乳无遮挡动漫免费视频| 在线观看www视频免费| 国内毛片毛片毛片毛片毛片| 淫妇啪啪啪对白视频| 欧美精品高潮呻吟av久久| 麻豆国产av国片精品| 亚洲自偷自拍图片 自拍| 亚洲视频免费观看视频| 男女高潮啪啪啪动态图| 亚洲国产精品sss在线观看 | 亚洲全国av大片| 亚洲熟女精品中文字幕| 在线十欧美十亚洲十日本专区| 高清视频免费观看一区二区| 日本wwww免费看| 精品久久久精品久久久| 在线永久观看黄色视频| 很黄的视频免费| 手机成人av网站| 欧美老熟妇乱子伦牲交| 国产成人一区二区三区免费视频网站| 美女视频免费永久观看网站| 亚洲五月天丁香| 国产99久久九九免费精品| 捣出白浆h1v1| 正在播放国产对白刺激| 男女午夜视频在线观看| 久久精品亚洲av国产电影网| 国产精品久久久久成人av| 日韩欧美一区二区三区在线观看 | 一二三四社区在线视频社区8| 久久久国产成人精品二区 | 久久午夜综合久久蜜桃| 无遮挡黄片免费观看| 欧美日韩视频精品一区| av国产精品久久久久影院| 亚洲男人天堂网一区| 大码成人一级视频| 中文字幕制服av| 国产精品成人在线| 在线观看66精品国产| 18禁美女被吸乳视频| 中文字幕av电影在线播放| 国产主播在线观看一区二区| 国产男女内射视频| 一级a爱片免费观看的视频| 午夜老司机福利片| 9热在线视频观看99| 老司机影院毛片| 久久久久久久久免费视频了| 亚洲成人免费电影在线观看| 亚洲国产欧美一区二区综合| 激情在线观看视频在线高清 | 一区二区三区国产精品乱码| 国产一区二区三区综合在线观看| xxx96com| 在线永久观看黄色视频| 国产精品欧美亚洲77777| 久热这里只有精品99| 精品午夜福利视频在线观看一区| 免费黄频网站在线观看国产| 操出白浆在线播放| 日韩中文字幕欧美一区二区| 亚洲黑人精品在线| 亚洲一区中文字幕在线| а√天堂www在线а√下载 | 久久精品91无色码中文字幕| 国产精品一区二区精品视频观看| 少妇裸体淫交视频免费看高清 | 窝窝影院91人妻| 99re6热这里在线精品视频| 国产亚洲欧美精品永久| 好男人电影高清在线观看| 免费看十八禁软件| 免费在线观看视频国产中文字幕亚洲| 99国产极品粉嫩在线观看| 久99久视频精品免费| 麻豆成人av在线观看| 亚洲一码二码三码区别大吗| 丁香六月欧美| 日韩欧美三级三区| 亚洲精品国产色婷婷电影| 巨乳人妻的诱惑在线观看| 亚洲五月婷婷丁香| 国产精品久久久久久精品古装| 日韩大码丰满熟妇| 国产成人av激情在线播放| 宅男免费午夜| 欧美国产精品一级二级三级| 亚洲少妇的诱惑av| 亚洲人成电影免费在线| 777久久人妻少妇嫩草av网站| 久久精品成人免费网站| 一区二区三区国产精品乱码| 91精品三级在线观看| 亚洲色图综合在线观看| 欧美黄色片欧美黄色片| www.精华液| 人人妻,人人澡人人爽秒播| 两性夫妻黄色片| 悠悠久久av| 两性夫妻黄色片| 午夜精品国产一区二区电影| 91国产中文字幕| 欧美成人午夜精品| 亚洲一区二区三区欧美精品| 天堂俺去俺来也www色官网| 成人手机av| 亚洲专区国产一区二区| 国产精品一区二区在线观看99| 校园春色视频在线观看| 欧美 日韩 精品 国产| 麻豆av在线久日| 久久国产精品影院| 亚洲欧美一区二区三区黑人| 色婷婷久久久亚洲欧美| 亚洲av熟女| 成人永久免费在线观看视频| 在线天堂中文资源库| 美女午夜性视频免费| 久久香蕉激情| 精品少妇一区二区三区视频日本电影| 亚洲精品美女久久久久99蜜臀| 91麻豆av在线| 亚洲精品久久成人aⅴ小说| 亚洲av日韩在线播放| 搡老熟女国产l中国老女人| 久久亚洲真实| 免费一级毛片在线播放高清视频 | 啪啪无遮挡十八禁网站| 好看av亚洲va欧美ⅴa在| 日韩欧美免费精品| 高清欧美精品videossex| 亚洲aⅴ乱码一区二区在线播放 | 激情在线观看视频在线高清 | 亚洲色图综合在线观看| 国产精品国产av在线观看| 50天的宝宝边吃奶边哭怎么回事| 捣出白浆h1v1| 精品少妇久久久久久888优播| 日韩熟女老妇一区二区性免费视频| 在线观看舔阴道视频| 满18在线观看网站| av电影中文网址| 亚洲国产精品一区二区三区在线| 麻豆乱淫一区二区| 久久香蕉激情| 一进一出抽搐gif免费好疼 | 一个人免费在线观看的高清视频| 久久影院123| 十八禁网站免费在线| 91在线观看av| 中文字幕人妻丝袜制服| 国产精品乱码一区二三区的特点 | 亚洲在线自拍视频| 国产三级黄色录像| 丝袜美足系列| 亚洲中文av在线|