#### Table of Content

18 May 2019, Volume 51 Issue 3
Research Review
A REVIEW ON THE DEFORMATION AND FAILURE OF SEALS IN OILFIELD1)
Zhengjin Wang, Tiejun Wang
2019, 51(3):  635-655.  DOI: 10.6052/0459-1879-19-061
Asbtract ( 1063 )   HTML ( 52)   PDF (45847KB) ( 558 )
Figures and Tables | References | Related Articles | Metrics

Oil and gas will still be the major energy resource of human in a long period. Sealing is one of the most important operations in oil and gas exploitation, which plays an essential role in horizontal drilling and hydraulic fracturing, booming the exploitation of unconventional oil and gas like shell gas. The quality of sealing is critical to the performance and safety of wells. The deformation and failure of seals may lead to severe environmental pollution, losses of life and property. Therefore, it is important to study the mechanical behavior of seals in oilfields systematically and deeply for safety assessment as well as better sealing of wells. Cement and swellable elastomers are the most widely used sealing materials in oilfield. Well cementing has been developed for more than 100 years. While swellable elastomeric packer is relatively new. In this review, we present the recent advances in the studies of the deformation and failure of cement sheath and swellable packers. The first part focuses on the cement sheath, including stress evolution during the solidification process and production stage, failure modes, and failure criteria. In the second part, emphasis is placed on the work principle, basic theory, test method, failure modes and instability of swellable packers.

RESEARCH PROGRESS IN AUXETIC MATERIALS AND STRUCTURES1)
Xin Ren, Xiangyu Zhang, Yimin Xie
2019, 51(3):  656-689.  DOI: 10.6052/0459-1879-18-381
Asbtract ( 2113 )   HTML ( 101)   PDF (1758KB) ( 1525 )
Figures and Tables | References | Related Articles | Metrics

Auxetic materials and structures have special mechanical properties. In contrast to traditional materials, auxetic materials contract (expand) under uniaxial compression (tension). Auxetic materials have many desirable properties and are superior to traditional materials in terms of shear capacity, fracture resistance, energy absorption capacity and indentation resistance. Because of these excellent properties, auxetic materials are very promising in the fields of medical equipment, intelligent filters, intelligent sensors, protective equipment, aviation, navigation and national defense engineering, but the application and popularization of auxetic materials still face some challenges. In this paper, the domestic and foreign research results of auxetic materials are extensively reviewed and the latest progress of auxetic materials is introduced. Auxetic materials are roughly classified into the following four categories: natural auxetic materials, cellular auxetic materials, metallic auxetic materials, multi-material auxetics and auxetic composites. Internal structure, auxetic deformation mechanism and mechanical properties of each kind are discussed. Some excellent mechanical properties of auxetic materials are summarized and new inventions and applications of auxetic materials in various industries are introduced as well. In view of the fact that the actual application of auxetic materials are much less than the research theories and experimental results, this paper analyzes some shortcomings and the current challenges of auxetic materials, including high manufacturing cost, high porosity, insufficient bearing capacity and only suitable for small strain conditions. To promote the application of auxetics, the metallic auxetic materials are introduced in detail. The steps and methods for designing and fabricating 3D metallic auxetic materials are proposed as well.

Theme Articles on
Preface
Zhang Xiong
2019, 51(3):  688-689.
Asbtract ( 174 )   PDF (63KB) ( 99 )
Related Articles | Metrics
A GRADIENT SMOOTHING GALERKIN MESHFREE METHOD FOR THIN PLATE ANALYSIS WITH LINEAR BASIS FUNCTION1)
Like Deng, Dongdong Wang, Jiarui Wang, Junchao Wu
2019, 51(3):  690-702.  DOI: 10.6052/0459-1879-19-004
Asbtract ( 577 )   HTML ( 14)   PDF (12603KB) ( 207 )
Figures and Tables | References | Related Articles | Metrics

The fourth order governing equation of thin plate necessitates the employment of C$^{1}$ continuous shape functions with a minimum degree of two in a Galerkin formulation. Thus at least a quadratic basis function should be utilized in meshfree approximation to enable the Galerkin meshfree thin plate analysis. However, due to the rational nature of reproducing kernel meshfree shape functions, the computation of the second order derivatives of meshfree shape functions is quite complex and costly, which also requires expensive high order Gauss quadrature rules to properly integrate the stiffness matrix. In this work, a gradient smoothing Galerkin meshfree method with particular reference to the linear basis function is proposed for thin plate analysis. The foundation of the present development is the construction of smoothed meshfree gradients with linear basis function, where the second order smoothed gradients are expressed as combinations of standard first order gradients and the computational burden is remarkably reduced. Furthermore, it is shown that the smoothed meshfree gradients with linear basis function satisfy both the linear and quadratic gradient consistency conditions and consequently they are adequate for thin plate analysis in the context of Galerkin formulation. An interpolation error study is given as well to validate the higher order consistency conditions and applicability of smoothed meshfree gradients for Galerkin analysis of thin plates. It turns out that efficient lower order Gauss integration rules now work well for the proposed method. Numerical results demonstrate that compared with the conventional Galerkin meshfree method with quadratic basis function, the proposed gradient smoothing Galerkin meshfree method with linear basis function yields similar convergence rates, but with better accuracy and less integration points for stiffness computation.

FREE ELEMENT METHOD AND ITS APPLICATION IN STRUCTURAL ANALYSIS1)
Xiaowei Gao, Bingbing Xu, Jun Lü, Haifeng Peng
2019, 51(3):  703-713.  DOI: 10.6052/0459-1879-19-011
Asbtract ( 359 )   HTML ( 10)   PDF (9187KB) ( 218 )
Figures and Tables | References | Related Articles | Metrics

