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

    Optimal contract wall for desired orientation of fibers and its effect on flow behavior*

    2017-06-07 08:22:46WeiYang楊煒
    水動力學研究與進展 B輯 2017年3期

    Wei Yang (楊煒)

    College of Mechanical Engineering, Ningbo University of Technology, Ningbo 315010, China, E-mail: wyang3715@sina.com

    Optimal contract wall for desired orientation of fibers and its effect on flow behavior*

    Wei Yang (楊煒)

    College of Mechanical Engineering, Ningbo University of Technology, Ningbo 315010, China, E-mail: wyang3715@sina.com

    2017,29(3):495-503

    The orientation of suspended fibers in the turbulent contraction is strongly related to the contraction ratio, which in some cases may be detrimental to the actual production. Here for a certain contraction ratio, the contraction geometry shape is optimized to obtain the desired fiber orientation. In view of the nonlinearity and the complexity of the turbulent flow equations, the parameterized shape curve, the dynamic mesh and a quasi-static assumption are used to model the contraction with the variable boundary and to search the optimal solution. Furthermore the Reynolds stress model and the fiber orientation distribution function are solved for various wall shapes. The fiber orientation alignment at the outlet is taken as the optimization objective. Finally the effect of the wall shape on the flow mechanism is discussed in detail.

    Fibers orientation, wall shape, turbulent contraction, dynamic mesh, nonlinear optimization

    Introduction

    The fiber filling can effectively improve the stiffness and the strength of manufactured goods, whose performances depend largely on the fiber orientation distribution. The contraction is one of common internal flows, mainly used for the speedup flow, such as the pulp flow box in the paper industry and the jet nozzle of the extrusion die. It is shown that the flow behaviors inside the contraction are more complex than in the channel due to the acceleration effect, which leads to a higher alignment of the fiber orientation under the same entrance conditions. And what’s more, the fiber orientation might very likely to follow the flow direction for a large contraction ratio C because the strong streamwise strain rate, rather than the stochastic turbulence effect, would play a major rolein the vector field at this moment. As the fiber orientation alignment may affect the overall properties of the industrial products significantly, it is very important to investigate how to control the fiber orientation under a given large contraction ratio.

    The dynamics of the fiber suspension and the orientation behavior of the fibers suspended in the planar contraction were much studied. Batchelor and Proudman[1]and Ribner and Tucker[2]considered the cases when C is very large and the time scale of the flow is much smaller than that of the vortex, and suggested that the development of the turbulent kinetic energy is dependent on the local flow strain rate which is related to the geometry of the flow field. Hussain and Ramjee[3]concluded that the flow shape and the Reynolds number in the axial symmetric contraction flow have little effect on the turbulent characteristics, but the inlet condition is the most important parameter. A series of studies by Lin et al.[4-8]show that with the increase of the Reynolds number more fibers will be aligned with the flow direction in the laminar flow regime while the fiber orientation distributions become more homogeneous in the turbulent regime. For the fibers with a high aspect ratio the mean velocity distribution has a more significant effect on the fiber orientation distribution than the fluctuating velocity. Olsonet al.[9,10]calculated the fiber orientation distribution with a constant diffusion coefficient in the contraction and it is shown that the influence of the contraction ratio on the fiber orientation is more significant than the entrance velocity. Parsheh et al.[11,12]employed a model with the rotational diffusion rate depending heavily on the turbulence of the entrance and decaying exponentially with the contraction ratio. Then four contraction wall shapes, including the constant mean rate of strain, the linear rate of strain, the quadratic rate of strain and the flat walls, were considered to analyze the influence of the contraction geometry on the orientation anisotropy. It is shown that the anisotropic orientation of fibers is dependent on the contraction wall shape while both the turbulence intensity and the rotational diffusion coefficient have little to do with the contraction shape. And what’s more, the rotational Peclet number, defined as the streamwise velocity gradient divided by the rotational diffusion coefficient, is shown to be the crucial quantity to determine whether the turbulent effect on the fiber orientation is dominant or not.

    Although the optimal shape is not given finally by Parsheh et al.[11], the importance of the flow morphology for the motion and orientation of fibers is shown by the fact that the contraction with flat walls has the smallest orientation anisotropy and the contraction with the constant rate of strain has the largest anisotropy at the outlet . Thus the optimization of the wall shape is the way to deal with the anisotropic distribution of fibers under a constant contraction ratio, just as did by M?kinen and H?m?l?inen[13], Nazemi and Farahi[14]and Farhadinia[15]. M?kenin and H?m?l?inen[13]discretized the contaction shape with the Bezier function and solved the diffusion–convection equation with the streamline upwind Petrov-Galerkin (SUPG) finite element method. Nazemi and Farahi[14]and Farhadinia[15]turned the shape optimization problem into the coefficient optimization of the state equation, which is related to the average flow velocity gradient and the fiber orientation angle. By using an embedding method, the optimal measure representing the optimal shape is approximated by the solution of a linear programming problem. But the assumption of the quasi one dimensional potential flow was made and the effect of turbulence was neglected in their optimization solutions.

    For the nonlinear and multi-variable coupling trait of the turbulence control equation, the latest optimization methods for the wall shape are mainly divided into two major categories. One is the reverse design method with the design variables obtained to meet a given objective function value by solving the control equations or by a progressive approach. But this method brings about many technical problems in practical applications: (1) It is difficult to obtain the relationship directly between the variables and the objective function, because a variable is often coupled with other variables in a group of turbulence dynamics equations. (2) The solution of the turbulent dynamics equation would be very sensitive to the mutative flow state or the boundary. Therefore, the inverse design method faces problems of sensitivity and uncertainty of the solution. (3) The method requires some relevant experience or parameters to initialize, which would make the solutions fall into a small range and make the whole optimization system lack of robustness.

    The other is the optimal design method. The optimal solution with respect to the objective function can be obtained by combining the direct solution with the optimization searching of the solution space with a good performance in the practice[16,17]. But this method often uses a deputy model or an approximate calculation model, associated with the number of parameters, the sample size, the variable range and the approximate precision, and the method still needs to be adjusted and tested repeatedly for the optimization accuracy and reliability.

    In order to achieve the accuracy of the dynamics solution, the deputy model is not used here. The dynamic mesh technology is adopted to model the variable flow boundary and to obtain directly the solution of the fiber orientation distribution in the turbulence. Then a parametric curve equation of the wall shape as Fig.1 is constructed at first. Based on the dynamic mesh model and the quasi-static assumption, the optimal parameters of the wall geometry can be obtained for the desired fiber orientation distribution. The effects of the wall shape parameters of the contraction, such as the inlet tangential angle a, the outlet tangential angle b and the ratioLR of the length L to the height0h, on the flow behavior are studied in detail.

    Fig.1 The 2-D model of plannar contraction

    1. Parametric optimal model

    1.1 Parametric equation of wall shape

    Based on the bicubic parametric curves, the quintic curve, and the Vyshinsky curve for the contraction section of the wind tunnel, the parametric cubic curve for the wall shape function is constructed as follows

    where0h andeh are the heights of the inlet and the outlet, respectively, and L is the length of the contraction, as shown in Fig.1. The parameters a, b, c and d can be determined to meet the boundary conditions and its first order derivative conditions:

    where C is the contraction ratio of h0to heand RL=L/h0.

    The coefficients a, b and c in the Eq.(2) depend on the parameters a, b andLR of the boundary conditions, so various shapes of the wall surface can be obtained by adjusting these three parameters.

    1.2 Optimization model

    The probability distribution f(x,j) of the fiber orientation at the center streamline of the contraction can be calculated as[18].

    where A is the positive definite matrix, =[0,1]A ,hb andhc are the coefficient matrices, j is the orientation angle of the fibers, andfpD is the orientation diffusion coefficient of the fibers.

    Optimization objective function is as follows:

    whereiX is the parameter variable to be optimized, such as a, b and c in Eq.(2).1111a is the anisotropic index of the fiber orientation, defined as1111=aHereip is the thi component of p =[cosj, sin j]T.

    1.3 Dynamic mesh model

    The body-fitted grid method, the arbitrary Lagrange-Euler method and the immersed boundary method are the common methods for the flow with the variable boundary, and a large number of flow problems are solved by these methods[19-21], such as those of the complex surface fluid, the dynamic boundary, the free or elastic surface and the flow field interface interaction. It is noted that these three methods are generally applied to the simulations of the low Reynolds number movement. For the turbulence problems with a high Reynolds number, these methods do not perform well due to the complexity of the coordinate transformation or the force source term for the turbulence model. In this paper the dynamic mesh technique is adopted to solve the RMS equation for the small deformation boundary.

    1.4 Spring-based smoothing method

    The spring-based smoothing method is utilized to realize the deformation of the dynamic mesh. At first the edges between any two mesh nodes are idealized as a network of interconnected springs, as shown in Fig.2. Hereijk is the spring constant between node i and its neighbor j,

    Fig.2 Spring-based smoothing algorithm model

    Based on the Hook?s Law and the force-balance principle, an iterative equation is applied as Eq.(5). The positions of the mesh nodes are updated as Eq.(6) until convergence.

    where x is the relaxation factor of the mesh nodes.x =0 means that the nodes in the mesh are not movable.

    1.5 Remeshing method

    When the node displacement is larger than thelocal cell sizes, the cell quality might deteriorate and the mesh distortion will even generate negative cell volumes, which would cause convergence problems. The local remeshing method is required to circumvent the problem.

    There are two criteria for the local remeshing. One is the cell distortion rate, the other is the cell size. Once the cell distortion rate is greater than a specified maximum value or the cell size exceeds the range between the specified minimum length scale and the maximum length scale these bad cells would be marked and be locally updated with the new cells.

    1.6 Quasi-static assumption

    The isotropic orientation distribution of the fibers is expected at the outlet of the contraction by optimizing the wall shape, which is essentially to obtain the optimal steady state solution of the fiber orientation distribution under the boundary conditions of different geometric shapes. In fact the boundary is not in motion during the search for the solution, so it is necessary to introduce the quasi-static hypothesis. According to the dynamic mesh and the conservation equations, the parameter m is used to control the boundary displacement to eliminate the influence of the non steady state term for the cell motion on the solution of the conservation equation. When the flow boundary changes from one curve to the other curve, it should undergo the process of a slow variation m times to reach the steady state, to meet the quasi-static condition.

    Fig.31111a under different parameters m

    Figure 3 shows the solution of the fiber orientation distribution1111a along the center streamline with the steady state model and the dynamic mesh model of m =20, m =1000, m =10000. It can be seen obviously that when =20m ,1111a oscillates up and down around its steady-state solution, especially in the front pathway with a great shape change. And when m takes a larger value, such as =m 1 000 and =m 10 000, the solution of1111a is basically consistent with the steady-state solution. That is to say, the effect of the boundary motion on the solution can be eliminated. In view of the iteration efficiency, we take =m 1 000.

    2. Optimal solution and discussion

    2.1 Impact of wall shape parameters on1111a

    In order to understand the impact of the parameters of the wall shape on the fiber orientation distribution a1111in the contract turbulence, three parameters take three levels, as shown in Table 1. According to the method of the orth4ogonal experiment design, an orthogonal table L9(3) with three factors and three levels is constructed as in Table 2. The fiber orientation distribution a1111is calculated using Eq.(3) for the nine cases, some results are shown in Fig.4.

    Table 1 The levels of parameters for wall shape

    Table 2 L9(34) of parameters for wall shape

    It can be seen clearly from Fig.4 that the profile of a1111changes correspondingly with the shape of the wall. The greater the b is, the larger the a1111will be. As known in Table 2,a1111,out- 1/p reaches the minimum value when a= -8o, b= 0oand RL=5.02. The influence order of the parameters on the optimal object can be obtained as in Table 3 by using the range analysis method.

    Fig.4 The distribution of1111a for different walls in L9(34)

    As shown in Table 3, the tangential angle b at the outlet is the leading factor to affect the optimal object a1111,out- 1/p, the next is the field length L and the tangential angle a at the inlet. The fiber orientation distribution at the outlet is the best when b= 0o. Here b is made identically equal to zero to simplify the number of parameters. As a result, the number of the optimal variables changes from three to two, i.e., the field length L and the inlet tangential angle a. Thus the variable range can be adjusted to 350mm£L£ 650mm , a would be in the appropriate range of [-4 0o,0o] according to the channel length L.

    Table 3 Analysis results ofa1 111,out - 1/p

    2.2 Optimal shape parameters

    According to the scope of the above parameters, various wall shapes are generated by using the dynamic mesh method and the fiber orientation distribution a1111is calculated by Eq.(3) using the inlet data in Table 4, which is the same as in Parsheh’s case. The distribution map of a1111is listed only in the field of L =350mm and -4 0o<a<0oas follows.

    Table 4 Inlet parameters of the contract flow

    It is known that from the trend of1111a in Fig.5 the largest contraction occurs at the inlet at first and then moves forward gradually to a place near the outlet when a changes from the maximum negative angle to zero and the field length L remains constant. Accordingly, the maximum alignment index1111a of the fiber orientation distribution also varies with the contraction trend of the wall shape. Although the flow length L and the tangent angle a at the inlet take different values in each numerical case, the maximum value of a1111no longer occurs at the outlet but in the front of the outlet without exception due to the tangent angle b= 0oat the outlet. And what’s more, a1111values are all close to the average probability distribution, i.e., 1/p, at the outlet which indicates that1111a is effectively controlled in the contraction. Of course in addition to the main effect of b on the fiber orientation, the parameters L and a still have a subtle influence on the maximum value a1111,maxof the fiber orientation distribution and the average value a1111,outat the outlet.

    The values of a1111,maxand a1111,out- 1/p are shown in Fig.6 and Fig.7 for the different cases. a1111,maxcan take the minimum when the parameters L =550mmand a= -1 0o.a1111,out- 1/p can take the minimum when the parameters L =350mm anda =0o. Thus it can be seen that when the wall shape is close to the linear contraction line, the fiber orientation a1111,maxwould be relatively smaller, that is to say, the wall shape near the linear contraction can reduce the fiber alignment in the contract flow. From Fig.7a1111,out- 1/p all can reach the lowest point at the outlet in the cases of a= 0o, which suggests that the fi ber orienta tio n cou ld exh ibit a better isotropy . In ordertoobtainagoodresultof a1111,out- 1/p,thelength of the flow does not seem to need too long, and a slightly shorter length seems more appropriate. If only for the purpose of a1111,out- 1/pmin, the wall curve with the parameters of L =350mm , a=0oand b= 0ois optimal.

    Fig.5 The distribution of a1 111for L =350mm and different a

    Fig.6 The variation of a1 111,maxfor different a and L

    Fig.7 The variation of a1 111,out - 1/p for different a and L

    Fig.8 Turbulent intensity along the center streamline for different a and L=550mm

    2.3 Flow mechanisms and shape parameters

    2.3.1 Impact of a on flow behavior

    It is easy to see from Eq.(3) that the fiber orientation distribution is mainly controlled by the evolution of the velocity gradient and the flow structure in the flow field. The variation of the boundary shapecan change the velocity gradient directly, and a large velocity gradient is one of the main factors to generate a large scale turbulence vortex. The curve of the turbulent structure along the center streamline is obtained by numerical simulations, as shown in Fig.8, when the channel length L =550mm and the inlet tangential angle a takes various values from a= -3 0otoa= 0o.

    The turbulent intensity is the ratio of the mean square root of the turbulent fluctuation velocity to the mean velocity of the flow field. The turbulence strength and the diffusion coefficient of the turbulent transport is generally related to this index. Usually the turbulence in the range of 5%-20% is called the higher turbulence, and that in the range of 1%-5% is called the moderate turbulence; that in the range less than 1% is called the low turbulence. The natural transition, as one of the most common form of the transition, occurs in the case of the low turbulence.

    From Fig.8, the turbulence intensity (T I) in general decreases at first and then increases along the center streamline. When the inlet tangential angle changes from a= -3 0oto a=0othe inflection point of different curves is more and more close to the outlet. When a= -3 0othe wall shape is transformed from the necking mouth to a tiny flaring mouth at the position of x=350mm-400mm , the turbulence intensity begins to oscillate and increase rapidly. When a approaches 0o, the curve changes gradually from oscillation to a stable state. The simulation results show that too large turbulence intensity is not appropriate for the fiber orientation distribution because it may produce a large lateral velocity and affect the flow quality.

    Fig.9 Turbulent vorticity along the center streamline for different a and L=550mm

    The turbulent vorticity ( )Vot is defined as the curl of the fluid velocity vector to describe usually the strength and the direction of the turbulent eddies. It is one of the most important physical quantity of the vortex motion. Figure 9 shows that the change of the vorticity with a is similar to the turbulence intensity, except that the vorticity under a linear wall is close to zero along the center streamline. It can be inferred that the tensile effect under the linear contract wall is stronger than the others, which make the vortex stretch, the vorticity be weakened and the fiber more easily align along the streamline. As a result, the phenomenon of the weak vorticity is common in the flow with a linear contract wall instead of the flow with a wall shape of large curvature variations. And due to the fact thatb=0oat the outlet, the curvature variation of the wall shape is small and the vorticity remains at the same energy level or lower on the whole.

    2.3.2 Impact of L on flow behavior

    Consider the inlet tangential angle a= 0oand the channel length L=350mm-650mm in the numerical simulation. The turbulence changes with the contraction ratio C along the center streamline as shown in Fig.10. Figure 10 shows that the shorter the channel length L is, the more quickly the turbulence intensity of the flow increases with the contraction ratio C at the same inlet tangential angle a, which maybe explains why the fiber orientation is inclined to the state of isotropy in a short contraction.

    Fig.10 Turbulent intensity changes with C along the center streamline

    Fig.11 Turbulent vorticity changes with C along the center streamline

    The maximum amplitude of the turbulent vorticity and its peak position vary with the channel length L, as shown in Fig.11. The maximum peak in the flow field of L =400mm occurs at C =9, its amplitude is about 8 s?1, then, the maximum peak of L= 500 mm occurs almost at C =7.5 and its amplitude is about 6.5 s?1. The maximum peaks of L =450mm, L =600mm and L =650mm occur near the place C =11, its amplitude is about 5 s?1. The maximum peak of L =350mmoccurs near the place C =8.5, its amplitude is also about 5 s?1. The maximum peak of L =550mm occurs near the place C =11, its amplitude is also about 4.5 s?1. Finally the maximum peak in the linear contraction of L =550mm occurs at the outlet, its amplitude is only about 2 s?1. Although from Fig.11 it is not sure about the certain relationship between the channel length L and the peak magnitude of the vorticity, it can be concluded that the maximum amplitude of the turbulence vorticity very likely occurs slightly away from the outlet in order to reduce the orientation alignment of the fibers.

    2.3.3 Impact of a and L on Peclet number

    Parsheh’s study[11]shows that the rotating Peclet number ( )Pe is effective to assess the impact of the turbulence on the fiber anisotropic orientation. It can be defined as the ratio of the local velocity gradient to the rotational diffusion coefficient. The influence of the turbulence on the fiber orientation can be ignored when Pe> 10 and the fiber orientation tends to be isotropic when C<4 and Pe is very small. Here the Peclet number is used to explain effectively the effect of the wall shape on the flow mechanisms and the fiber orientation behavior.

    Fig.12 The variation of Pe along the cental streamline for different a and L=550mm

    Figure 12 shows that the Peclet number evolves along the center streamline for different inlet tangential angle a for a given channel length L =550mm . When a changes from a= -3 0oto a= 0o, the wall shape turns from a concave curve to a first-convexthen-concave curve, which directly varies the velocity gradient of the flow. In the first half of the contract field, the Peclet number is large when a = -3 0oand it is small when a=0owhile the situation is opposite in the second half field. Finally the Peclet numbers all remain at a low level to decrease the anisotropic orientation of the fibers significantly due to the fact thatb= 0oat the outlet except for the linear contraction. The fluctuations of the Peclet number curve just reflect the effect of the wall curvature on the velocity gradient in the flow field, so it can be inhibited to a certain extent by controlling the curvature of the wall geometry.

    The Peclet number varies similarly with the contraction ratio C, but much less than the linear wall in Fig.13, where the inlet tangential angle a retains zero degree and the channel length L takes different values from 350 mm to 550 mm. The evolution of the Peclet number with the contraction ratio reflects more intuitively that the velocity gradient in the flow has the trend of first increasing and then decreasing with C. As a result, the anisotropic orientation distribution of the fibers in all the parametric contractions is smaller than in the linear contraction.

    Fig.13 The variation of Pe with C along the cental streamline for different L and a=0o

    3. Conclusions

    The contraction ratio C is an important factor to determine the orientation distribution of suspending fibers, especially in case of a large speedup needed by the actual manufacturing. Though at present there are many methods to optimize the anisotropic orientation of the fibers, an effective method is still not easy to find to obtain a desired boundary shape based on the flow mechanism and the dynamics precision. Therefore a dynamic mesh model is established for the turbulent contraction using an adaptive grid algorithm and the quasi static hypothesis. Then the dynamics model of the turbulent flow with a variable boundary can be used by optimization objectives and the optimalparameters of the wall shape can be searched and obtained. Finally, the influence of the wall shape parameters including the inlet tangential angle, the channel length and the outlet tangential angle on the flow behavior and the fiber orientation distribution in the turbulence is analyzed and discussed. Several meaningful conclusions are drawn as follows:

    (1) The essence of an optimal wall shape is the optimization of the distribution of the local velocity and velocity gradient in the flow field, to affect the diffusion and orientation behavior of suspending fibers and then consequently to obtain the desired fiber orientation distribution. The inlet tangential angle a, the channel length L and the outlet tangential angle bare three parameters to control the contract shape of the wall. When b= 0oand the wall shape is very close to the linear contraction, a1111,maxachieves the minimum but a1111,out-1/p is not optimal. Whenb =0o, a = 0oand the channel length L is slightly short,a1111,out-1/p assumes the minimum value and is optimal. How to explain these results? The outlet is similar to a short straight channel at b= 0oand the tensile effect of the straight channel flow is much smaller than that of the contraction so that the fibers are no longer strongly under the stretching action of the flow field but take various orientations under the turbulent effect. The inlet tangential angle a and the channel length L can control the position of the maximum velocity gradient in the flow field, so they can adjust the value and the occurrence location of a1111,max.

    (2) The velocity gradient increases almost linearly with the contraction ratio and the tensile effect of the mean flow is much stronger than the random effect of the turbulence in the linear contract flow field, which make the fibers easy to align with the streamline consequently. In contrast, the variation of the velocity gradient with the boundary curvature is not monotonic in the flow field with the parametric surface. The non monotonic variation of the curvature would have the random effects of the turbulence enhanced to a certain extent and the fibers would have more chances of random orientations. Finally the orientation alignment of the fibers at the outlet can be improved in the contraction.

    [1] Batchelor G. K., Proudman I. The effect of rapid distortion of a fluid in turbulent motion [J]. Quarterly Journal of Mechanics and Applied Mathematics, 1954, 7(1): 83-103.

    [2] Ribner H., Tucker M. Spectrum of turbulence in a contracting stream [R]. Technical Report Archive and Image Library, 1953, 99-115.

    [3] Hussain A. K., Ramjee V. Effects of the axisymmetric contraction shape on incompressible turbulent flow [J]. Journal of Fluids Engineering, 1976, 98(1):58-68.

    [4] Lin J. Z., Zhang L. X. On the structural features of fiber suspensions in converging channel flow [J]. Journal of Zhejiang University-Science A, 2003, 4(4): 400-406.

    [5] Lin J., Shi X., You Z. Effects of the aspect ratio on the sedimentation of a fiber in Newtonian fluids [J]. Journal of Aerosol Science, 2003, 34(7): 909-921.

    [6] Lin J., Shi X., Yu Z. The motion of fibers in an evolving mixing layer [J]. International Journal of Multiphase Flow, 2003, 29(8): 1355-1372.

    [7] Lin J., Zhang W., Yu Z. Numerical research on the orientation distribution of fibers immersed in laminar and turbulent pipe flows [J]. Journal of Aerosol Science, 2004, 35(1): 63-82.

    [8] Lin J., Liang X., Zhang S. Fibre orientation distribution in turbulent fibre suspensions flowing through an axisymmetric contraction [J]. Canadian Journal of Chemical Engineering, 2011, 89(12): 1416-1425.

    [9] Olson J. A. The motion of fibre in turbulent flow, stochastic simulation of isotropic homogeneous turbulence [J]. International Journal of Multiphase Flow, 2001, 27(12): 2083-2103.

    [10] Olson J. A. Frigaard I., Chan C. et al. Modeling a turbulent fibre suspension flowing in a planar contraction: The onedimensional headbox [J]. International Journal of Multiphase Flow, 2004, 30(1): 51-66.

    [11] Parsheh M., Matthew L. B., Cyrusk A. On the orientation of stiff fibres suspended in turbulent flow in a planar contraction [J]. Journal of Fluid Mechanics, 2005, 545: 245-269.

    [12] Parsheh M., Brown M. L., Aidun C. K. Variation of fiber orientation in turbulent flow inside a planar contraction with different shapes [J]. International Journal of Multiphase Flow, 2006, 32(12): 1354-1369.

    [13] M?kinen R. A. E., H?m?l?inen J. Optimal control of a turbulent fibre suspension flowing in a planar contraction [J]. Communications in Numerical Methods in Engineering, 2010, 22: 567-575.

    [14] Nazemi A. R., Farahi M. H. Control the fibre orientation distribution at the outlet of contraction [J]. Acta Applicandae Mathematicae, 2009, 106(2): 279-292.

    [15] Farhadinia B. Planar contraction design with specified fibre orientation distribution at the outlet [J]. Structural and Multidisciplinary Optimization, 2012, 46(4): 503-511.

    [16] Ye S., Xiong J. Geometric optimization design system incorporating hybrid GRECO-WM scheme and genetic algorithm [J]. Chinese Journal of Aeronautics, 2009, 22(6): 599-606.

    [17] Zhang Z., Ding Y., LIU Y. Shape optimization design method for the conceptual design of reentry vehicles [J]. Acta Aeronautica Et Astronautica Sinica, 2011, 32(11): 1971-1979.

    [18] Yang W. A new method for fiber orientation distribution in a planar contracting turbulent flow [J]. Thermal Science, 2012, 16(5): 1367-1371.

    [19] Yang J. Sharp interface direct forcing immersed boundary methods: A summary of some algorithms and applications [J]. Journal of Hydrodynamics, 2016, 28(5): 713-730.

    [20] Feng Z. R., Wang L. X., Cui J. H. Research on the generation of two-dimensional fitted-body grid by adjacent boundary nodes sliding method [J]. Advanced Materials Research, 2013, 651: 835-840.

    [21] Xu W., Chen G., Deng Z. et al. The dynamical design of a steel frame based on the arbitrary Lagrange-Euler method [J]. Structures Under Shock and Impact XIII, 2014, 141: 51-57.

    10.1016/S1001-6058(16)60761-8

    February 27, 2016, Revised May 18, 2016)

    * Project supported by the National Natural Science Foundation of China (Grant No. 11302110), the Public Project of Science and Technology Department of Zhejiang Province (Grant No. 2015C31152), the Natural Science Foundation of Ningbo (Grant No. 2014A610086) and “Wang Weiming”Entrepreneurship Supporting Fund.

    Biography:Wei Yang (1975-), Female, Ph. D.,

    Associate Professor

    卡戴珊不雅视频在线播放| 美女大奶头视频| 午夜精品国产一区二区电影 | 午夜福利高清视频| 亚洲av二区三区四区| 中文字幕久久专区| 又黄又爽又刺激的免费视频.| 18+在线观看网站| 久久久久久久久久久免费av| 久久久久久久亚洲中文字幕| 亚洲综合色惰| 一级毛片电影观看| 国产免费视频播放在线视频 | 国产精品熟女久久久久浪| 三级毛片av免费| 久久精品国产亚洲网站| 国产v大片淫在线免费观看| 欧美xxxx性猛交bbbb| 乱系列少妇在线播放| 亚洲国产精品国产精品| 中国美白少妇内射xxxbb| 欧美日韩综合久久久久久| 免费观看av网站的网址| 国产午夜福利久久久久久| 日韩一本色道免费dvd| 国产有黄有色有爽视频| 欧美不卡视频在线免费观看| 特级一级黄色大片| 免费看光身美女| 国内精品一区二区在线观看| 国产黄片美女视频| 国产一区二区三区综合在线观看 | 国产 一区精品| 日韩不卡一区二区三区视频在线| 少妇裸体淫交视频免费看高清| 亚洲最大成人手机在线| 免费人成在线观看视频色| 免费观看av网站的网址| 丰满人妻一区二区三区视频av| 欧美成人a在线观看| 免费看日本二区| 中文资源天堂在线| 中文欧美无线码| 97在线视频观看| 国产片特级美女逼逼视频| 丝袜喷水一区| .国产精品久久| eeuss影院久久| 国产精品女同一区二区软件| 国产亚洲一区二区精品| 久久亚洲国产成人精品v| 中文字幕av成人在线电影| 久久久国产一区二区| 18禁裸乳无遮挡免费网站照片| 丝瓜视频免费看黄片| 不卡视频在线观看欧美| 在线观看一区二区三区| 嫩草影院入口| 十八禁国产超污无遮挡网站| 一级毛片黄色毛片免费观看视频| 日韩欧美三级三区| 亚洲欧美中文字幕日韩二区| 亚洲av中文字字幕乱码综合| 国产女主播在线喷水免费视频网站 | 一二三四中文在线观看免费高清| 丝瓜视频免费看黄片| 99久久九九国产精品国产免费| 亚洲精品成人久久久久久| 国产综合精华液| 三级国产精品片| 国产成人午夜福利电影在线观看| 91aial.com中文字幕在线观看| a级毛片免费高清观看在线播放| 51国产日韩欧美| 久久精品熟女亚洲av麻豆精品 | 亚洲国产精品成人综合色| 干丝袜人妻中文字幕| 91狼人影院| 国产精品福利在线免费观看| 91久久精品国产一区二区成人| 国产精品伦人一区二区| av福利片在线观看| 777米奇影视久久| 在现免费观看毛片| a级一级毛片免费在线观看| 精品人妻偷拍中文字幕| 搡老妇女老女人老熟妇| 亚洲图色成人| 亚洲av日韩在线播放| 国产男人的电影天堂91| 亚洲精品日本国产第一区| 日韩视频在线欧美| 久久热精品热| 亚洲av电影在线观看一区二区三区 | 免费在线观看成人毛片| 国产探花极品一区二区| av免费观看日本| 丝袜美腿在线中文| 国产亚洲av片在线观看秒播厂 | 搡老妇女老女人老熟妇| 免费av毛片视频| 99九九线精品视频在线观看视频| 女的被弄到高潮叫床怎么办| 99久久九九国产精品国产免费| 国内揄拍国产精品人妻在线| 最后的刺客免费高清国语| av黄色大香蕉| 美女脱内裤让男人舔精品视频| 菩萨蛮人人尽说江南好唐韦庄| 亚洲婷婷狠狠爱综合网| 欧美97在线视频| 午夜福利高清视频| 欧美日韩视频高清一区二区三区二| 啦啦啦啦在线视频资源| 不卡视频在线观看欧美| 爱豆传媒免费全集在线观看| 国产麻豆成人av免费视频| 日韩在线高清观看一区二区三区| 亚洲精品,欧美精品| 一级二级三级毛片免费看| 国产亚洲一区二区精品| 欧美区成人在线视频| 亚洲av免费在线观看| 在现免费观看毛片| 91狼人影院| 3wmmmm亚洲av在线观看| 久久国内精品自在自线图片| 人人妻人人看人人澡| a级毛色黄片| 成年女人在线观看亚洲视频 | 亚洲内射少妇av| 亚洲精品第二区| 亚洲精华国产精华液的使用体验| 国产激情偷乱视频一区二区| 国产亚洲5aaaaa淫片| 久久综合国产亚洲精品| 九九爱精品视频在线观看| .国产精品久久| 免费不卡的大黄色大毛片视频在线观看 | 国内精品宾馆在线| kizo精华| 草草在线视频免费看| 亚洲国产最新在线播放| 久久久欧美国产精品| 日韩欧美精品v在线| 欧美zozozo另类| 国产乱人偷精品视频| 免费观看av网站的网址| 免费看美女性在线毛片视频| 国产大屁股一区二区在线视频| 天堂中文最新版在线下载 | 国产精品女同一区二区软件| videos熟女内射| 国产精品一区二区性色av| 精品一区在线观看国产| 久久精品国产亚洲av天美| 波野结衣二区三区在线| 男女边吃奶边做爰视频| 久久精品久久久久久久性| www.av在线官网国产| 成人国产麻豆网| 男插女下体视频免费在线播放| 免费少妇av软件| 亚洲av中文av极速乱| 免费观看性生交大片5| 精品亚洲乱码少妇综合久久| 天天躁夜夜躁狠狠久久av| 噜噜噜噜噜久久久久久91| 韩国高清视频一区二区三区| 国产高清不卡午夜福利| 亚洲自偷自拍三级| 免费看av在线观看网站| 肉色欧美久久久久久久蜜桃 | 中文乱码字字幕精品一区二区三区 | 日韩制服骚丝袜av| 最近2019中文字幕mv第一页| 99久久九九国产精品国产免费| 97超视频在线观看视频| 99热这里只有是精品在线观看| www.av在线官网国产| 日韩欧美三级三区| 男人狂女人下面高潮的视频| 国产成人精品久久久久久| 久久久精品94久久精品| 亚洲精品456在线播放app| 大陆偷拍与自拍| 亚洲熟妇中文字幕五十中出| 午夜福利成人在线免费观看| 成人毛片a级毛片在线播放| 美女xxoo啪啪120秒动态图| 久久久久久久国产电影| 久久久久久久久久人人人人人人| 国产亚洲最大av| 欧美 日韩 精品 国产| .国产精品久久| 有码 亚洲区| 特大巨黑吊av在线直播| 国产爱豆传媒在线观看| 一夜夜www| 国产高清三级在线| 看非洲黑人一级黄片| 久久午夜福利片| 80岁老熟妇乱子伦牲交| 亚洲国产精品专区欧美| 三级国产精品片| av天堂中文字幕网| 色综合站精品国产| 中文在线观看免费www的网站| 22中文网久久字幕| 性插视频无遮挡在线免费观看| 又黄又爽又刺激的免费视频.| 亚洲av二区三区四区| 亚洲av日韩在线播放| 欧美97在线视频| 99久久九九国产精品国产免费| av播播在线观看一区| 97在线视频观看| 亚洲精品日本国产第一区| 伊人久久国产一区二区| 亚洲av中文av极速乱| 亚洲aⅴ乱码一区二区在线播放| 国产色爽女视频免费观看| 国产精品国产三级专区第一集| av卡一久久| 国产免费又黄又爽又色| av国产免费在线观看| 春色校园在线视频观看| 亚洲最大成人中文| 亚洲一级一片aⅴ在线观看| 久久99热6这里只有精品| h日本视频在线播放| 国产精品三级大全| 天堂影院成人在线观看| 自拍偷自拍亚洲精品老妇| 一区二区三区免费毛片| 熟妇人妻久久中文字幕3abv| 床上黄色一级片| 18禁动态无遮挡网站| 亚洲电影在线观看av| 欧美不卡视频在线免费观看| 爱豆传媒免费全集在线观看| 99久国产av精品| 国内精品美女久久久久久| 超碰av人人做人人爽久久| 又爽又黄a免费视频| 亚洲精品久久久久久婷婷小说| 97热精品久久久久久| a级一级毛片免费在线观看| 久久这里只有精品中国| 成年免费大片在线观看| 91午夜精品亚洲一区二区三区| 国产成人a∨麻豆精品| 亚洲人成网站在线观看播放| 久久精品久久久久久噜噜老黄| 成人av在线播放网站| 波野结衣二区三区在线| 91在线精品国自产拍蜜月| 精品一区二区三卡| 亚洲欧美成人综合另类久久久| 少妇熟女欧美另类| 日韩欧美 国产精品| 亚洲婷婷狠狠爱综合网| 少妇熟女aⅴ在线视频| 一级片'在线观看视频| 狠狠精品人妻久久久久久综合| 日韩欧美 国产精品| 性插视频无遮挡在线免费观看| 18+在线观看网站| 色吧在线观看| 蜜桃亚洲精品一区二区三区| 日韩欧美 国产精品| 精品久久久久久久末码| 国产午夜精品论理片| 中文精品一卡2卡3卡4更新| 如何舔出高潮| 精品人妻偷拍中文字幕| 五月天丁香电影| 网址你懂的国产日韩在线| 淫秽高清视频在线观看| 中文字幕久久专区| 亚洲一级一片aⅴ在线观看| 久热久热在线精品观看| 亚洲激情五月婷婷啪啪| 一本久久精品| 亚洲色图av天堂| 综合色av麻豆| 亚洲人成网站高清观看| 久久久久久久久大av| 亚洲精品中文字幕在线视频 | 丰满人妻一区二区三区视频av| 久热久热在线精品观看| av播播在线观看一区| 在线 av 中文字幕| 亚洲国产高清在线一区二区三| 久久99热这里只频精品6学生| 一级毛片aaaaaa免费看小| 精品久久久久久电影网| 99久国产av精品| 亚洲激情五月婷婷啪啪| 乱系列少妇在线播放| 亚洲电影在线观看av| 两个人视频免费观看高清| 中文字幕av在线有码专区| 日韩大片免费观看网站| 尤物成人国产欧美一区二区三区| 亚洲欧美成人综合另类久久久| 色视频www国产| 1000部很黄的大片| 日本欧美国产在线视频| 国产有黄有色有爽视频| 大片免费播放器 马上看| 国产成人a区在线观看| 日韩av在线免费看完整版不卡| 免费观看性生交大片5| 免费黄频网站在线观看国产| 国产女主播在线喷水免费视频网站 | 超碰97精品在线观看| 我的老师免费观看完整版| 丝袜美腿在线中文| 久久久久国产网址| 日本爱情动作片www.在线观看| 久久久欧美国产精品| 中文精品一卡2卡3卡4更新| 亚洲国产精品国产精品| 日韩制服骚丝袜av| 亚洲图色成人| 草草在线视频免费看| 蜜桃亚洲精品一区二区三区| 国产精品无大码| 国产中年淑女户外野战色| 国产午夜福利久久久久久| 一区二区三区四区激情视频| 十八禁国产超污无遮挡网站| 天天躁夜夜躁狠狠久久av| 一夜夜www| 内地一区二区视频在线| 日韩在线高清观看一区二区三区| 看非洲黑人一级黄片| 精品一区二区三区视频在线| 纵有疾风起免费观看全集完整版 | 黄色配什么色好看| 观看美女的网站| 免费av观看视频| 久久精品综合一区二区三区| 日韩一区二区视频免费看| 亚洲自偷自拍三级| 免费观看精品视频网站| 一边亲一边摸免费视频| 成年女人在线观看亚洲视频 | 免费观看的影片在线观看| av天堂中文字幕网| 亚洲国产最新在线播放| 九九在线视频观看精品| videossex国产| 国产精品.久久久| 国产高清有码在线观看视频| 九草在线视频观看| 免费观看无遮挡的男女| 亚洲av二区三区四区| 国精品久久久久久国模美| 亚洲,欧美,日韩| 成人漫画全彩无遮挡| 免费av观看视频| 成人特级av手机在线观看| 欧美潮喷喷水| 一级a做视频免费观看| 成人一区二区视频在线观看| 麻豆国产97在线/欧美| 亚洲av成人精品一区久久| 免费电影在线观看免费观看| 国产一区有黄有色的免费视频 | 中文字幕制服av| 日韩中字成人| av线在线观看网站| av免费在线看不卡| av在线老鸭窝| 亚洲18禁久久av| 国产 一区 欧美 日韩| 美女黄网站色视频| 亚洲在线自拍视频| 一边亲一边摸免费视频| 成人午夜精彩视频在线观看| 亚洲最大成人av| 亚洲精品一二三| 久久久精品94久久精品| 精品熟女少妇av免费看| 午夜精品在线福利| 亚洲第一区二区三区不卡| 亚洲欧美清纯卡通| av在线蜜桃| 毛片女人毛片| 国内精品美女久久久久久| 一本一本综合久久| 国产免费一级a男人的天堂| 禁无遮挡网站| 欧美xxⅹ黑人| 亚洲欧美中文字幕日韩二区| 神马国产精品三级电影在线观看| 国产精品一二三区在线看| 国产极品天堂在线| 亚洲精品456在线播放app| 91av网一区二区| 午夜精品国产一区二区电影 | 丰满人妻一区二区三区视频av| 国产综合懂色| 高清av免费在线| 可以在线观看毛片的网站| 麻豆国产97在线/欧美| a级一级毛片免费在线观看| 久久久精品94久久精品| 伦精品一区二区三区| 一夜夜www| 搡老妇女老女人老熟妇| 亚洲av二区三区四区| 国模一区二区三区四区视频| 你懂的网址亚洲精品在线观看| 国产亚洲av嫩草精品影院| 欧美三级亚洲精品| 国产淫片久久久久久久久| ponron亚洲| 国产成人aa在线观看| 免费高清在线观看视频在线观看| 国产单亲对白刺激| 深夜a级毛片| 久久久久久久大尺度免费视频| 日韩中字成人| 天堂影院成人在线观看| 亚洲内射少妇av| 亚洲人成网站在线观看播放| 嫩草影院精品99| 国产伦理片在线播放av一区| 女人十人毛片免费观看3o分钟| 国产成人精品福利久久| 丰满少妇做爰视频| 人妻夜夜爽99麻豆av| 午夜免费激情av| 亚洲电影在线观看av| 中文乱码字字幕精品一区二区三区 | 日韩av在线大香蕉| 人体艺术视频欧美日本| 国产精品一区www在线观看| 日韩av不卡免费在线播放| 99热这里只有是精品在线观看| 久久国产乱子免费精品| 直男gayav资源| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 日韩视频在线欧美| or卡值多少钱| 欧美不卡视频在线免费观看| 婷婷色综合www| 身体一侧抽搐| 免费在线观看成人毛片| 看非洲黑人一级黄片| 日韩欧美国产在线观看| 有码 亚洲区| 久久亚洲国产成人精品v| 国产黄片视频在线免费观看| 久久99精品国语久久久| 能在线免费观看的黄片| 国产亚洲5aaaaa淫片| 欧美成人精品欧美一级黄| 亚洲精品自拍成人| 精品一区二区三区视频在线| 欧美日韩综合久久久久久| 一二三四中文在线观看免费高清| 日韩一区二区视频免费看| 我的女老师完整版在线观看| 中文字幕制服av| 国产一区二区三区av在线| 麻豆精品久久久久久蜜桃| 精品酒店卫生间| 超碰av人人做人人爽久久| 国产成人午夜福利电影在线观看| 国产黄a三级三级三级人| 丝瓜视频免费看黄片| 看非洲黑人一级黄片| 一级黄片播放器| 1000部很黄的大片| 亚洲av一区综合| 国产精品不卡视频一区二区| 国国产精品蜜臀av免费| 大又大粗又爽又黄少妇毛片口| 菩萨蛮人人尽说江南好唐韦庄| 99热这里只有精品一区| 看非洲黑人一级黄片| 一级毛片久久久久久久久女| 一级毛片aaaaaa免费看小| 久久这里只有精品中国| 欧美极品一区二区三区四区| 日韩中字成人| 精品久久久久久成人av| 26uuu在线亚洲综合色| 国产v大片淫在线免费观看| 免费高清在线观看视频在线观看| 国产美女午夜福利| 三级国产精品欧美在线观看| 午夜福利视频1000在线观看| 三级经典国产精品| 69av精品久久久久久| 国产成人freesex在线| 日本免费a在线| 国产激情偷乱视频一区二区| 少妇的逼好多水| 老师上课跳d突然被开到最大视频| 日韩伦理黄色片| 午夜精品国产一区二区电影 | 亚洲无线观看免费| 麻豆国产97在线/欧美| 老女人水多毛片| 男女下面进入的视频免费午夜| 女的被弄到高潮叫床怎么办| 深夜a级毛片| 国产亚洲一区二区精品| 久久久午夜欧美精品| 一本久久精品| 18+在线观看网站| 嫩草影院精品99| kizo精华| 嫩草影院精品99| 午夜免费观看性视频| 99久久精品热视频| 免费看av在线观看网站| 一个人看的www免费观看视频| 五月玫瑰六月丁香| 成人亚洲精品av一区二区| 精品人妻视频免费看| 精品欧美国产一区二区三| 99久久九九国产精品国产免费| 成人午夜高清在线视频| 亚洲,欧美,日韩| 国产黄a三级三级三级人| 最近最新中文字幕大全电影3| 精品一区二区免费观看| 国产久久久一区二区三区| 深夜a级毛片| 欧美丝袜亚洲另类| 久久久久久久久中文| 亚洲精品成人久久久久久| 久久精品久久久久久久性| 国产精品无大码| 国产乱来视频区| 亚洲色图av天堂| 天天躁日日操中文字幕| 色网站视频免费| a级毛片免费高清观看在线播放| 亚洲三级黄色毛片| 欧美精品国产亚洲| 精品久久久久久久久久久久久| 国产 亚洲一区二区三区 | 91在线精品国自产拍蜜月| 国产成人freesex在线| 女人十人毛片免费观看3o分钟| 国产白丝娇喘喷水9色精品| 嫩草影院入口| 亚洲成色77777| 啦啦啦中文免费视频观看日本| 99久久九九国产精品国产免费| 日韩av在线大香蕉| 国产成人精品福利久久| 亚洲精品影视一区二区三区av| 免费大片黄手机在线观看| av在线天堂中文字幕| 国产激情偷乱视频一区二区| 国产成人精品婷婷| 成年av动漫网址| 26uuu在线亚洲综合色| 亚洲欧美清纯卡通| 人妻夜夜爽99麻豆av| 亚洲成人久久爱视频| 国产爱豆传媒在线观看| 人妻系列 视频| 亚洲成人中文字幕在线播放| 亚洲不卡免费看| 亚州av有码| 久久久精品免费免费高清| 欧美性猛交╳xxx乱大交人| 日韩,欧美,国产一区二区三区| 久久韩国三级中文字幕| 欧美日韩综合久久久久久| 国产亚洲精品久久久com| 纵有疾风起免费观看全集完整版 | 国产欧美日韩精品一区二区| 91久久精品国产一区二区三区| 欧美人与善性xxx| 久久韩国三级中文字幕| 亚洲国产色片| 国产高潮美女av| 久久人人爽人人片av| 亚洲美女视频黄频| 国产精品不卡视频一区二区| 大香蕉97超碰在线| 欧美bdsm另类| 一级毛片aaaaaa免费看小| 嫩草影院入口| 国产精品福利在线免费观看| 国产午夜精品一二区理论片| 美女被艹到高潮喷水动态| 九草在线视频观看| 一边亲一边摸免费视频| 亚洲最大成人中文| 一级毛片电影观看| 国产伦理片在线播放av一区| 日韩av在线大香蕉| 人妻夜夜爽99麻豆av| 欧美极品一区二区三区四区| 国产成年人精品一区二区| 国产伦一二天堂av在线观看| 国产在视频线精品| 亚洲熟女精品中文字幕| 欧美一级a爱片免费观看看| 18禁在线播放成人免费| 国产熟女欧美一区二区| 久久99热6这里只有精品| 久久久成人免费电影| 国产又色又爽无遮挡免|