By absorbing advantages of the finite element and meshless methods, a new numerical method, free element method, is proposed in the paper. In the discretization, the isoparametric elements as used in FEM are employed to represent the geometry and interpolate physical variables; and in the algorithm, the point collocation technique using elements is employed to generate the system of equations point by point. The main feature of the method is that only one independent element formed by freely selecting surrounding points is required for each collocation point, without need to consider the connective relationship between adjacent elements and the continuity of physical variables and their spatial derivatives at interfaces of the connected elements. Two types of free element methods, the strong-form method and weak form method, will be described in the paper. The former directly generates the system of equations from the governing equations and the Neumann boundary conditions, while the latter establishes the weak-form integral expression of the governing equations by the weighted residual technique over the free element first and then generates the system of equations through an integration process similar to that employed in the standard FEM. The method proposed in the paper is an element collocation method. To achieve highly accurate spatial derivatives for internal collocation points of the computational domain, isoparametric elements with at least one internal node are required. For this purpose, apart from the arbitrary order quadrilateral Lagrange elements, a new seven-node triangle element is constructed in the paper, which can be used to model problems with complex geometries.

NUMERICAL SIMULATION OF LIQUID SLOSHING IN LNG TANK USING GPU-ACCELERATED MPS METHOD1)
Xiang Chen, Decheng Wan
2019, 51(3):  714-729.  DOI: 10.6052/0459-1879-18-410
Asbtract ( 397 )   HTML ( 14)   PDF (35394KB) ( 157 )
Figures and Tables | References | Related Articles | Metrics

Liquid sloshing is a common phenomenon induced in partially filled tanks under external excitations, which may destroy the tank structure and vessel stability. Moving particle semi-implicit (MPS) method is a typical meshfree method which can effectively simulate violent liquid sloshing problem. However, the low computational efficiency of MPS is the bottleneck of its application in large-scale three-dimensional problems. In the past years, GPU parallel acceleration technique has been widely used in the field of scientific computing. In this work, GPU parallel acceleration technique is introduced into MPS method and an in-house solver MPSGPU-SJTU is developed by using CUDA language. Then this solver is used to simulate 3-D liquid sloshing in liquefied natural gas (LNG) tank. The convergent validation of particle spacing is carried out to verify the accuracy of present solver. The maximum particle number of simulation model is over two million particles. MPSGPU-SJTU solver can accurately predict the impact pressures by comparing with other results. In addition, the violent flow phenomena such as large deformation and nonlinear fragmentation of free surface can be observed in these simulations. The comparison of computation time between GPU and CPU solvers demonstrates GPU parallel acceleration technique can significantly reduce the computation time and improve the computational efficiency of MPS. The phenomena of liquid sloshing in LNG tank and rectangular tank are compared. The results show that LNG tank can reduce the sloshing amplitude and impact pressure in high filling level. However, the sloshing is more violent and the free surface presents three-dimensional feature in LNG tank with middle and low filling level. Finally, the investigation of the effect of different fluids such as water and LNG on sloshing phenomena is also conducted in this paper. It shows that the flow fields of both liquids are almost similar and the impact pressure is proportional to the liquid density.

AN IMPROVED SPH ALGORITHM FOR LARGE DENSITY RATIOS MULTIPHASE FLOWS BASED ON RIEMANN SOLUTION1)
Qiuzu Yang, Fei Xu, Lu Wang, Yang Yang
2019, 51(3):  730-742.  DOI: 10.6052/0459-1879-18-451
Asbtract ( 401 )   HTML ( 13)   PDF (34100KB) ( 179 )
Figures and Tables | References | Related Articles | Metrics

The discontinuities of physical fields (such as density, viscosity and so on) exist in the different phase interface of multiphase flow problems, and the numerical simulation method using the traditional SPH model is liable to cause spurious~oscillations in the pressure and velocity field at the interface, which is a big problem in the application of multiphase flows. An improved SPH model based on Riemann solution on dealing with the abrupt physical quantities of multiphase flows with large density radios is presented. Using the advantage of Riemann solution in dealing with the contact discontinuity problems, we introduce it into the SPH multiphase flow model. For the sake of accurately calculating the physical viscosity of multiphase fluid and decreasing the Riemann dissipation, the SPH momentum equation of the Riemann solution form is improved. In the new model, we combine the Adami fixed particle wall-boundary with the one-sided Riemann problem to impose solid boundary of the SPH multiphase flow, and consider the influence of the surface tension on the small-scale interface. The new model without adding any artificial viscosity, artificial dissipation and non-physical treatment technology can simulate the real physical viscousity and the physical evolution process of multiphase flow problems. In order to verify the ability of the improved model in dealing with the multiphase flow problems with the discontinuous interface and the convergence of the model particle spacing, firstly the squared droplet oscillating problem is simulated under different discrete particle spacing. Afterwards, the multiphase flow problems of the Rayleigh-Taylor instability, the single bubble buoyancy and the double bubble buoyancy are simulated. The interface is clearly capture and the results are good in agreement with the literature, which proved that the improved multiphase flow SPH model can stably and effectively deal with the problems of the multiphase flow large with density ratio and large viscosity ratio.

FINTE SUBDOMAIN RADIAL BASIS COLLOCATION METHOD FOR THE LARGE DEFORMATION ANALYSIS1)
Lihua Wang, Yiming Li, Fuyun Chu
2019, 51(3):  743-753.  DOI: 10.6052/0459-1879-19-005
Asbtract ( 385 )   HTML ( 13)   PDF (5440KB) ( 200 )
Figures and Tables | References | Related Articles | Metrics

The meshfree methods can avoid grid distortion problems because it does not need to be meshed, which make them widely used in large deformations and other complicated problems. Radial basis collocation method (RBCM) is a typical strong form meshfree method. This method has the advantages of no need for any mesh, simple solution process, high precision, good convergence and easy expansion to high-dimensional problems. Since the global shape function is used, this method has the disadvantages of low precision and poor representation to local characteristics when solving high gradient problems. To resolve this issue, this paper introduces finite subdomain radial basis collocation method to solve the large deformation problem with local high gradients. Based on the total Lagrangian formulation, the Newton iteration method is used to establish the incremental solution scheme of the FSRBCM in large deformation analysis. This method partitions the solution domain into several subdomains according to its geometric characteristics, then constructs radial basis function interpolation in the subdomains, and imposes all the interface continuous conditions on the interfaces, which results in a block sparse matrix for the numerical solution. The proposed method has super convergence and transforms the original full matrix into a sparse matrix, which reduces the storage space and improves the computational efficiency. Compared to the traditional RBCM and finite element method (FEM), this method can better reflect local characteristics and solve high gradient problems. Numerical simulations show that the method can effectively solve the large deformation problems with local high gradients.

Fluid Mechanics
LARGE-EDDY SIMULATION OF FLOWS WITH COMPLEX GEOMETRIES BY USING THE SLIP-WALL MODEL1)
Beiji Shi, Guowei He, Shizhao Wang
2019, 51(3):  754-766.  DOI: 10.6052/0459-1879-19-033
Asbtract ( 491 )   HTML ( 21)   PDF (10807KB) ( 249 )
Figures and Tables | References | Related Articles | Metrics

A slip-wall model is combined with the immersed boundary method for large-eddy simulation of flows with complex geometries. Firstly, large eddy simulation of flows over periodic hills are conducted to evaluate the effects of the tangential pressure gradient in wall model. Both the equilibrium stress balance model which based on the assumption of an equilibrium boundary-layer and the non-equilibrium wall model, in which the pressure gradient blend into the simplified thin boundary-layer equations and the RANS-like eddy viscosity in both the procedure of computing the wall-shear stress and reconstructing the wall-slip velocity, are utilized for comparison. The numerical results show that the pressure coefficient is not sensitive to the types of wall model which we considered, especially the strong pressure gradient in front of the hill crest is well catched by both models. However, when taking into account the tangential pressure gradient, the non-equilibrium wall model is superior to the equilibrium one for its ability to improve the prediction of the wall-shear stress and flow separation. When the equilibrium stress balance model is used, the wall-shear stress is heavily under-predicted and remarkable discrepancies of the mean velocity profiles can also be seen in the recirculation region. By comparison, the correction of the non-equilibrium wall model is proportional to both the tangential pressure gradient and the normal distance away from the wall, thus the hydrodynamic coefficients and the mean flow statistics are all in good agreement with the references even on very coarse grids. Secondly, large-eddy simulation of flow around an axisymmetric body is conducted to assess the applicability of current method when applied to high Reynolds number wall-bounded turbulent flows. The flow structures and the hydrodynamic characteristics are well predicted by the non-equilibrium wall model. This work confirms that the immersed boundary method in combination with the non-equilibrium slip-wall model is a possible and promising way to deal with turbulent flows which have complex geometries.

THE EFFECT OF PERIODIC PERTURBATION ON MULTI SCALES IN A TURBULENT BOUNDARY LAYER FLOW UNDER DRAG REDUCTION1)
Shuaijie Wang, Xiaotong Cui, Jianxia Bai, Zhanqi Tang, Nan Jiang
2019, 51(3):  767-774.  DOI: 10.6052/0459-1879-18-409
Asbtract ( 459 )   HTML ( 9)   PDF (6838KB) ( 106 )
Figures and Tables | References | Related Articles | Metrics

This study reports the drag reduction of a turbulent boundary layer flow under the active control by a pair of wall-mounted piezoelectric oscillators with different working frequency and amplitude. Under the maximum local skin drag reduction case of 80V and 160 Hz, a miniature boundary layer hot-wire probe of single sensor was used to measure velocity signals at the streamwise location 2 mm downstream of the oscillators. The time series of streamwise velocity signals at different wall normal positions were obtained. By comparing the turbulent fluctuations at different scales in the perturbed cases and canonical turbulent boundary layer flow, it was found that the perturbation produced by piezoelectric oscillators have an effect on the near-wall region, which attenuates the large-scale intensity and enhances the small-scale turbulence. Meanwhile, the outer layer is insensitive to the oscillator perturbation. Furthermore, conditional average of small-scale fluctuations shows that the impact of piezoelectric oscillators on small-scale fluctuations is not uniform in the condition of the large scales. The small-scale amplitude is increased more obviously as the large-scale fluctuations are positive than that the large-scale fluctuations are negative. It indicates that the current periodic perturbation dominantly modifies the high skin friction events by breaking the large-scale high-speed sweeping fluids into small-scale structures, which results in the local skin drag reduction.

AN INTERFACE TRACKING METHOD OF COUPLED YOUNGS-VOF AND LEVEL SET BASED ON GEOMETRIC RECONSTRUCTION1)
Manman Zhang, Jiao Sun, Wenyi Chen
2019, 51(3):  775-786.  DOI: 10.6052/0459-1879-18-439
Asbtract ( 568 )   HTML ( 5)   PDF (18627KB) ( 113 )
Figures and Tables | References | Related Articles | Metrics

Abstract The Lagrangian method and the Euler-Lagrangian method in the interface tracking method have low computational efficiency. They are not suitable for large deformation, and can not be applied to three-dimensional numerical calculation models. An Euler motion interface tracking method is applied to calculate the migration characteristics of gas-liquid two-phase interface with high efficiency and clear interface and suitable for three-dimensional models. This method uses the "米" shaped adjacent grid Youngs method for interface reconstruction. The Youngs-VOF and level set are coupled by geometric methods, so it can improve the accuracy of the interface, and overcome the defects of the VOF method and level set method, and avoid to solve the level set convection equation and the distance function equation with its own stability of the high-order derivation. The "米" shaped adjacent grid Youngs method is used to avoid the situation that the captured interface was blurred due to numerical dissipation and numerical dispersion, as well as nonlinear effects. The Youngs-VOF coupled level set method not only ensures the stability of the computing interface, but improves the computational efficiency compared with the Lagrangian method. The Youngs-VOF coupled level set method and VOF method is designed to simulate the rising process of a single bubble in water. After compared, the Youngs-VOF coupled level set method is more efficient and sharper than the VOF method. By using the Youngs-VOF coupled level set method and the VOF method to numerically calculate the circular motion interface model in the classical shear flow field, it is verified that the Youngs-VOF coupled level set method can better calculate the interface curvature compared with the VOF method. Through the numerical calculation of the break dam-free surface flow process and compared with the experiment, the stability of the Youngs-VOF coupled level set method are verified. And the Youngs-VOF coupled level set method can be applied to the three-dimensional numerical models.

STUDY OF FLOW-INDUCED MOTION CHARACTERISTICS OF THREE TANDEM CIRCULAR CYLINDERS IN PLANAR SHEAR FLOW1)
Jiahuang Tu, Xiaoling Tan, Xuhui Deng, Xiaogang Guo, Jingqun Liang, Ping Zhang
2019, 51(3):  787-802.  DOI: 10.6052/0459-1879-18-348
Asbtract ( 288 )   HTML ( 7)   PDF (37993KB) ( 173 )
Figures and Tables | References | Related Articles | Metrics

Flow-induced motions of three elastically mounted circular cylinders subjected to the planar shear flow in tandem arrangement at $Re =100$ were numerically investigated by using semi-implicit characteristics-based split (CBS) finite element algorithm. Firstly, the results are compared with the existing data in the literature to verify the accurateness of the method. Then, the influence of some key parameters, such as shear ratio ($k$), natural frequency ratio ($r$) and reduced velocity ($U_{\rm r}$) on the characteristics of flow-induced dynamic responses and flow field of three circular cylinders with tandem arrangement was analyzed. The numerical results show that the shear ratio, natural frequency ratio and reduced velocity play an important role in the amplitude and trajectory of motion. The change of maximum amplitude of the upstream cylinder is similar to that of an isolated cylinder case with the increasing of $k$. The maximum amplitudes of the midstream and downstream cylinders increase, and the dual resonance is observed. Meanwhile, the increasing of $k$ causes the extension of the resonance region. With the increasing of $r$, the lock-in region of the upstream cylinder reduces in the in-line direction. However, for the midstream and downstream ones, the dual-lock-in region broadens. On the other hand, the prominent $X-Y$ trajectories of three circular cylinder resemble figure-eight shape and irregular one in the uniform flow. With the change of $k$, the motion orbit transforms from figure-eight to the raindrop shape in the lock-in region. In the larger shear ratio and higher natural frequency ratio case, the dual-raindrop is observed in the $X-Y$ vibrating trajectory of the midstream cylinder. Finally, the flow-induced motion mechanism of three tandem circular cylinders exposed to planar shear flow is revealed according to the analysis of the characteristic of the flow field.

NUMERICAL INVESTIGATION ON CAVITY STRUCTURES AND HYRODYNAMICS OF THE VEHICLE DURING VERTICAL WATER-ENTRY1)
Jiayue Zhang, Daqin Li, Qin Wu, Biao Huang, Ying Liu
2019, 51(3):  803-812.  DOI: 10.6052/0459-1879-18-364
Asbtract ( 412 )   HTML ( 4)   PDF (17227KB) ( 124 )
Figures and Tables | References | Related Articles | Metrics

The study of the water-entry behavior of the vehicle with the tail downward is of great significance to many engineering problems, such as recycle of unpowered vehicles and missiles. In this paper, the VOF homogeneous flow model is combined with dynamic mesh technique to study the vertical water-entry process of the vehicle. Good agreement has been obtained between the experimental and numerical results on the water-entry velocity and trajectory. The evolution of the hydrodynamic characteristics, the cavity patterns and the flow structures during the vertical water entry process is analyzed. The results show that the whole water entry process can be divided into four stages: the moment of contact, the open cavity stage, the surface seal stage and the deep seal stage. The pressure drag plays a major role during the water entering process, and the drag coefficient reaches to the maximum in the moment of contact stage when the vehicle touches the free surface. With time evolution, the drag coefficient is gradually decreasing, and tends to be stable in final. A slight fluctuation occurs when the cavity is collapsing. The influence of water-entry velocities on the hydrodynamics and cavity patterns is also studied. During the vertical water-entry of the vehicle. With the increase of the water-entry velocity, the peak of the drag coefficient increases and the maximum dimensionless cavity diameter and the cavity shrinking rate increase. Moreover, the dimensionless moments of the surface seal and deep seal during the vertical water-entry process are almost the same with the different Froude number, as well as the ratio of the depth of the cavity pinch-off to the depth of the vehicle in the deep seal.

VERIFICATION AND VALIDATION OF HYPERFLOW SOLVER FOR SUBSONIC AND TRANSONIC TURBULENT FLOW SIMULATIONS ON UNSTRUCTURED/HYBRID GRIDS1)
Nianhua Wang, Xinghua Chang, Rong Ma, Laiping Zhang
2019, 51(3):  813-825.  DOI: 10.6052/0459-1879-18-331
Asbtract ( 412 )   HTML ( 2)   PDF (2725KB) ( 130 )
Figures and Tables | References | Related Articles | Metrics

Computational fluid dynamics (CFD) is playing a more and more important role in aerospace and relevant industries. While the credibility of CFD numerical simulation requires continuous verification and validation. This paper proposes for CFD solvers a complete verification and validation procedure including rigorous accuracy tests based on the method of manufactured solutions (MMS), simple and complex turbulent flow simulation and grid convergence study. The procedure is implemented on the in-house CFD solver HyperFLOW on verifying and validating its ability on subsonic and transonic turbulent flow simulation. We firstly verified that HyperFLOW can achieve designed second-order accuracy on arbitrary unstructured grids via numerical accuracy tests based on the MMS of the Euler equation and the diffusion equation. Then, simple subsonic turbulent flow cases, such as the turbulent flat plate, the 2D Airfoil Near-Wake and the 2D Bump cases from NASA Turbulence Modeling Resources are taken into consideration for grid convergence tests. Observed accuracy order and grid convergence index are calculated and compared with the results obtained by CFL3D and FUN3D, which verified and validated the accuracy and grid convergence performance of HyperFLOW in simple case simulation. Finally, the NASA Common Research Model in Drag Prediction Workshop VI is simulated with fixed lift coefficient, and grid convergence study is carried out on a series of refined grids. Drag polar is predicted by HyperFLOW and compared with other verified and validated results. It demonstrated that HyperFLOW has good accuracy and grid convergence performance in complex configuration high Reynolds subsonic and transonic turbulent flow simulation on unstructured/hybrid grids.

DIMPLE'S VORTEX ANALYSIS BASED ON EIGENVALUE OF VELOCITY GRADIENT TENSOR1)
Jing Liu, Jie Li, Heng Zhang
2019, 51(3):  826-834.  DOI: 10.6052/0459-1879-18-257
Asbtract ( 434 )   HTML ( 3)   PDF (8374KB) ( 77 )
Figures and Tables | References | Related Articles | Metrics

As a new internal cooling technology for high performance turbine blades, dimple which belongs to vortex generator has the advantage of small flow resistance and high comprehensive heat transfer performance. The qualitative and quantitative analysis of the vortex is an important basis for the optimization design of the dimple, due to the fact that the stronger vortex cause greater heat transfer enhancement. Aiming at the problem that the vortex strength can not be quantitatively analyzed in different dimple models under the variation of the vortex structure, separation mode and background pressure, this paper proposes a method by using vortex core velocity and vortex core velocity gradient tensor eigenvalue. Using the velocity vector and velocity gradient tensor eigenvalue expressed by the local coordinate system at the vortex core, the axial velocity, radial velocity, swirling strength, axial acceleration and radial acceleration of the vortex core are obtained, and quantitative analysis method of vortex strength as a combination of maximum axial velocity, maximum axial acceleration, and maximum swirling strength is concluded. The vortex structure induced by different depth to print diameter ratio dimple is analyzed by this method. As the increase of the depth to print diameter ratio, the maximum axial velocity, maximum axial acceleration and maximum swirling strength all increased obviously, so the vortex strength increases. The research shows that this method has the advantages of simple data processing, wide application range, not limited by separation mode, not limited by background pressure, and the extracted data has clear physical meaning, thus this method is suitable for all kinds of vortex quantitative analysis.

Solid Mechanics
NUMERICAL SIMULATION OF MICRO-ABLATION BEHAVIOR FOR CARBON FIBER REINFORCED COMPOSITES1)
Wei Li, Guodong Fang, Weijie Li, Bing Wang, Jun Liang
2019, 51(3):  835-844.  DOI: 10.6052/0459-1879-18-351
Asbtract ( 491 )   HTML ( 12)   PDF (22787KB) ( 550 )
Figures and Tables | References | Related Articles | Metrics

Carbon fiber reinforced composites are widely used in ablative thermal protection systems (TPS). The microstructures of the composites are greatly related with the ablation behavior. Thus, the study of micro-ablation mechanisms of the composites is significant for the material design and manufacturing. A mathematic micro-ablation model is developed by using the finite-volume-method (FVM) combing with the piecewise linear interface calculation (PLIC) method. Comparing with the analytical results, the numerical method is validated by calculating the ablation morphologies for single fiber embedded in carbon matrix. The effect of carbon fiber inclination on the microscopic ablation behavior is studied for the unidirectional carbon reinforced composites with different carbon fiber inclination. The results are found that the ablation morphology for carbon fiber is bamboo shoot shape if the oxygen diffusion rate is far greater than that of carbon oxygen reaction when the oxidation resistance of carbon fibers is stronger than that of matrix. If the carbon oxygen reaction rate is far greater than that of oxygen diffusion, carbon fiber and matrix will be ablated at the same rate. When the oxidation resistance of carbon fiber is stronger than that of the matrix, the inclination angle of the carbon fiber has a great influence on the micro-ablation behavior of the material. On the contrary, the effect is not significant.

STUDY ON THE INFLUENCE OF COMPRESSIVE STRESS ON THE COMPRESSION SHEAR CRACK PROPAGATION1)
Bo Wang, Yangbo Zhang, Hong Zuo, Houjin Wang
2019, 51(3):  845-851.  DOI: 10.6052/0459-1879-18-369
Asbtract ( 554 )   HTML ( 11)   PDF (3565KB) ( 244 )
Figures and Tables | References | Related Articles | Metrics

MULTIPHASE MATERIAL LAYOUT OPTIMIZATION CONSIDERING EMBEDDING MOVABLE HOLES1)
Xuan Wang, Ping Hu, Kai Long
2019, 51(3):  852-862.  DOI: 10.6052/0459-1879-18-327
Asbtract ( 480 )   HTML ( 14)   PDF (10063KB) ( 97 )
Figures and Tables | References | Related Articles | Metrics

In structural engineering design, it is often necessary to embed one or more fixed-shaped holes to meet certain functional or manufacturing design requirements. To effectively solve the multi-phase material layout optimization problem of continuum structure with embedded movable holes, it is usually necessary to simultaneously optimize the position and orientation of these embedded holes and the topology configuration of the multi-phase material structure to improve the overall performance of the structure. To this end, parameterized level set functions are used to describe the geometry of the embedded holes. The material densities defining the structural topology of multiphase materials, and the geometric parameters used to describe the position and orientation of the embedded holes, are considered as design variables of the optimization problem considered here. To avoid the cumbersome of re-meshing the grids caused by the movement of holes and improve the efficiency of computation, the embedded holes are mapped into a density field on a fixed grid using a smoothed Heaviside function. Meanwhile, a SIMP-like material interpolation invoked at the finite element level is introduced for material parameterization of the optimization problem, and then the simultaneous optimization of the topology configuration of the multi-phase material structure and the position and orientation of the embedded hole can be realized. The material interpolation scheme supports full analytical sensitivity analysis, which allows the current optimization problem to be solved using gradient-based optimization algorithms. Numerical examples illustrate that the proposed method can effectively deal with the layout optimization problem of multiphase material embedded with multiple embedded holes.

A MULTIAXIAL LOW-CYCLE FATIGUE MODEL CONSIDERING NON-PROPORTIONAL ADDITIONAL DAMAGE1)
Xiangyang Cui, Kecheng Hong
2019, 51(3):  863-872.  DOI: 10.6052/0459-1879-18-347
Asbtract ( 361 )   HTML ( 2)   PDF (965KB) ( 160 )
Figures and Tables | References | Related Articles | Metrics

INVESTIGATION ON CRUSHING BEHAVIOR OF METAL HONEYCOMB-LIKE HIERARCHICAL STRUCTURES1)
Yuanming Xia, Wei Zhang, Tianning Cui, Jianxun Zhang, Binwen Wang, Xiaochuan Liu, Chunyu Bai, Qinghua Qin
2019, 51(3):  873-883.  DOI: 10.6052/0459-1879-18-434
Asbtract ( 379 )   HTML ( 9)   PDF (32952KB) ( 164 )
Figures and Tables | References | Related Articles | Metrics

Crushing behavior of metal honeycomb-like hierarchical structures of perforations on single walls (HHSPSW) and double walls (HHSPDW) was systematically investigated by using experiment and numerical simulation methods. Effects of specimen size, perforation location, perforation offsets and perforation gradient on the mechanical properties of honeycomb-like hierarchical structures were analyzed. The results show that the crushing process of honeycomb-like hierarchical structures can be divided into three deformation stages: elastic deformation, buckling deformation and densification. The deformation mode of honeycomb-like hierarchical structures of perforations on single walls is a progressive concave mode, while honeycomb-like hierarchical structures of perforations on double walls deform by axial crushing mode. Specimen sizes have significant influences on mechanical behavior of honeycomb-like hierarchical structures. Mechanical properties are almost independent of honeycomb numbers attaining a certain number. The peak stress of the honeycomb-like hierarchical structures of perforations on single walls is greater than that of honeycomb-like hierarchical structures of perforations on double walls, while its mean crushing stress is less than that of honeycomb-like hierarchical structures of perforations on double walls. Comparing to the traditional honeycombs, the design of perforations on honeycomb walls reduces the specific energy absorption of the honeycombs. Perforation offsets result in decreasing the peak stress of the honeycomb-like hierarchical structures of perforations on single walls, while the mean crushing stress firstly decreases and then increases with increasing the perforation offset. Honeycomb-like hierarchical structures with positive gradient perforations reduce the peak stresses while improve the mean crushing stress comparing to honeycomb-like hierarchical structures with the uniform perforations. The design of muti-gradient perforation distributions on honeycomb walls has small influence on the peak stress and the mean crushing stress.

TOPOLOGY OPTIMIZATION ANALYSIS OF ADHESIVE SOUND ABSORBING MATERIALS STRUCTURE WITH SUBDIVISION SURFACE BOUNDARY ELEMENT METHOD
Leilei Chen, Chuang Lu, Yanming Xu, Wenchang Zhao, Haibo Chen
2019, 51(3):  884-893.  DOI: 10.6052/0459-1879-18-354
Asbtract ( 473 )   HTML ( 10)   PDF (8748KB) ( 119 )
Figures and Tables | References | Related Articles | Metrics

Isogeometric analysis (IGA) realizes the seamless integration of CAD and CAE by using spline basis functions to represent geometric models and implement the numerical analysis, and is widely used in elastic mechanics, electromagnetic field, potential problem and other fields. However, it is difficult to construct directly a complex structure by using IGA. The subdivision surface method can be used to subdivide and fit the discrete mesh of a traditional model in order to construct smooth surfaces. Such a method is suitable for complicated problems. This method has the following advantages: (1) It is suitable for any topological structure; (2) The numerical calculation is stable; (3) It is simple to implement; (4) Local refinement and continuity control. Because of its flexibility and convenience in the construction of complex structural models, this method has been widely used in aerospace, automobile, animation, game making and other modeling fields. The subdivision surface method (SSM) and the boundary element method (BEM) were integrated to perform structural acoustic analysis. The box-spline interpolation of SSM was used for both geometric interpolation and physical interpolation, achieving high-order approximation of the structural surface and physical field. The acoustic scattering of the structure of adhesive sound absorption materials was taken as an example to test the effectiveness of the algorithm. The above analysis was combined with the method of moving asymptotes (MMA) algorithm to conduct a topological optimization of the distribution of sound absorption materials. In this study, the adjoint variable method and the acoustic BEM were used to analyze the sensitivity of the distribution of sound absorption materials on the surface of the structure. Each update of design variables brings small changes in the layout of sound absorbing materials, and ultimately the optimal solution is obtained. The resulting solvers provide an efficient computational tool for topology optimization design. The proposed algorithm is then applied to some numerical examples to illustrate the potential for engineering optimization design.

Dynamics, Vibration and Control
PARAMETERS OPTIMIZATION OF A DYNAMIC VIBRATION ABSORBER WITH AMPLIFYING MECHANISM AND NEGATIVE STIFFNESS1)
Zhaoyang Xing, Yongjun Shen, Haijun Xing, Shaopu Yang
2019, 51(3):  894-903.  DOI: 10.6052/0459-1879-18-375
Asbtract ( 477 )   HTML ( 8)   PDF (1543KB) ( 198 )
Figures and Tables | References | Related Articles | Metrics

Mechanical vibration is detrimental in most engineering situations, and it may not only generate noise but also reduce the operational accuracy and working life of the equipment. It is generally difficult for vibration absorption and isolation systems with positive stiffness characteristics to achieve satisfactory performance, which becomes noticeable especially in low-frequency vibration control systems. Amplifying mechanism and negative stiffness element both show good performance in the field of vibration control, but the dynamic vibration absorber with both amplifying mechanism and negative stiffness element is rarely studied. Based on the Voigt type dynamic vibration absorber, a dynamic vibration absorber with negative stiffness element using amplifying mechanism is presented, and the optimal system parameters are studied in detail. Firstly, the differential equation of motion is established and the analytic solution of the system is obtained, and it is found that there are two fixed points independent of damping ratio in the amplitude-frequency curves of the primary system. The optimal frequency ratio of the dynamic vibration absorber is obtained based on the fixed-point theory. According to the characteristics of negative stiffness, the optimal negative stiffness ratio is founded under the premise of ensuring the system stability. A simple method is used to derive the approximate optimal damping ratio of the system. The correctness of the analytical results is verified by the comparison with the results by numerical simulation. Compared with other dynamic vibration absorbers under harmonic and random excitations, it could be found that the model with optimal parameters in this paper can greatly reduce the resonance amplitude, broaden the vibration band, and lower the resonance frequency of the primary system. These results may provide theoretical basis for the optimal design of similar dynamic vibration absorbers.

POSITIVE AND NEGATIVE PULSE-SHAPED EXPLOSION AS WELL AS BURSTING OSCILLATIONS INDUCED BY IT1)
Mengke Wei, Xiujing Han, Xiaofang Zhang, Qinsheng Bi
2019, 51(3):  904-911.  DOI: 10.6052/0459-1879-18-319
Asbtract ( 495 )   HTML ( 3)   PDF (1282KB) ( 111 )
Figures and Tables | References | Related Articles | Metrics

Bursting oscillations, characterized by the alternation between large-amplitude oscillations and small-amplitude oscillations, are complex behaviors of dynamical systems with multiple time scales and have become one of the hot subjects in nonlinear science. Up to now, various underlying mechanisms of bursting oscillations as well as the classifications have been investigated intensively. Recently, a sharp transition behavior, called the "pulsed-shaped explosion (PSE)", was uncovered based on nonlinear oscillators of Rayleigh's type. PSE is characterized by pulse-shaped sharp quantitative changes appearing in the branches of equilibrium point and limit cycle. However, the previous work related to the PSE merely focused on the sharp transitions of unidirectional PSE, and more complex forms of PSE which may lead to more complicated bursting patterns need to be further investigated. Taking a parametrically and externally excited Rayleigh system as an example, we reveal different expression of PSE as well as bursting patterns induced by it. According to the frequency relationship between the two slow excitations, the fast subsystem and the slow variable are obtained by means of frequency-transformation fast-slow analysis. Bifurcation behaviors of the fast subsystem show that, there exist two originally disunited branches of equilibrium point and limit cycle, which extend steeply in different directions and inosculate as a integral structure according to the variation of system parameters. This inosculated integral structure inherits the "steep" properties in different directions from the originally disunited branches, and PSE is thus generated. Unlike the PSE phenomena studied in the previous works, the PSE reported here contains two different peaks in positive and negative directions, which can be named as "positive and negative PSE". Note that the positive and negative PSE essentially complicates bursting dynamics and plays a critical role in bursting. On the other hand, only with the properly chosen parameters could it be created. Based on this, two different types of bursting, i.e., point-point type and cycle-cycle type, are obtained. Subsequently, the transformed phase portraits are introduced to explore dynamical mechanisms of the bursting patterns. We show that, with the variation of the slow parameter, the trajectory may undergo sharp transitions between the rest and active states by positive and negative PSE, and therein lies the generation of bursting. Our results demonstrate the diversity of PSE and give a complement to the underlying mechanisms of bursting.

NONLINEAR VIBRATIONS OF COMPOSITE CANTILEVER PLATE IN SUBSONIC AIR FLOW1)
Gen Liu, Wei Zhang
2019, 51(3):  912-921.  DOI: 10.6052/0459-1879-18-336
Asbtract ( 358 )   HTML ( 6)   PDF (1182KB) ( 248 )
Figures and Tables | References | Related Articles | Metrics

With the development of materials science, more and more new materials have been applied to engineering practice. Under the action of airflow excitation, the nonlinear dynamics of plate and shell structures with composite materials based on aerospace engineering is still a hot research topic. In this work, the nonlinear vibrations and responses of a laminated composite cantilever plate under subsonic air flow are investigated. According to the flow condition of ideal incompressible fluid and the Kutta--Joukowski lift theorem, the subsonic aerodynamic lift on the three-dimensional finite length flat wing is calculated by using the Vertex Lattice (VL) method, which is applied to the cantilever plate. The finite length flat wing is modeled as a laminated composite cantilever plate based on the Reddy's third-order shear deformation plate theory, moreover the von Karman geometry nonlinearity is introduced. The nonlinear partial differential governing equations of motion for the laminated composite cantilever plate subjected to the subsonic aerodynamic force are established via Hamilton's principle. The partial differential equations are separated into two nonlinear ordinary differential equations via Galerkin method. Through comparing the natural frequencies of the system with different material and geometry parameters, the 1:2 internal resonance is considered here using multiple scales method. Corresponding to several selected parameters, the frequency-response characteristics are obtained. The hardening-spring-type behaviors and jump phenomena are exhibited.

DIMENSION REDUCTION OF FPK EQUATION FOR VELOCITY RESPONSE ANALYSIS OF STRUCTURES SUBJECTED TO ADDITIVE NONSTATIONARY EXCITATIONS1)
Zhenmei Rui, Jianbing Chen
2019, 51(3):  922-931.  DOI: 10.6052/0459-1879-18-411
Asbtract ( 369 )   HTML ( 5)   PDF (8144KB) ( 93 )
Figures and Tables | References | Related Articles | Metrics

The nonlinear response analysis of structures subjected to random excitation is a highly challenging problem. The solution of FPK equation provides exact solutions to these problems. However, for nonlinear multi-degree-of-freedom systems, the direct solutions of the FPK equation is prohibitively difficult. Actually, the numerical solutions are strictly limited by the dimension of the equation, while the analytical solutions are available only for very few specific systems, and most of them are steady-state solution. Therefore, reducing the dimension of the FPK equation is an important way to solve the high-dimensional nonlinear dynamic response analysis problem. In the present paper, for the nonstationary response analysis of multi-degree-of-freedom nonlinear structures subjected to amplitude-modulated additive white noise, the high-dimensional FPK equation in terms of the joint probability density function is reduced in dimension. For the probability density function of velocity response, it is converted to a one-dimensional FPK-like equation through introducing the equivalent drift coefficient. The method of conditional mean function estimate is suggested to construct the equivalent drift coefficient. Afterwards, the numerical results of probability density function of velocity can be obtained by applying the path integration method to solve the dimension-reduced FPK-like equation. The accuracy and efficiency of the proposed method are discussed and verified through the numerical examples, including the non-stationary response analysis of velocity of a single-degree-of-freedom Rayleigh oscillator, a ten-story linear shear frame structure and a nonlinear shear frame structure subjected to amplitude-modulated additive white noise.

Biomechanics, Engineering and Interdiscipliary Mechanics
STUDY ON MOISTURE TRANSPORT CHARACTERISTICS OF SHALE BASED ON ISOTHERMAL ADSORPTION1)
Weijun Shen, Xizhe Li, Xiaobing Lu, Yujin Wan, Wei Guo, Luo Zuo
2019, 51(3):  932-939.  DOI: 10.6052/0459-1879-18-229
Asbtract ( 494 )   HTML ( 59)   PDF (2515KB) ( 240 )
Figures and Tables | References | Related Articles | Metrics

The study on moisture transport characteristics of shales is critical, which is not only helpful to understand the physical and chemical properties in shales, but also to evaluate the adsorption, diffusion and flow ability of shale gas. In this study, the experimental device of moisture transport in shales was designed and the shale samples from Woodford in USA and Longmaxi Formation in Southern China were used. The moisture transport in shales was carried out at different temperatures and humidities, and the transport characteristics and the effects in shales were investigated. The results indicate that moisture adsorption isotherms of shales belong to type II curve, including the monolayer, multilayer adsorption and capillary condensation, and the GAB model can be used to describe the moisture adsorption process of shale rocks. With the increasing of relative pressure, the moisture adsorption of shales increases. The content of organic carbon and temperature strengthen the moisture adsorption in shales while calcite will inhibit the process. The moisture diffusion coefficient in shales initially increases, then decreases and finally increases with relative pressure, and the value ranges between 8.73$\times$10$^{ - 9}$ m$^{2}$/s and 5.95$\times$10$^{ - 8 }$ m$^{2}$/s. The isothermal heat of moisture adsorption in Woodford shale is higher than that of Longmaxi Formation, which is related to shale maturity. These results provide some reference basis for understanding the physical and chemical properties in shales and evaluating the adsorption and flow capacity of shale gas.

STUDY ON STRENGTH CHARACTERISTICS OF MICROPOROUS CLAY IN SHALE BASED ON HOMOGENIZATION THEORY1)
Qiang Han, Zhan Qu, Zhengyin Ye
2019, 51(3):  940-948.  DOI: 10.6052/0459-1879-18-214
Asbtract ( 333 )   HTML ( 9)   PDF (301KB) ( 97 )
Figures and Tables | References | Related Articles | Metrics

As one of the basic parameters necessary for shale oil and development, the analysis of shale strength is carried out in the whole process of drilling and hydraulic fracturing. Macroscopic experiments have problems such as sample preparation and time consuming. Limited by downhole conditions, not only the quality of data obtained by geophysical method is not good enough for mechanical analysis, but also it increases the risk of equipment stuck and buried in downhole. In this paper, the strength evaluation method of microporous clay in shale was proposed based on the homogenization theory. The composition and mechanical analysis of porous clay was carried out. Based on dissipative energy principle and Drucker-Prager criterion, the strength evaluation of porous clay was transformed into a solution to the strain of the microscopic $\pi$ function. The mechanical properties of the intergranular pores of clay were discussed and the homogenization strain energy of porous clay was established. The microscopic nonlinear function was constructed based on the strength homogenization theory. A homogenization $\pi$ function was established in relation to parameters such as the composition of porous clay, coefficient of friction and cohesion. Based on nanomechanical experiments, dimensional analysis and finite element simulation, the intrinsic relationship between hardness, strength and composition of porous clay was evaluated. The results show that the elastic modulus and hardness of microporous clay in shale are positively correlated with the packing density of shale. The ratio of hardness to cohesion coefficient exhibits a nonlinear increase with increasing friction coefficient when the clay packing density is constant. The $\pi$ function of porous clay with respect to hardness, strength and clay packing density is solved by dimensional analysis and finite element simulation. The composition and mechanical relationship of shale microporous clay are described. It lays a foundation for further research on shale meso-strength parameters and macro-strength prediction.

THE ACOUSTIC CHARACTERISTICS OF DOLOMITE UNDER COMPRESSION AND ITS APPLICATION IN THE CRACK RESEARCH1)
Honglin Zhu, Qiao Chen, Xiaohu Xu, Hong Liu, Xiangyu Wen, Jilong Chen, Wei Peng, Binling Zhao
2019, 51(3):  949-960.  DOI: 10.6052/0459-1879-18-328
Asbtract ( 416 )   HTML ( 4)   PDF (11430KB) ( 54 )
Figures and Tables | References | Related Articles | Metrics

The variation characteristics of acoustic velocity, amplitude and spectral characteristics of dolomite under triaxial and uniaxial loading are studied by means of the experimental method. The following conclusions can be drawn. (1) During the process of compression deformation of the dolomite, the change of acoustic wave velocity well reflects the closure, generation, expansion and penetration of cracks in the rock. In general, the shear acoustic wave velocity is better than the longitudinal acoustis wave to predict the generation of cracks, while the longitudinal acoustic wave has a more sensitive reflection on the unstable expansion of the cracks. (2) With the axial load increases, the scattered wave signal appears at the end of the longitudinal acoustic wave and transverse wave pattern. And the "fishtail" scattered wave signal at the end of the shear waveform is more obvious (at this point, the axial stress reaches 60% of the limit strength), which indicates the generation and stability expansion of the crack in the rock. (3) The spectral curve is also a good reflection to the deformation of the internal structure in the dolomite rock. As the axial load increases, the amplitude of the spectrum curve shows an increasing trend, which indicates the compaction stage of the rock; when the low frequency signal is more active than the high frequency one in frequency spectrum curve, it marks the occurrence of the crack; what is more, low frequency will replace the high frequency to become the main frequency. (4) In the stage of crack closure, both the first arrival wave amplitude and the principal amplitude tend to increase; while in the stage of unsteady expansion of the crack, the principal amplitude tends to increase more slowly than the first arrival wave amplitude; after the amplitude curve reaches the peak there is a sudden drop point, which indicates the destruction of the rock. The variable characteristic of acoustic propagation displayed during the process of rock pressed makes a good sense to the prediction of rock crack and the evaluation of rock stability. This study has important theoretical reference value for the dynamic long-term monitoring to the rock deformation and the stability evaluation of engineering rock mass.

Science Foundation
A BRIEF INTRODUCTION OF COMPLETED KEY PROGRAM PROJECTS ON MECHANICS IN 2018
Kuncao Bai, Dongxing Cao, Gang Wang, Tiangang Lei
2019, 51(3):  965-969.  DOI: 10.6052/0459-1879-19-097
Asbtract ( 444 )   HTML ( 13)   PDF (171KB) ( 207 )
Figures and Tables | References | Related Articles | Metrics

The article briefly introduces the completion and evaluation of 12 NSFC key program projects on mechanics in 2018. The project list and the evaluation opinions of expert committee are given in detail.