Prof. Chenfeng Li FLSW

Director, Zienkiewicz Institute for Modelling, Data and AI Post: Faculty of Science and Engineering, Swansea University Bay Campus, Swansea SA1 8EN, UK.
Non-executive Director, Temporary Works Forum Tel: +44-01792-602256
Editor-in-Chief, Engineering Computations Email:


  • Computational Engineering and Data Analytics
  • Uncertainty Quantification, Reliability and Risk Assessment
  • Civil, Construction and Geotechnical Engineering
  • Concrete, Composite, Heterogeneous and Porous Materials

  • Characterization and modelling of heterogeneous and porous materials
  • Numerical simulation of multiphase and non-Newtonian fluids
  • Uncertainty quantification, reliability and risk assessment
  • Real-time simulation of deformation, fracture, energy and mass transport

  • JOURNAL Papers


    An efficient 3-node triangular plate element for static and dynamic analyses of microplates based on modified couple stress theory with micro-inertia
    Engineering with Computers, 39 (5) (2023) 3061-3084.
    Y.H. Mao, Y. Shang, S. Cen and C.F. Li

    Within the modified couple stress elasticity, a novel nonconforming 3-node triangular plate element is derived from the d'Alembert–Lagrange principle for simulating the size-dependent static and dynamic bending behaviors of thin microplates and studying the effects of micro-inertia. This is accomplished via two main steps. First, the Trefftz functions that are derived from the governing equations of the problem concerned are adopted as the basis functions for constructing the element's displacement interpolation. Second, according to the generalized conforming theory, the SemiLoof constraints are used to enforce the C1 compatibility requirement for guaranteeing the computation convergence. The benchmark tests are carried out and the results reveal that the new element exhibits satisfactory numerical accuracy and captures the size dependences effectively in the static, free vibration and forced vibration analyses. Moreover, the findings also show that the micro-inertia affects the dynamic response of the plate mainly through the natural frequency. In general, the influences on higher order modes are more obvious than that on lower order modes.

    Penalty C0 8-node quadrilateral and 20-node hexahedral elements for consistent couple stress elasticity based on the unsymmetric finite element method
    Engineering Analysis with Boundary Elements, 147 (2023) 302-319.
    H.P. Wu, Y. Shang, S. Cen and C.F. Li

    In this paper, the penalty unsymmetric finite element framework for the consistent couple stress theory is derived from the virtual work principle. The C1 continuity requirement is satisfied in weak form by using the penalty function method to constrain the independently introduced rotations for approximating the mechanical rotations, enabling the utilization of C? continuous interpolations for designing the element displacement without the loss of convergence property. Within the proposed framework, 8-node quadrilateral element and 20-node hexahedral solid element are constructed for analyzing the size-dependent mechanical responses of consistent couple stress elasticity materials. In these developments, the quadratic serendipity isoparametric shape functions are enriched by the rotation degrees of freedom for determining the test functions, whilst the metric stress functions that are derived from the concerned equilibrium equations are used to design the trial functions. A series of numerical benchmarks are examined for verifying their effectiveness and accuracy. It is shown that the elements can efficiently capture the size dependences, exhibiting good accuracy and low susceptibility to mesh distortion.

    High-performance unsymmetric 8-node hexahedral element in modeling nearly-incompressible soft tissues
    International Journal of Mechanical Sciences, 260 (2023) 108647.
    Y.F. Wang, S. Cen, C.F. Li and Q. Zhang

    In biomechanics problems, the biological soft tissues are usually treated as anisotropic nearly incompressible hyperelastic materials, but such complicated nonlinear material models often cause challenging problems of severe volumetric locking and instabilities in numerical simulations. In this paper, the recent unsymmetric 8-node, 24-DOF hexahedral solid element US-ATFH8 with different test and analytical trial functions (ATFs) is modified for the analysis of the anisotropic nearly-incompressible hyperelastic soft tissues. Unlike the original formulation, the linear analytical general solutions for anisotropic elasticity and the consistent tangent modulus are firstly employed for formulating the trial functions, and are used to construct the incremental displacement fields that result in the incremental deformation gradient. The total deformation gradient is obtained by multiplying the incremental deformation gradient by the deformation gradient, after which the Cauchy stresses can be directly calculated from a total-form constitutive equation relating to the deformation gradient. Numerical tests, including commonly used benchmarks and cardiac examples, demonstrate attractive properties of the proposed finite element formulation in modeling nearly-incompressible anisotropic hyperelastic materials. It is free of various locking and quite insensitive to mesh distortions, and provides high accuracy with faster convergence rates when compared with other existing models.

    Physics-data combined machine learning for parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes
    Computer Methods in Applied Mechanics and Engineering, 404 (2023) 115771.
    J.L. Fu, D.H. Xiao, R. Fu, C.F. Li, C.H. Zhu, R. Arcucci, I.M. Navon

    Repeatedly solving nonlinear partial differential equations with varying parameters is often an essential requirement to characterise the parametric dependences of dynamical systems. Reduced-order modelling (ROM) provides an economical way to construct low-dimensional parametric surrogates for rapid predictions of high-dimensional physical fields. This paper presents a physics-data combined machine learning (PDCML) method for non-intrusive parametric ROM in small-data regimes. Proper orthogonal decomposition (POD) is adopted for dimension reduction by deriving basis functions from a limited number of highfidelity snapshots, and parametric ROM is thus transformed into establishing reliable mappings between the system parameters and the POD coefficients. To overcome labelled data scarcity, a physics-data combined ROM framework is developed to jointly integrate the physical principle and the small labelled data into feedforward neural networks (FNN) via a step-by-step training scheme. Specifically, a preliminary FNN model is firstly fitted via data-driven training, and then the governing physical rules are embedded into the loss function to improve the model interpolation and extrapolation performances through physics-guided training constrained by the labelled data. During the constrained optimisation procedure, dynamic weighting factors are used to adjust the physics-data proportion of the loss functions, aiming at continuously highlighting the physics loss as the primary optimisation objective and keeping the data loss as the constraint. This new PDCML method is tested on a series of nonlinear problems with different numbers of physical variables, and it is also compared with the data-driven ROM, the physics-guided ROM and the traditional projection-based ROM methods. The results demonstrate that the proposed method provides a costeffective way for non-intrusive parametric ROM via machine learning, and it possesses good characteristics of high prediction accuracy, strong generalisation capability and small data requirement. ©2022 Elsevier B.V. All rights reserved.


    Extension of the unsymmetric 8-node hexahedral solid element US-ATFH8 to 3D hyper-elastic finite deformation analysis
    International Journal for Numerical Methods in Engineering, 123 (23) (2022) 5749-5778.
    R.X. Ma, S. Cen, Y. Shang and C.F. Li

    This work extends the recent US-ATFH8 element to 3D hyper-elastic finite deformation analysis. Using two sets of shape functions, the new 3D element comprises of 8 nodes and 24 DOFs. The first set of shape functions represent the test functions that come from the conventional isoparametric interpolation, and the second set, representing the trial functions, are constructed from the homogenous solutions for linear elasticity governing equations, termed analytical trial functions (ATFs). This study considers finite deformation for hyper-elastic materials, but it is assumed that the analytical solutions associated with hyper-elastic materials can be updated to hold approximately in each incremental step. Moreover, the deformation information required for stress computation is updated by using the incremental deformation gradient, which is constructed from the updated ATFs. Numerical examples show that without additional pressure DOF, the element US-ATFH8 still behaves well in nearly incompressible hyper-elastic 3D problems with finite deformation, even when the meshes are extremely distorted.

    Investigation on energy conversion instability of pump mode in hydro-pneumatic energy storage system
    Journal of Energy Storage, 53 (2022) 105079.
    C.Y. Wang, F.J. Wang, C.F. Li, W.H. Chen, H. Wang and L. Lu

    The pump mode of hydro-pneumatic energy storage (HPES) system often experiences off-design conditions due to the boundary pressure rises, and the resultant energy conversion instability has an adverse effect on the system operation. However, the evolutionary process of this instability and the corresponding flow events are still not fully understood. Experimental and numerical simulation studies of a centrifugal pump with vaned diffuser were conducted in a wide operating range, and the following notable results were obtained. The energy conversion of pump mode shows three typical stages as the pump head increases, stable work?stumbled work?shaky work, corresponding to the inner flow pattern development of stall-free?stall inception?stall deepening. In this evolutionary process, the attenuation of rotor blade loading gradually spreads from the inlet shroud side to the entire blade. This strengthens the natural adverse potential rothalpy gradient and induces the energetic vortical structures including leading backflow vortices, dynamic stall cells and trailing helical vortices. A strong coupling between large-scale vortical motions and energy conversion is observed and this yields a Logistic growth of the shaft power and pump head fluctuations. The critical condition for the energy conversion instability is the diagnostic function SI = 1, and this is recommended for determining the stability limit of the pump mode.

    How opening windows and other measures decrease virus concentration in a moving car
    Engineering Computations, 39 (6) (2022) 2350-2366.
    S. Shu, T.E. Mitchell, M.R.R. Wiggins, S.Z. You, H.R. Thomas and C.F. Li

    Purpose – Due to the ongoing Covid-19 pandemic, ventilation in a small cabin where social distancing cannot be guaranteed is extremely important. This study aims to find out the best configuration of open and closed windows in a moving car at varying speeds to improve the ventilation efficiency. The effectiveness of other mitigation measures including face masks, taxi screens and air conditioning (AC) systems are also evaluated. Design/methodology/approach – Each window is given three opening levels: fully open, half open and fully closed. For a car with four windows, this yields 81 different configurations. The location of virus source is also considered, either emitting from the driver or from the rear seat passenger. Then three different travelling speeds, 5 m/s, 10 m/s and 15 m/s, are examined for the window opening/closing configurations that provide the best ventilation effect. A study into the effectiveness of face masks is realised by adjusting virus injection amounts; and the simulation of taxi screens and AC system simply requires a small modification to the car model. Findings – The numerical studies identify the top window opening/closing configurations that provide the most efficient ventilation at different moving speeds, along with a comprehensive ranking list. The results show that fully opening all windows is not always the best choice. Simulations evaluating other mitigation measures confirm good effect of face masks and poor performance of taxi screens and AC systems. Originality/value – This work is the first large-scale numerical simulation and parametric study about different window opening/closing configurations of a moving car. The results provide useful guides for travellers in shared cars to mitigate Covid-19 transmission risks. The findings are helpful to both individuals' health and society's recovery in the Covid-19 era and they also provide useful information to protect people from other respiratory infectious diseases such as influenza.

    An efficient 4-node facet shell element for the modified couple stress elasticity
    International Journal for Numerical Methods in Engineering, 123 (2022) 992-1012.
    Y. Shang, H.P. Wu, S. Cen and C.F. Li

    This work proposes a simple but robust 4-node 24-DOF facet shell element for static analysis of small-scale thin shell structures. To accommodate the size effects, the modified couple stress theory is employed as the theoretical basis. The element is constructed via two main innovations. First, the trial functions that can a priori satisfy related governing differential equations are adopted as the basic functions for formulating the element interpolations. Second, the generalized conforming theory and the penalty function method are employed to meet the C1 continuity requirement inweak sense for ensuring the computation convergence property. Several benchmarks of shells with different geometries are tested to assess the new facet shell element's capability. The numerical results reveal that the element can effectively simulate the size-dependent mechanical behaviors of small-scale thin shells, exhibiting satisfactory numerical accuracy and low susceptibility to mesh distortion. Moreover, as the shell element uses only six degrees of freedom per node, it can be incorporated into the commonly available finite element programs very readily.

    Simulating fractures with bonded discrete element method
    IEEE Transactions on Visualization and Computer Graphics, 28 (12) (2022) 4810-4824.
    J.M. Lu, C.F. Li, G.C. Cao and S.M. Hu

    Along with motion and deformation, fracture is a fundamental behaviour for solid materials, playing a critical role in physically-based animation. Many simulation methods including both continuum and discrete approaches have been used by the graphics community to animate fractures for various materials. However, compared with motion and deformation, fracture remains a challenging task for simulation, because the material's geometry, topology and mechanical states all undergo continuous (and sometimes chaotic) changes as fragmentation develops. Recognizing the discontinuous nature of fragmentation, we propose a discrete approach, namely the Bonded Discrete Element Method (BDEM), for fracture simulation. The research of BDEM in engineering has been growing rapidly in recent years, while its potential in graphics has not been explored.We also introduce several novel changes to BDEM to make it more suitable for animation design. Compared with other fracture simulation methods, the BDEM has some attractive benefits, e.g., efficient handling of multiple fractures, simple formulation and implementation, and good scaling consistency. But it also has some critical weaknesses, e.g., high computational cost, which demand further research. A number of examples are presented to demonstrate the pros and cons, which are then highlighted in the conclusion and discussion.

    A review of hydraulic fracturing simulation
    Archives of Computational Methods in Engineering, 29 (2022) 1-58.
    B. Chen, B.R. Barboza, Y.N. Sun, J. Bai, H.R. Thomas, M. Dutko, M. Cottrell and C.F. Li

    Along with horizontal drilling techniques, multi-stage hydraulic fracturing has improved shale gas production significantly in past decades. In order to understand the mechanism of hydraulic fracturing and improve treatment designs, it is critical to conduct modelling to predict stimulated fractures. In this paper, related physical processes in hydraulic fracturing are firstly discussed and their effects on hydraulic fracturing processes are analysed. Then historical and state of the art numerical models for hydraulic fracturing are reviewed, to highlight the pros and cons of different numerical methods. Next, commercially available software for hydraulic fracturing design are discussed and key features are summarised. Finally, we draw conclusions from the previous discussions in relation to physics, method and applications and provide recommendations for further research.

    Stochastic reconstruction of 3D microstructures from 2D cross-sectional images using machine learning-based characterization
    Computer Methods in Applied Mechanics and Engineering, 390 (2022) 114532.
    J.L. Fu, D.H. Xiao, D.F. Li, H.R. Thomas and C.F. Li

    The availability of high-quality three-dimensional (3D) microstructures is an essential prerequisite to understanding the macroscopic behaviours of heterogeneous media. Considering the high cost of 3D microscopy imaging, it is an economical way to derive large numbers of virtual 3D microstructures from limited morphological information of 2D cross-sectional images. This paper presents an efficient method that incorporates machine learning-based characterization of 2D images to statistically reconstruct 3D microstructures. The latent stochasticity of 2D images is mastered by fitting supervised machine learning models, which essentially characterize the local morphological statistics. A morphology integration scheme is developed to project the 2D morphological statistics into the 3D space, and new equivalent 3D microstructures can then be synthesized by sequentially generating voxel values from probability sampling. The new method is tested on a series of heterogeneous media with distinct morphologies, and it is also compared with two classical methods (i.e. stochastic optimization-based reconstruction and Gaussian random field transformation) in terms of reconstruction accuracy and efficiency. Besides, various microstructural descriptors are used to quantify the discrepancies between the reconstructed and target microstructures. The results confirm that the proposed method provides a highly cost-effective and widely applicable way to reproduce 3D realizations that precisely preserve the statistical characteristics, geometrical irregularities, long-distance connectivity and anisotropy that exist in 2D cross-sectional images.

    A data-driven approach to predicting the anisotropic mechanical behaviour of voided single crystals
    Journal of Mechanics and Physics of Solids, 159 (2022) 104700.
    H.J. Guo, C. Ling, D.F. Li, C.F. Li and E.P. Busso

    Machine learning techniques are increasingly used to extract important physical information from a broad range of materials and to identify their processing-structure–property relationships. In this work, a neural network framework is coupled with a crystal plasticity finite element based approach to predict the deformation and ductile failure behaviour of porous FCC single crystals when subjected to multi-axial loading conditions. The work relies on 3D unit cells with a centrally located spherical void to represent the microstructure of the porous single crystal material, and on a crystallographic slip-based crystal plasticity constitutive model to describe its deformation behaviour. Stress–strain data generated by unit cell finite element simulations are relied upon to construct the neural network model. Different strategies for the neural network input and output variables and parameters are first explored so as to optimise performance and accuracy. Both proportional and non-proportional loading conditions resulting from a constant and a varying stress triaxiality during deformation, respectively, are considered. The optimum neural network strategy is shown to successfully predict the behaviour of the porous single crystal under both proportional and non-proportional loading, albeit with the void behaviour at high stress triaxialities better described than that at low triaxialities. The results also reveal that the use of tensorial quantities for both stresses and strains as input and output neural network quantities is more suitable as a form of data representation for multiaxial loading conditions than uniaxial equivalent stress and strain quantities. It was also found that the inclusion of prior knowledge as neural network input quantities in the form of reference stress–strain solutions for the void-free single crystal considerably improves the predictive capabilities of the proposed data-driven approach, even when only a very limited number of training cases was used.

    Incompressibility enforcement for multiple-fluid SPH using deformation gradient
    IEEE Transactions on Visualization and Computer Graphics, 28 (10) (2022) 3417-3427.
    B. Ren, W. He, C.F. Li and X. Chen

    To maintain incompressibility in SPH fluid simulations is important for visual plausibility. However, it remains an outstanding challenge to enforce incompressibility in such recent multiple-fluid simulators as the mixture-model SPH framework. To tackle this problem, we propose a novel incompressible SPH solver, where the compressibility of fluid is directly measured by the deformation gradient. By disconnecting the incompressibility of fluid from the conditions of constant density and divergence-free velocity, the new incompressible SPH solver is applicable to both single- and multiple-fluid simulations. The proposed algorithm can be readily integrated into existing incompressible SPH frameworks developed for single-fluid, and is fully parallelizable on GPU. Applied to multiple-fluid simulations, the new incompressible SPH scheme significantly improves the visual effects of the mixture-model simulation, and it also allows exploitation for artistic controlling.


    Investigation of hydraulic fracture branching in porous media with a hybrid finite element and peridynamics approach
    Theoretical and Applied Fracture Mechanics, 116 (2021) 103133.
    Y.N. Sun, B. Chen, M.G. Edwards and C.F. Li

    Simulation of complex fracture patterns in porous media can help understand and improve hydraulic fracturing processes, with potential for significant impact on enhancing oil and gas recovery. In this paper, a fully coupled hydraulic fracture propagation simulation method employing a hybrid finite element method (FEM) and peridynamic (PD) approach is presented. Considering the ability of PD in solving discontinuous problems, the area where cracks can potentially occur is discretised by PD and the crack-free area is discretised by FEM. The solid deformation and fracture propagation are captured by PD and FEM, while the fluid flow in both the reservoir and fracture is simulated with FEM. The whole process is solved in a monolithic way with an implicit scheme. The presented method demonstrates the capability of modelling complex dynamic crack propagation via benchmark examples. Branching phenomenon is then investigated with the proposed model. It is found that faster loading rate, lower-energy release rate, and more brittle and impermeable media will cause crack branching more easily.

    Extension of the unsymmetric 8-node hexahedral solid element US-ATFH8 to geometrically nonlinear analysis
    Engineering Computations, 38 (8) (2021) 3219-3253.
    Z. Li, S. Cen and C.F. Li

    Purpose – The purpose of this paper is to extend a recent unsymmetric 8-node, 24-DOF hexahedral solid element US-ATFH8 with high distortion tolerance, which uses the analytical solutions of linear elasticity governing equations as the trial functions (analytical trial function) to geometrically nonlinear analysis. Design/methodology/approach – Based on the assumption that these analytical trial functions can still properly work in each increment step during the nonlinear analysis, the present work concentrates on the construction of incremental nonlinear formulations of the unsymmetric element US-ATFH8 through two different ways: the general updated Lagrangian (UL) approach and the incremental corotational (CR) approach. The key innovation is how to update the stresses containing the linear analytical trial functions. Findings – Several numerical examples for 3D structures show that both resulting nonlinear elements, US-ATFH8-UL and US-ATFH8-CR, perform very well, no matter whether regular or distorted coarse mesh is used, and exhibit much better performances than those conventional symmetric nonlinear solid elements. Originality/value – The success of the extension of element US-ATFH8 to geometrically nonlinear analysis again shows the merits of the unsymmetric finite element method with analytical trial functions, although these functions are the analytical solutions of linear elasticity governing equations.

    A state-of-the-art review of crack branching
    Engineering Fracture Mechanics, 257 (2021) 108036.
    Y.N. Sun, M.G. Edwards, B. Chen and C.F. Li

    Crack branching has important theoretical and practical significance in many natural phenomena and practical engineering problems. At present, the field of crack branching is still at an exploration stage, lacking a unified explanation of the underlying mechanisms and an effective method to predict crack branching in practical materials. This paper provides a state-of-the-art review of crack branching, including experimental observations, physics, fracture models and associated numerical methods. The experimental observations are first summarized, followed by the physics of crack branching. Then, the crack models including discrete crack models and smeared crack models are discussed, highlighting their key features, advantages and limitations. Next, a number of numerical methods that have been used to simulate crack branching are reviewed in detail, including the finite element method (FEM), extended finite element method (XFEM), boundary element method (BEM), meshfree methods (MMs), peridynamics (PD) and discrete element method (DEM). Finally, based on the information reviewed above, the future research directions of crack branching modelling are discussed.

    Efficient structural reliability analysis based on adaptive Bayesian support vector regression
    Computer Methods in Applied Mechanics and Engineering, 387 (2021) 114172.
    J.S. Wang, C.F. Li, G.J. Xu, Y.L. Li and A. Kareem

    To reduce the computational burden for structural reliability analysis involving complex numerical models, many adaptive algorithms based on surrogate models have been developed. Among the various surrogate models, the support vector machine for regression (SVR) which is derived from statistical learning theory has demonstrated superior performance to handle nonlinear problems and to avoid overfitting with excellent generalization. Therefore, to take the advantage of the desirable features of SVR, an Adaptive algorithm based on the Bayesian SVR model (ABSVR) is proposed in this study. In ABSVR, a new learning function is devised for the effective selection of informative sample points following the concept of the penalty function method in optimization. To improve the uniformity of sample points in the design of experiments (DoE), a distance constraint term is added to the learning function. Besides, an adaptive sampling region scheme is employed to filter out samples with weak probability density to further enhance the efficiency of the proposed algorithm. Moreover, a hybrid stopping criterion based on the error-based stopping criterion using the bootstrap confidence estimation is developed to terminate the active learning process to ensure that the learning algorithm stops at an appropriate stage. The proposed ABSVR is easy to implement since no embedded optimization algorithm nor iso-probabilistic transformation is required. The performance of ABSVR is evaluated using six numerical examples featuring different complexity, and the results demonstrate the superior performance of ABSVR for structural reliability analysis in terms of accuracy and efficiency.

    The correlation between statistical descriptors of heterogeneous materials
    Computer Methods in Applied Mechanics and Engineering, 384 (2021) 113948.
    S.Q. Cui, J.L. Fu, S. Cen, H.R. Thomas and C.F. Li

    Heterogeneous materials such as rocks and composites are comprised of multiple material phases of different sizes and shapes that are randomly distributed through the medium. The random microstructure is typically described by using various statistical descriptors, which include volume fraction, two-point correlation function, and tortuosity, to name a few. Capturing different morphological features, a large number of statistical descriptors are proposed in different research fields, such as material science, geoscience and computational engineering. It is well known that these statistical descriptors are not independent from each other, but until recently it remains unclear what descriptors are more similar or more different. In particular, it is extremely difficult to look for quantified relations between various descriptors, since they are often defined in very different formats. The lack of quantified understanding of descriptors' relations can cause uncertainties or even systematic errors in heterogeneous materials studies. To address this issue, we propose a novel and generic correlation analysis strategy and establish, for the first time, the quantified relations between various statistical descriptors for heterogeneous materials. Based on data science techniques, our approach consists of three operational steps: data regularization, dimension reduction and correlation analysis. A total of 41 statistical descriptors are collected and analysed in this study, which is readily extensible to include other new descriptors. The generic and quantified correlation results are compared with other established descriptor relations that are obtained from analytical analysis or physical intuition, and good agreements are observed in all cases. The quantified relations between various descriptors are summarized in a single correlation graph, which provides useful guiding information for the characterization, reconstruction and property prediction of heterogeneous materials.

    Unified particle model for multiple-fluid flow and porous material
    ACM Transactions on Graphics, 40 (4) (2021) 118.
    B. Xu, C.F. Li and B. Ren

    Porous materials are common in daily life. They include granular material (e.g. sand) that behaves like liquid flow when mixed with fluid and foam material (e.g. sponge) that deforms like solid when interacting with liquid. The underlying physics is further complicated when multiple fluids interact with porous materials involving coupling between rigid and fluid bodies, which may follow different physics models such as the Darcy's law and the multiple-fluid Navier-Stokes equations. We propose a unified particle framework for the simulation of multiple-fluid flows and porous materials. A novel virtual phase concept is introduced to avoid explicit particle state tracking and runtime particle deletion/insertion. Our unified model is flexible and stable to cope with multiple fluid interacting with porous materials, and it can ensure consistent mass and momentum transport over the whole simulation space.

    Hermite polynomial normal transformation for structural reliability analysis
    Engineering Computations, 38 (8) (2021) 3193-3218.
    J.S. Wang, M. Aldosary, S. Cen and C.F. Li

    Purpose – Normal transformation is often required in structural reliability analysis to convert the nonnormal random variables into independent standard normal variables. The existing normal transformation techniques, for example, Rosenblatt transformation and Nataf transformation, usually require the joint probability density function (PDF) and/or marginal PDFs of non-normal random variables. In practical problems, however, the joint PDF and marginal PDFs are often unknown due to the lack of data while the statistical information is much easier to be expressed in terms of statistical moments and correlation coefficients. This study aims to address this issue, by presenting an alternative normal transformation method that does not require PDFs of the input random variables. Design/methodology/approach – The new approach, namely, the Hermite polynomial normal transformation, expresses the normal transformation function in terms of Hermite polynomials and it works with both uncorrelated and correlated random variables. Its application in structural reliability analysis using different methods is thoroughly investigated via a number of carefully designed comparison studies. Findings – Comprehensive comparisons are conducted to examine the performance of the proposed Hermite polynomial normal transformation scheme. The results show that the presented approach has comparable accuracy to previous methods and can be obtained in closed-form. Moreover, the new scheme only requires the first four statistical moments and/or the correlation coefficients between random variables, which greatly widen the applicability of normal transformations in practical problems. Originality/value – This study interprets the classical polynomial normal transformation method in terms of Hermite polynomials, namely, Hermite polynomial normal transformation, to convert uncorrelated/ correlated random variables into standard normal random variables. The new scheme only requires the first four statistical moments to operate, making it particularly suitable for problems that are constraint by limited data. Besides, the extension to correlated cases can easily be achieved with the introducing of the Hermite polynomials. Compared to existing methods, the new scheme is cheap to compute and delivers comparable accuracy.

    Shape-free arbitrary polygonal hybrid stress/displacement function flat shell element for linear and geometrically nonlinear analyses
    International Journal for Numerical Methods in Engineering, 122 (2021): 4172-4218.
    C.J. Wu, S. Cen, R.X. Ma and C.F. Li

    A high-performance shape-free arbitrary polygonal hybrid stress/displacement-function flat shell finite element method is proposed for linear and geometrically nonlinear analyses of shells. First, an arbitrary polygonalMindlin–Reissner plate element and an arbitrary polygonal membrane element with drilling degrees of freedom are constructed based on hybrid displacement-function and hybrid stress-function methods, respectively. Both elements have only two corner nodes along each edge. Second, by assembling the plate and the membrane elements, an arbitrary polygonal flat shell element is constructed. Third, based on the corotational method, a proper best-fit corotated frame for geometrically nonlinear polygonal elements is designed. By updating the analytical trial functions of shell element in each increment step, the original linear flat shell element is generalized to a geometrically nonlinear model. Numerical examples show that the new element possesses excellent performance for both linear and geometrically nonlinear analyses, and possesses outstanding flexibility in dealing with complex loading distributions and mesh shapes.

    A review on proppant transport modelling
    Journal of Petroleum Science and Engineering, 204 (2021) 108753.
    B.R. Barboza, B. Chen and C.F. Li

    Proppant transport is a critical physical process in hydraulic fracturing which has been extensively used for reservoir stimulation in petroleum engineering. Proppants injected together with fracturing fluid provide structural support to the stimulated fracture network and prevent them from closing after flowback. Hence, the final proppant distribution in fracture networks affects directly the effectiveness of hydraulic fracturing. Owing to the limitation and high cost of well logging, computational modelling has been increasingly used to study proppant transport, where different assumptions and numerical models have been employed often without rigorous validation or justification. This work presents a comprehensive review on proppant transport modelling, from relevant physics to numerical approaches, aiming to provide an unbiased global picture of the-state-of-theart studies, inspire new insights, and promote the development of innovative and reliable computational models for proppant transport.

    Generalized conforming Trefftz element for size-dependent analysis of thin microplates based on the modified couple stress theory
    Engineering Analysis with Boundary Elements, 125 (2021) 46-58.
    Y. Shang, Y.H. Mao, S. Cen and C.F. Li

    In this work, a new displacement-based Trefftz plate element is developed for size-dependent bending analysis of the thin plate structures in the context of the modified couple stress theory. This is achieved via two steps. First, the Trefftz functions, that are derived by introducing the thin plate bending assumptions into the three- dimensional governing equations of the modified couple stress elasticity, are adopted as the basis functions for designing the element's displacement interpolations. Second, the generalized conforming theory is employed to meet the interelement compatibility requirements in weak sense for ensuring the convergence property. The resulting 4-node displacement-based plate element performs like nonconforming models on coarse meshes and gradually converges into a conforming one with the mesh refinement. Numerical tests reveal that the new element can efficiently capture the size-dependent mechanical responses of thin microplates and exhibits satisfactory numerical accuracy and distortion tolerance. Moreover, as the element has only three degrees of freedom (DOF) per node, it can be easily incorporated into the commonly available finite element programs.

    Tortuosity of porous media: image analysis and physical simulation
    Earth Science Reviews, 212 (2021) 103439.
    J.L. Fu, H.R. Thomas and C.F. Li

    Tortuosity is widely used as a critical parameter to predict transport properties of porous media, such as rocks and soils. But unlike other standard microstructural properties, the concept of tortuosity is vague with multiple definitions and various evaluation methods introduced in different contexts. Hydraulic, electrical, diffusional, and thermal tortuosities are defined to describe different transport processes in porous media, while geometrical tortuosity is introduced to characterize the morphological property of porous microstructures. In particular, the rapid development of microscopy imaging techniques has made digital microstructures of porous media increasingly accessible, from which geometrical and physical tortuosities can be evaluated using various image analysis and numerical simulation methods. These tortuosities are defined differently and can differ greatly in value, but in many works of literature, they are used interchangeably. To address this situation, we systematically examine geometrical, hydraulic, electrical, diffusional, and thermal tortuosities from the viewpoints of the definition and evaluation method. For the same porous medium, visible discrepancies are found in the evaluated geometrical and physical tortuosities, depending on the specific definition and the evaluation method adopted. This observation makes it questionable to directly use the geometrical tortuosity as a substitute for physical tortuosities, a common practice in the literature. Thus, the correlations between geometrical and physical tortuosities are further analyzed, which also takes into account the influence of both image size and resolution. From the correlation analysis, phenomenological relations between geometrical and physical tortuosities are established, so that the latter can be accurately predicted by using the former which is much cheaper to evaluate from digital microstructures.

    Statistical characterization and reconstruction of heterogeneous microstructures using deep neural network
    Computer Methods in Applied Mechanics and Engineering, 373 (2021) 113516.
    J.L. Fu, S.Q. Cui, S. Cen and C.F. Li

    Heterogeneous materials, whether natural or artificial, are usually composed of distinct constituents present in complex microstructures with discontinuous, irregular and hierarchical characteristics. For many heterogeneous materials, such as porous media and composites, the microstructural features are of fundamental importance for their macroscopic properties. This paper presents a novel method to statistically characterize and reconstruct random microstructures through a deep neural network (DNN) model, which can be used to study the microstructure–property relationships. In our method, the digital microstructure image is assumed to be a stationary Markov random field (MRF), and local patterns covering the basic morphological features are collected to train a DNN model, after which statistically equivalent samples can be generated through a DNN-guided reconstruction procedure. Furthermore, to overcome the short-distance limitation associated with the MRF assumption, a multi-level approach is developed to preserve the long-distance morphological features of heterogeneous microstructures. A large number of tests have been conducted to compare the reconstructed and target microstructures in both morphological characteristics and physical properties, and good agreements are observed in all test cases. The proposed method is efficient, accurate, versatile, and especially beneficial to the statistical reconstruction of 2D/3D microstructures with long-distance correlations.


    A moving least square reproducing kernel particle method for unified multiphase continuum simulation
    ACM Transactions on Graphics, 39 (6) (2020) 176.
    X.S. Chen, C.F. Li, G.C. Cao, Y.T. Jiang and S.M. Hu

    In physically based-based animation, pure particle methods are popular due to their simple data structure, easy implementation, and convenient parallelization. As a pure particle-based method and using Galerkin discretization, the Moving Least Square Reproducing Kernel Method (MLSRK) was developed in engineering computation as a general numerical tool for solving PDEs. The basic idea of Moving Least Square (MLS) has also been used in computer graphics to estimate deformation gradient for deformable solids. Based on these previous studies, we propose a multiphase MLSRK framework that animates complex and coupled fluids and solids in a unified manner. Specifically, we use the Cauchy momentum equation and phase field model to uniformly capture the momentum balance and phase evolution/interaction in a multiphase system, and systematically formulate the MLSRK discretization to support general multiphase constitutive models. A series of animation examples are presented to demonstrate the performance of our new multiphase MLSRK framework, including hyperelastic, elastoplastic, viscous, fracturing and multiphase coupling behaviours etc.

    Resolution effect: an error correction model for intrinsic permeability of porous media estimated from lattice Boltzmann method
    Transport in Porous Media, 132 (3) (2020) 627-656.
    J.L. Fu, J.B. Dong, Y.L. Wang, Y. Ju, D.R.J. Owen and C.F. Li

    In digital rock physics, the intrinsic permeability of a porous rock sample can be evaluated from its micro-computed tomography ( µ-CT) image through lattice Boltzmann method (LBM) simulation. The LBM permeability evaluation has been increasingly adopted by the oil and gas industries, especially when the access to core samples is limited. In order to accurately evaluate the permeability of porous media, this digital approach requires highquality µ-CT images with sufficient resolution and size. In practice, however, the LBM simulation is often performed using images of reduced resolution, due to limitations in computing power and simulation time. As a result, the permeability results obtained are often compromised with significant errors, known as the resolution effect. In this study, the resolution effect is quantitatively investigated to identify the primary causes of error, based on which an error correction model for the LBM permeability evaluation is proposed. The model uses such geometric attributes as connected porosity, specific surface area and diffusion tortuosity to quantify the resolution effect and achieve error correction. Demonstrated on various types of porous media including sandstone, carbonate rock, sand pack, synthesis silica, etc., the proposed error correction model can effectively correct the errors in LBM permeability evaluation due to the resolution effect. Our error correction model makes image resolution reduction more meaningful and creditable for LBM permeability evaluation of porous media, thereby supporting its adoption in practical applications.

    Hyperelastic finite deformation analysis with the unsymmetric finite element method containing homogeneous solutions of linear elasticity
    International Journal for Numerical Methods in Engineering, 121 (16) (2020) 3702-3721.
    Z. Li, S. Cen, J.B. Huang and C.F. Li

    A recent unsymmetric 4-node, 8-DOF plane finite element US-ATFQ4 is generalized to hyperelastic finite deformation analysis. Since the trial functions of US-ATFQ4 contain the homogenous closed analytical solutions of governing equations for linear elasticity, the key of the proposed strategy is how to deal with these linear analytical trial functions (ATFs) during the hyperelastic finite deformation analysis. Assuming that the ATFs can properly work in each increment, an algorithm for updating the deformation gradient interpolated by ATFs is designed. Furthermore, the update of the corresponding ATFs referred to current configuration is discussed with regard to the hyperelastic material model, and a specified model, neo-Hookean model, is employed to verify the present formulation of US-ATFQ4 for hyperelastic finite deformation analysis. Various examples show that the present formulation not only remain the high accuracy andmesh distortion tolerance in the geometrically nonlinear problems, but also possess excellent performance in the compressible or quasi-incompressible hyperelastic finite deformation problems where the strain is large.

    8-node hexahedral unsymmetric element with rotation DOFs for modified couple stress elasticity
    International Journal for Numerical Methods in Engineering, 121 (12) (2020) 2683-2700.
    Y. Shang, C.F. Li and K.Y. Jia

    A new C0 8-node 48-DOF hexahedral element is developed for analysis of size-dependent problems in the context of the modified couple stress theory by extending the methodology proposed in our recent work (Shang et al., Int J NumerMethods Eng 119(9): 807-825, 2019) to the three-dimensional (3D) cases. There are two major innovations in the present formulation. First, the independent nodal rotation degrees of freedom (DOFs) are employed to enhance the standard 3D isoparametric interpolation for obtaining the displacement and strain test functions, as well as to approximatively design the physical rotation field for deriving the curvature test function. Second, the equilibrium stress functions instead of the analytical functions are used to formulate the stress trial function whilst the couple stress trial function is directly obtained from the curvature test function by using the constitutive relationship. Besides, the penalty function is introduced into the virtual work principle for enforcing the C1 continuity condition in weak sense. Several benchmark examples are examined and the numerical results demonstrate that the element can simulate the size-dependent mechanical behaviors well, exhibiting satisfactory accuracy and low susceptibility to mesh distortion.

    A modified STRUCT model for efficient engineering computations of turbulent flows in hydro-energy machinery
    International Journal of Heat and Fluid Flow, 85 (2020) 108628.
    C.Y. Wang, F.J. Wang, C.F. Li, C.L. Ye, T.T. Yan and Z.C. Zou

    A modified STRUCT (MST) turbulence model for efficient engineering computations of turbulent flows in hydro-energy machinery is proposed in this paper. The MST model switches between URANS and LES-like modes using a new damping function to adjust the turbulent viscosity. Compared with the original STRUCT method, the modifications are as follows: (1) the BSL k-? model with the Spalart-Shur correction is chosen as the new baseline to improve the sensitivity to rotation and curvature; (2) a new adaptive time-scale ratio is proposed to avoid the arbitrariness of geometric averaging operation in the original method; (3) the normalized helicity is introduced into the new damping function to detect the energy backscatter phenomenon. Five classical high Reynolds number flow cases are tested. The results show that the turbulent viscosity of the MST model is reasonably reduced in the massively separated regions and LES-like mode is activated, which captures more turbulent vortices and fluctuations on the URANS grids. With high efficiency and robustness, the MST model inherits the advantages of the original STRUCT method and improves the prediction accuracy of the turbulence with rotation and curvature, which enables efficient engineering computations of turbulent flows in hydro-energy machinery.

    Phase-field simulation of hydraulic fracturing with a revised fluid model and hybrid solver
    Engineering Fracture Mechanics, 229 (2020) 106928.
    B. Chen, Y.N. Sun, B.R. Barboza, A.R. Barron and C.F. Li

    With an intrinsic advantage in describing complex fracture networks, the phase field method has demonstrated promising potential for the simulation of hydraulic fracturing processes in recent literatures. We critically examine the existing phase-field hydraulic fracturing models, and propose a hybrid solution scheme with a revised fluid model. Specifically, the formation deformation and phase field are solved using the finite element method (FEM), while the fluid flows are solved using the finite volume method (FVM). The proposed hybrid scheme is validated with the analytical solution for the toughness-dominated fracture propagation and is tested on the complex hydraulic fracturing process in a naturally fractured formation. Demonstrated by numerical examples, the proposed hybrid phase-field framework has several advantages: (1) it captures the effect of fluid pressure inside the fracture and reservoir more accurately than existing models; (2) it provides a sharper capture of formation fractures; (3) it avoids the nonphysical oscillation of fluid pressure when using a pure FEM solver; and (4) it has a superior performance in mesh and time step convergence.

    A divergence-free mixture model for multiphase fluids
    Computer Graphics Forum, 39 (8) 2020 69-77.
    Y.T. Jiang, C.F. Li and S.M. Hu

    We present a novel divergence free mixture model for multiphase flows and the related fluid-solid coupling. The new mixture model is built upon a volume-weighted mixture velocity so that the divergence free condition is satisfied for miscible and immiscible multiphase fluids. The proposed mixture velocity can be solved efficiently by adapted single phase incompressible solvers, allowing for larger time steps and smaller volume deviations. Besides, the drift velocity formulation is corrected to ensure mass conservation during the simulation. The new approach increases the accuracy of multiphase fluid simulation by several orders. The capability of the new divergence-free mixture model is demonstrated by simulating different multiphase flow phenomena including mixing and unmixing of multiple fluids, fluid-solid coupling involving deformable solids and granular materials.


    Development of a multi-continuum quadruple porosity model to estimate CO2 storage capacity and CO2 enhanced shale gas recovery
    Journal of Petroleum Science and Engineering, 178 (2019) 964-974.
    M.L. Wu, M.C. Ding, J. Yao, C.F. Li, X. Li and J.M. Zhu

    Geologic storage of CO2 in shale formation not only enhances natural gas recovery, but also sequestrates CO2 effectively. According to this technology, a multi-continuum quadruple porosity binary component gas model is developed to investigate carbon dioxide storage capacity and CO2 enhanced shale gas recovery, which is based on multiple flow mechanisms, including dissolution, adsorption/desorption, viscous flow, diffusion, slip flow and stress sensitivity of hydraulic fractures. This fully coupled model is divided into quadruple media, including organic matters, organic pore system, matrix system and natural fracture system. The matrix-fracture transfer flow is simulated by modified multiple interacting continua (MINC) method. Embedded discreate fracture model (EDFM) is introduced to describe the gas flow in hydraulic fractures and the transfer flow between hydraulic fractures and natural fractures. Finite difference method (FDM) and quasi-Newton iterative method are applied to solve this model. The reliability and practicability of this model is validated by matching the production history of a fractured horizontal well in shale gas reservoir. The effects of relevant parameters on production curves are analyzed, including adsorption parameters, dissolution parameters, well production pressure, injection pressure, volumetric fraction of kerogen and injection opportunity. The result shows that the model in this work is reliable and practicable, and the model presented here can be used to investigate the injectivity of CO2 and CO2 enhanced shale gas recovery.

    An unsymmetric 8-node hexahedral solid-shell element with high distortion tolerance: Geometric nonlinear formulations
    International Journal for Numerical Methods in Engineering, 120 (5) (2019) 580-606.
    Z. Li, J.B. Huang, S. Cen and C.F. Li

    Arecent distortion-tolerant unsymmetric 8-node hexahedral solid-shell element US-ATFHS8, which takes the analytical solutions of linear elasticity as the trial functions, is successfully extended to geometric nonlinear analysis. This extension is based on the corotational (CR) approach due to its simplicity and high efficiency, especially for geometric nonlinear analysis where the strain is still small. Based on the assumption that the analytical trial functions can properly work in each increment during the nonlinear analysis, the incremental corotational formulations of the nonlinear solid-shell element US-ATFHS8 are derived within the updated Lagrangian (UL) framework, in which an appropriate updated strategy for linear analytical trial functions is proposed. Numerical examples show that the present nonlinear element US-ATFHS8 possesses excellent performance for various rigorous tests no matter whether regular or distorted mesh is used. Especially, it even performs well in some situations that other conventional elements cannot work.

    A simple unsymmetric 4-node 12-DOF membrane element for the modified couple stress theory
    International Journal for Numerical Methods in Engineering, 119 (2019) 807-825.
    Y. Shang, Z.H. Qian, S. Cen and C.F. Li

    In this work, the recently proposed unsymmetric 4-node 12-DOF (degree-of-freedom) membrane element (Shang and Ouyang, Int J Numer Methods Eng 113(10): 1589-1606, 2018), which has demonstrated excellent performance for the classical elastic problems, is further extended for the modified couple stress theory, to account for the size effect of materials. This is achieved via two formulation developments. Firstly, by using the penalty function method, the kinematic relations between the element's nodal drilling DOFs and the true physical rotations are enforced. Consequently, the continuity requirement for themodified couple stress theory is satisfied in weak sense, and the symmetric curvature test function can be easily derived from the gradients of the drilling DOFs. Secondly, the couple stress field that satisfies a priori the related equilibrium equations is adopted as the energy conjugate trial function to formulate the element for themodified couple stress theory. As demonstrated by a series of benchmark tests, the new element can efficiently capture the size-dependent responses of materials and is robust to mesh distortions. Moreover, as the new element uses only three conventional DOFs per node, it can be readily incorporated into the standard finite element program framework and commonly available finite element programs.

    Some advances in high-performance finite element methods
    Engineering Computations, 36 (8) (2019) 2811-2834.
    S. Cen, C.J. Wu, Z. Li, Y. Shang and C.F. Li

    Purpose – The purpose of this paper is to give a review on the newest developments of high-performance finite element methods (FEMs), and exhibit the recent contributions achieved by the authors' group, especially showing some breakthroughs against inherent difficulties existing in the traditional FEMfor a long time. Design/methodology/approach – Three kinds of new FEMs are emphasized and introduced, including the hybrid stress-function element method, the hybrid displacement-function element method for Mindlin– Reissner plate and the improved unsymmetric FEM. The distinguished feature of these three methods is that they all apply the fundamental analytical solutions of elasticity expressed in different coordinates as their trial functions. Findings – The new FEMs show advantages from both analytical and numerical approaches. All the models exhibit outstanding capacity for resisting various severe mesh distortions, and even perform well when other models cannot work. Some difficulties in the history of FEM are also broken through, such as the limitations defined by MacNeal's theorem and the edge-effect problems ofMindlin–Reissner plate. Originality/value – These contributions possess high value for solving the difficulties in engineering computations, and promote the progress of FEM.

    A novel displacement-based Trefftz plate element with high distortion tolerance for orthotropic thick plates
    Engineering Analysis with Boundary Elements, 106 (2019) 452-461.
    Y. Shang, C.F. Li and M.J. Zhou

    In this work, a simple but robust 4-node displacement-based Trefftz plate element is proposed for efficiently analysis of orthotropic plate structures. This is achieved via two steps. First, the element's deflection and rota- tion fields are directly formulated at the base of the Trefftz functions of the orthotropic Mindlin–Reissner plate. Second, to ensure the convergence property and further enhance the element performance, the generalized con- forming theory is employed to satisfy the requirement for interelement compatibility in weak sense. The resulting displacement-based element belongs to the nonconforming models on coarse meshes and gradually converges into a conforming one with the mesh refinement. Numerical benchmarks reveal that the new element can produce high precision results for both displacement and stress resultant. Moreover, it is free of shear locking and has good tolerances to mesh distortions.

    A rigging-skinning scheme to control fluid simulation
    Computer Graphics Forum, 38 (7) (2019) 501-512.
    J.M. Lu, X.S. Chen, X. Yan, C.F. Li, M.C. Lin and S.M. Hu

    Inspired by skeletal animation, a novel rigging-skinning flow control scheme is proposed to animate fluids intuitively and efficiently. The new animation pipeline creates fluid animation via two steps: fluid rigging and fluid skinning. The fluid rig is defined by a point cloud with rigid-body movement and incompressible deformation, whose time series can be intuitively specified by a rigid body motion and a constrained free-form deformation, respectively. The fluid skin generates plausible fluid flows by virtually fluidizing the point-cloud fluid rig with adjustable zero- and first-order flow features and at fixed computational cost. Fluid rigging allows the animator to conveniently specify the desired low-frequency flow motion through intuitive manipulations of a point cloud, while fluid skinning truthfully and efficiently converts the motion specified on the fluid rig into plausible flows of the animation fluid, with adjustable fine-scale effects. Besides being intuitive, the rigging-skinning scheme for fluid animation is robust and highly efficient, avoiding completely iterative trials or time-consuming nonlinear optimization. It is also versatile, supporting both particle- and grid- based fluid solvers. A series of examples including liquid, gas and mixed scenes are presented to demonstrate the performance of the new animation pipeline.

    Predicting the effective mechanical property of heterogeneous materials by image based modeling and deep learning
    Computer Methods in Applied Mechanics and Engineering, 347 (2019) 735-753.
    X. Li, Z.L. Liu, S.Q. Cui, C.C. Luo, C.F. Li and Z. Zhuang

    In contrast to the composition uniformity of homogeneous materials, heterogeneous materials are normally composed of two or more distinctive constituents. It is usually recognized that the effective material property of a heterogeneous material is related to the mechanical property and the distribution pattern of each forming constituent. However, to establish an explicit relationship between the macroscale mechanical property and the microstructure appears to be complicated. On the other hand, machine learning methods are broadly employed to excavate inherent rules and correlations based on a significant amount of data samples. Specifically, deep neural networks are established to deal with situations where input–output mappings are extensively complex. In this paper, a method is proposed to establish the implicit mapping between the effective mechanical property and the mesoscale structure of heterogeneous materials. Shale is employed in this paper as an example to illustrate the method. At the mesoscale, a shale sample is a complex heterogeneous composite that consists of multiple mineral constituents. The mechanical properties of each mineral constituent vary significantly, and mineral constituents are distributed in an utterly random manner within shale samples. Large quantities of shale samples are generated based on mesoscale scanning electron microscopy images using a stochastic reconstruction algorithm. Image processing techniques are employed to transform the shale sample images to finite element models. Finite element analysis is utilized to evaluate the effective mechanical properties of the shale samples. A convolutional neural network is trained based on the images of stochastic shale samples and their effective moduli. The trained network is validated to be able to predict the effective moduli of real shale samples accurately and efficiently. Not limited to shale, the proposed method can be further extended to predict effective mechanical properties of other heterogeneous materials. ?c 2019 ElsevierB.V.All rights reserved.

    Production-performance analysis of composite shale-gas reservoirs by the boundary-element method
    SPE Reservoir Evaluation & Engineering, 22 (1) (2019) 238-252.
    M.L. Wu, M.C. Ding, J. Yao, C.F. Li, Z.Q. Huang and S.N. Xu

    A shael-gas reservoir with a multiple-fractured horizontal well (MFHW) is divided into two regions: The inner region is defined as stimulated reservoir volume (SRV), which is interconnected by the fracture network after fracturing, while the outer region is called unstimulated reservoir volume (USRV), which has not been stimulated by fracturing. Considering an arbitrary interface boundary between SRV and USRV, a composite model is presented for MFHWs in shale-gas reservoirs, which is based on multiple flow mechanisms, inlcuding adsorption/desorption, viscous flow, diffusive flow, and stress sensitivity of natural fractures. The boundary-element method (BEM) is applied to solve the production of MFHWs in shale-gas reservoirs. The accuracy of this model is validated by comparing its production solution with the result derived from an analyticl method and the reservoir simulator. Furthermore, the practicability of this model is validated by matching the production history of the MFHW in a shale-gas revervoir. The result shows that the model in this work is reliable and practicable. The effects of relevant parameters on production curves are analyzed, including Langmuir volume, Langmuir presure, hydraulic-fracutre width, hydraulic-fracture premeability, natural-fracture permeability, matrix permeability, diffusion coefficient, stress-sensitivity coefficient, and the shape of the SRV. The model presented here can be used for production abalysis for shale-gas-reservoir development.


    MPM simulation of interacting fluids and solids
    Computer Graphics Forum, 37 (8) (2018) 183-193.
    X. Yan, C.F. Li, X.S. Chen and S.M. Hu

    The material point method (MPM) has attracted increasing attention from the graphics community, as it combines the strengths of both particle- and grid-based solvers. Like the smoothed particle hydrodynamics (SPH) scheme, MPM uses particles to discretize the simulation domain and represent the fundamental unknowns. This makes it insensitive to geometric and topological changes, and readily parallelizable on a GPU. Like grid-based solvers, MPM uses a background mesh for calculating spatial derivatives, providing more accurate and more stable results than a purely particle-based scheme. MPM has been very successful in simulating both fluid flow and solid deformation, but less so in dealing with multiple fluids and solids, where the dynamic fluid-solid interaction poses a major challenge. To address this shortcoming of MPM, we propose a new set of mathematical and computational schemes which enable efficient and robust fluid-solid interaction within the MPM framework. These versatile schemes support simulation of both multiphase flow and fully-coupled solid-fluid systems. A series of examples is presented to demonstrate their capabilities and performance in the presence of various interacting fluids and solids, including multiphase flow, fluid-solid interaction, and dissolution.

    High-performance unsymmetric 3-node triangular membrane element with drilling DOFs can correctly undertake in-plane moments
    Engineering Computations, 35 (7) (2018) 2543-2556.
    Y. Shang, S. Cen, Z.H. Qian and C.F. Li

    Purpose – This paper aims to propose a simple but robust three-node triangular membrane element with rational drilling DOFs for efficiently analyzing plane problems. Design/methodology/approach – This new element is developed within the general framework of unsymmetric FEM. The element test functions are determined by using a conforming displacement field which is slightly different with the classical Allman's interpolations, while a self-equilibrated stress field formulated based on the analytical airy stress solutions is adopted as the trial functions. To ensure the correctness between the drilling DOFs and the true rotations in elasticity, reasonable constraints are introduced through the penalty function method. Moreover, the special quadrature strategy is used for operating related integrations for future enrichment of element behavior. Findings – Numerical benchmark tests reveal that this new triangular membrane element has exceptional prediction capabilities. In particular, this element can correctly reproduce a rigid body rotation motion and correctly undertake the external in-plane twisting moments; thus, it is a reasonable choice for being used to formulate flat shell elements or to be connected with other kind of elements with physical rotational DOFs. Originality/value – This work provides a new approach for developing high-performance lower-order elements with simple formulations and good numerical accuracies.

    Visual simulation of multiple fluids in computer graphics: a state-of-the-art report
    Journal of Computer Science and Technology, 33 (3) (2018) 431-451.
    B. Ren, X.Y. Yang, M.C. Lin, N. Thuerey, M. Teschner and C.F. Li

    Realistic animation of various interactions between multiple fluids, possibly undergoing phase change, is a challenging task in computer graphics. The visual scope of multi-phase multi-fluid phenomena covers complex tangled surface structures and rich color variations, which can greatly enhance visual effect in graphics applications. Describing such phenomena requires more complex models to handle challenges involving the calculation of interactions, dynamics and spatial distribution of multiple phases, which are often involved and hard to obtain real-time performance. Recently, a diverse set of algorithms have been introduced to implement the complex multi-fluid phenomena based on the governing physical laws and novel discretization methods to accelerate the overall computation while ensuring numerical stability. By sorting through the target phenomena of recent research in the broad subject of multiple fluids, this state-of-the-art report summarizes recent advances on multi-fluid simulation in computer graphics.

    An unsymmetric 8-node hexahedral solid-shell element with high distortion tolerance: linear formulations
    International Journal for Numerical Methods in Engineering, 116 (2018) 759-783.
    J.B. Huang, S. Cen, Z. Li and C.F. Li

    A locking-free unsymmetric 8-node solid-shell element with high distortion tolerance is proposed for general shell analysis, which is equipped with translational degrees of freedom only. The prototype of this newmodel is a recent solid element US-ATFH8 developed by combining the unsymmetric finite element method and the analytical solutions in three-dimensional local oblique coordinates. By introducing proper shell assumptions and assumed natural strain modifications for transverse strains, the new solid-shell element US-ATFHS8 is successfully formulated. This element is able to give highly accurate predictions for shells with different geometric features and loading conditions and is quite insensitive to mesh distortions. In particular, the excellent performance of US-ATFH8 under membrane load is well inherited, which is an outstanding advantage over other shell elements.

    Adaptive finite element analysis for damage detection of non-uniform Euler–Bernoulli beams with multiple cracks based on natural frequencies
    Engineering Computations, 35 (3) (2018) 1203-1229.
    Y.L. Wang, Y. Ju, Z Zhuang and C.F. Li

    Purpose – This study aims to develop an adaptive finite element method for structural eigenproblems of cracked Euler–Bernoulli beams via the superconvergent patch recovery displacement technique. This research comprises the numerical algorithm and experimental results for free vibration problems (forward eigenproblems) and damage detection problems (inverse eigenproblems). Design/methodology/approach – The weakened properties analogy is used to describe cracks in this model. The adaptive strategy proposed in this paper provides accurate, efficient and reliable eigensolutions of frequency and mode (i.e. eigenpairs as eigenvalue and eigenfunction) for Euler–Bernoulli beams with multiple cracks. Based on the frequency measurement method for damage detection, using the difference between the actual and computed frequencies of cracked beams, the inverse eigenproblems are solved iteratively for identifying the residuals of locations and sizes of the cracks by the Newton–Raphson iteration technique. In the crack detection, the estimated residuals are added to obtain reliable results, which is an iteration process that will be expedited by more accurate frequency solutions based on the proposed method for free vibration problems. Findings – Numerical results are presented for free vibration problems and damage detection problems of representative non-uniform and geometrically stepped Euler–Bernoulli beams with multiple cracks to demonstrate the effectiveness, efficiency, accuracy and reliability of the proposed method. Originality/value – The proposed combination of methodologies described in the paper leads to a very powerful approach for free vibration and damage detection of beams with cracks, introducing the mesh refinement, that can be extended to deal with the damage detection of frame structures.

    High-performance geometric nonlinear analysis with the unsymmetric 4-node, 8-DOF plane element US-ATFQ4
    International Journal for Numerical Methods in Engineering, 114 (2018) 931-954.
    Z. Li, S. Cen, C.J. Wu, Y. Shang and C.F. Li

    A recent unsymmetric 4-node, 8-DOF plane element US-ATFQ4, which exhibits excellent precision and distortion-resistance for linear elastic problems, is extended to geometric nonlinear analysis. Since the original linear element US-ATFQ4 contains the analytical solutions for plane pure bending, how to modify such formulae into incremental forms for nonlinear applications and design an appropriate updated algorithm become the key of the whole job. First, the analytical trial functions should be updated at each iterative step in the framework of updated Lagrangian formulation that takes the configuration at the beginning of an incremental step as the reference configuration during that step. Second, an appropriate stress update algorithm in which the Cauchy stresses are updated by the Hughes-Winget method is adopted to estimate current stress fields. Numerical examples show that the new nonlinear element US-ATFQ4 also possesses amazing performance for geometric nonlinear analysis, no matter whether regular or distorted meshes are used. It again demonstrates the advantages of the unsymmetric finite element method with analytical trial functions.

    Numerical investigation of the fluid lag during hydraulic fracturing
    Engineering Computations, 35 (5) (2018) 2050-2077.
    B. Chen, S. Cen, A.R. Barron, D.R.J. Owen and C.F. Li

    Purpose – The purpose of this paper is to systematically investigate the fluid lag phenomena and its influence in the hydraulic fracturing process, including all stages of fluid-lag evolution, the transition between different stages and their coupling with dynamic fracture propagation under common conditions. Design/methodology/approach – A plane 2D model is developed to simulate the complex evolution of fluid lag during the propagation of a hydraulic fracture driven by an impressible Newtonian fluid. Based on the finite element method, a fully implicit solution scheme is proposed to solve the strongly coupled rock deformation, fluid flow and fracture propagation. Using the proposed model, comprehensive parametric studies are performed to examine the evolution of fluid lag in various geological and operational conditions. Findings – The numerical simulations predict that the lag ratio is around 5%or even lower at the beginning stage of hydraulic fracture under practical geological conditions. With the fracture propagation, the lag ratio keeps decreasing and can be ignored in the late stage of hydraulic fracturing for typical parameter combinations. On the numerical aspect, whether the fluid lag can be ignored depends not only on the lag ratio but also on the minimum mesh size used for fluid flow. In addition, an overall mixed-mode fracture propagation factor is proposed to describe the relationship between diverse parameters and fracture curvature. Research limitations/implications – In this study, relatively simple physical models such as linear elasticity for solid, Newtonian model for fluid and linear elasticity fracture mechanics for fracture are used. The current model does not account for such effects like leak off, poroelasticity and softening of rock formations, which may also visibly affect the fluid lag depending on specific reservoir conditions. Originality/value – This study helps to understand the effect of fluid lag during hydraulic fracturing processes and provides numerical experience in dealing with the fluid lag with finite element simulation.

    Propagation of a plane strain hydraulic fracture with a fluid lag in permeable rock
    Journal of Applied Mechanics, 85(9) (2018), 091003.
    B. Chen, A.R. Barron, D.R.J. Owen and C.F. Li

    Based on the KGD scheme, this paper investigates, with both analytical and numerical approaches, the propagation of a hydraulic fracture with a fluid lag in permeable rock. On the analytical aspect, the general form of normalized governing equations is first formulated to take into account both fluid lag and leak-off during the process of hydraulic fracturing. Then a new self-similar solution corresponding to the limiting case of zero dimensionless confining stress (T ¼0) and infinite dimensionless leak-off coefficient (L ¼1) is obtained. A dimensionless parameter R is proposed to indicate the propagation regimes of hydraulic fracture in more general cases, where R is defined as the ratio of the two time-scales related to the dimensionless confining stress T and the dimensionless leak-off coefficient L. In addition, a robust finite element-based KGD model has been developed to simulate the transient process from L ¼ 0 to L ¼1underT ¼0, and the numerical solutions converge and agree well with the self-similar solution atT ¼0 and L ¼1. More general processes fromT ¼0 and L ¼ 0 toT ¼1and L ¼1for three different values of R are also simulated, which proves the effectiveness of the proposed dimensionless parameter R for indicating fracture regimes.

    Structural reliability and stochastic finite element methods: state-of-the-art review and evidence-based comparison
    Engineering Computations, 35 (6) (2018) 2165-2214.
    M. Aldosary, J.S. Wang and C.F. Li

    Purpose – This paper aims to provide a comprehensive review of uncertainty quantification methods supported by evidence-based comparison studies. Uncertainties are widely encountered in engineering practice, arising from such diverse sources as heterogeneity of materials, variability in measurement, lack of data and ambiguity in knowledge. Academia and industries have long been researching for uncertainty quantification (UQ) methods to quantitatively account for the effects of various input uncertainties on the system response. Despite the rich literature of relevant research, UQ is not an easy subject for novice researchers/practitioners, where many different methods and techniques coexist with inconsistent input/ output requirements and analysis schemes. Design/methodology/approach – This confusing status significantly hampers the research progress and practical application of UQ methods in engineering. In the context of engineering analysis, the research efforts of UQ are most focused in two largely separate research fields: structural reliability analysis (SRA) and stochastic finite element method (SFEM). This paper provides a state-of-the-art review of SRA and SFEM, covering both technology and application aspects. Moreover, unlike standard survey papers that focus primarily on description and explanation, a thorough and rigorous comparative study is performed to test all UQ methods reviewed in the paper on a common set of reprehensive examples. Findings – Over 20 uncertainty quantification methods in the fields of structural reliability analysis and stochastic finite element methods are reviewed and rigorously tested on carefully designed numerical examples. They include FORM/SORM, importance sampling, subset simulation, response surface method, surrogate methods, polynomial chaos expansion, perturbation method, stochastic collocation method, etc. The review and comparison tests comment and conclude not only on accuracy and efficiency of each method but also their applicability in different types of uncertainty propagation problems. Originality/value – The research fields of structural reliability analysis and stochastic finite element methods have largely been developed separately, although both tackle uncertainty quantification in engineering problems. For the first time, all major uncertainty quantification methods in both fields are reviewed and rigorously tested on a common set of examples. Critical opinions and concluding remarks are drawn from the rigorous comparative study, providing objective evidence-based information for further research and practical applications.

    New analytical calculation models for the compressive arch action in reinforced concrete structures
    Engineering Structures, 168 (2018) 721-735.
    X.Z. Lu, K.Q. Lin, C.F. Li and Y. Li

    Research challenges associated with progressive collapse of reinforced concrete (RC) structures have attracted growing attention from researchers and industries worldwide, since the 1995 explosion at the Murrah Federal Building in Oklahoma City. The compressive arch action (CAA), as a favorable mechanism to provide the structural resistance to progressive collapse under a column removal scenario, has been extensively studied using both experimental and theoretical approaches. However, the existing prediction models for the CAA resistance are either too complicated or in need of additional information like the peak deformation of the specimen. Another major weakness in the previous CAA calculation models is the negligence of the slab effect, which can contribute significantly to the structural resistance. In this study, based on the finite element analysis of 50 progressive collapse tests reported in the literature and 217 newly designed beam-slab substructures, explicit and easy-to-use CAA calculation models are developed for RC frame beams with and without slabs. The proposed models are validated against both experimental and numerical results with a mean absolute error being less than 10%. The findings from this study can serve to provide a quantitative reference for practical design of RC frame structures against progressive collapse.

    Real-time high-fidelity surface flow simulation
    IEEE Transactions on Visualization and Computer Graphics, 24 (8) (2018) 2411-2423.
    B. Ren, T.L. Yuan, C.F. Li K. Xu and S.M. Hu

    Surface flow phenomena, such as rain water flowing down a tree trunk and progressive water front in a shower room, are common in real life. However, compared with the 3D spatial fluid flow, these surface flow problems have been much less studied in the graphics community. To tackle this research gap, we present an efficient, robust and high-fidelity simulation approach based on the shallow-water equations. Specifically, the standard shallow-water flow model is extended to general triangle meshes with a featurebased bottom friction model, and a series of coherent mathematical formulations are derived to represent the full range of physical effects that are important for real-world surface flow phenomena. In addition, by achieving compatibility with existing 3D fluid simulators and by supporting physically realistic interactions with multiple fluids and solid surfaces, the new model is flexible and readily extensible for coupled phenomena. A wide range of simulation examples are presented to demonstrate the performance of the new approach.


    Improved hybrid displacement function (IHDF) element scheme for analysis of Mindlin-Reissner plate with edge effect
    International Journal for Numerical Methods in Engineering, 111 (12) (2017) 1120-1169.
    Y. Shang, S. Cen, Z. Li and C.F. Li

    For a Mindlin–Reissner plate subjected to transverse loadings, the distributions of the rotations and some resultant forces may vary very sharply within a narrow district near certain boundaries. This edge effect is indeed a great challenge for conventional finite element analysis. Recently, an effective hybrid displacement function (HDF) finite element method was successfully developed for solving such difficulty [1, 2]. Although good performances can be obtained in most cases, the distribution continuity of some resulting resultants is destroyed when coarse meshes are employed. Moreover, an additional local coordinate system must be used for avoiding a singular problem in matrix inversion, which makes the derivations more complicated. Based on amodified complementary energy functional containing Lagrangian multipliers, an improved HDF (IHDF) element scheme is proposed in this work. And two new special IHDF elements, named by IHDF-P4-Free and IHDF-P4-SS1, are constructed for modeling plate behaviors near free and soft simply supported boundaries, respectively. The present modeling scheme not only greatly improves the precision of the numerical results but also avoids usage of the additional local Coordinate system. The numerical tests demonstrate that the new IHDF element scheme is an effective way for solving the challenging edge effect problem in Mindlin–Reissner plates.

    Domain decomposition approach for parallel improvement of tetrahedral meshes
    Journal of Parallel and Distributed Computing, 107 (2017) 101-113.
    J.J. Chen, D.W. Zhao, Y. Xu, Y. Zheng, C.F. Li and J.J. Zheng

    Presently, a tetrahedral mesher based on the Delaunay triangulation approach may outperform a tetrahedral improver based on local smoothing and flip operations by nearly one order in terms of computing time. Parallelization is a feasible way to speed up the improver and enable it to handle large-scale meshes. In this study, a novel domain decomposition approach is proposed for parallel mesh improvement. It analyses the dual graph of the input mesh to build an inter-domain boundary that avoids small dihedral angles and poorly shaped faces. Consequently, the parallel improver can fit this boundary without compromising the mesh quality. Meanwhile, the new method does not involve any interprocessor communications and therefore runs very efficiently. A parallel pre-processing pipeline that combines the proposed improver and existing parallel surface and volume meshers can prepare a quality mesh containing hundreds of millions of elements in minutes. Experiments are presented to show that the developed system is robust and applicable to models of a complication level experienced in industry.

    A new triangular hybrid displacement function element for static and free vibration analyses of Mindlin-Reissner plate
    Latin American Journal of Solids and Structures, 14 (5) (2017) 765-804.
    J.B. Huang, S. Cen, Y. Shang and C.F. Li

    A new 3-node triangular hybrid displacement function Mindlin-Reissner plate element is developed. Firstly, the modified variational functional of complementary energy for Mindlin-Reissner plate, which is eventually expressed by a so-called displacement function F, is proposed. Secondly, the locking-free formulae of Timoshenko's beam theory are chosen as the deflection, rotation, and shear strain along each element boundary. Thirdly, seven fundamental analytical solutions of the displacement function F are selected as the trial functions for the assumed resultant fields, so that the assumed resultant fields satisfy all governing equations in advance. Finally, the element stiffness matrix of the new element, denoted by HDF-P3-7ß, is derived from the modified principle of complementary energy. Together with the diagonal inertia matrix of the 3-node triangular isoparametric element, the proposed element is also successfully generalized to the free vibration problems. Numerical results show that the proposed element exhibits overall remarkable performance in all benchmark problems, especially in the free vibration analyses.

    Distortion-resistant and locking-free eight-node elements effectively capturing the edge effects of Mindlin-Reissner plates
    Engineering Computations, 34 (2) (2017) 548-586.
    Y. Bao, S. Cen and C.F. Li

    Purpose – A simple shape-free high-order hybrid displacement function element method is presented for precise bending analyses of Mindlin–Reissner plates. Three distortion-resistant and locking-free eight-node plate elements are proposed by utilizing thismethod. Design/methodology/approach – This method is based on the principle of minimum complementary energy, in which the trial functions for resultant fields are derived from two displacement functions, F and f, and satisfy all governing equations. Meanwhile, the element boundary displacements are determined by the locking-free arbitrary order Timoshenko's beam functions. Then, three locking-free eight-node, 24-DOF quadrilateral plate-bending elements are formulated: HDF-P8-23b for general cases, HDF-P8-SS1 for edge effects along soft simply supported (SS1) boundary and HDF-P8-FREE for edge effects along free boundary. Findings – The proposed elements can pass all patch tests, exhibit excellent convergence and possess superior precision when compared to all other existing eight-node models, and can still provide good and stable results even when extremely coarse and distorted meshes are used. They can also effectively solve the edge effect by accurately capturing the peak value and the dramatical variations of resultants near the SS1 and free boundaries. The proposed eight-node models possess potential in engineering applications and can be easily integrated into commercial software. Originality/value – This work presents a new scheme, which can take the advantages of both analytical and discrete methods, to develop high-order mesh distortion-resistant Mindlin–Reissner plate-bending elements.

    A constrained optimization solution for Caughey damping coefficients in seismic analysis
    Engineering Computations, 34 (3) (2017) 682-708.
    D.G. Pan and C.F. Li

    Purpose – Extended from the classic Rayleigh damping model in structural dynamics, the Caughey damping model allows the damping ratios to be specified in multiple modes while satisfying the orthogonality conditions. Despite these desirable properties, Caughey damping suffers from a few major drawbacks: depending on the frequency distribution of the significant modes, it can be difficult to choose the reference frequencies that ensure reasonable values for all damping ratios corresponding to the significant modes; it cannot ensure all damping ratios are positive. This paper aims to present a constrained quadratic programming approach to address these issues. Design/methodology/approach – The new method minimizes the error of the structural displacement peak based on the response spectrum theory, while all modal damping ratios are constrained to be greater than zero. Findings – Several comprehensive examples are presented to demonstrate the accuracy and effectiveness of the proposed method, and comparisons with existing approaches are provided whenever possible. Originality/value – The proposed method is highly efficient and allows the damping ratios to be conveniently specified for all significant modes, producing optimal damping coefficients in practical applications.

    Automatic sizing functions for unstructured surface mesh generation
    International Journal for Numerical Methods in Engineering, 109 (4) (2017) 577-608.
    J.J. Chen, Z.F. Xiao, Y. Zheng, J.J. Zheng, C.F. Li and K.W. Liang

    Accurate sizing functions are crucial for efficient generation of high-quality meshes, but to define the sizing function is often the bottleneck in complicated mesh generation tasks because of the tedious user interaction involved. We present a novel algorithm to automatically create high-quality sizing functions for surface mesh generation. First, the tessellation of a Computer Aided Design (CAD) model is taken as the background mesh, in which an initial sizing function is defined by considering geometrical factors and user-specified parameters. Then, a convex nonlinear programming problem is formulated and solved efficiently to obtain a smoothed sizing function that corresponds to a mesh satisfying necessary gradient constraint conditions and containing a significantly reduced element number. Finally, this sizing function is applied in an advancing front mesher. With the aid of a walk-through algorithm, an efficient sizing-value query scheme is developed. Meshing experiments of some very complicated geometry models are presented to demonstrate that the proposed sizing-function approach enables accurate and fully automatic surface mesh generation.

    Finite element analysis for inclined wellbore stability of transversely iso-tropic rock with HMCD coupling based on weak plane strength criterion
    Science China , 59 (2017) 1-14.
    Y.L. Wang, Z. Zhuang, Z.L. Liu, H.L. Yang and C.F. Li

    The finite element analysis (FEA) technology by hydraulic-mechanical-chemical-damage (HMCD) coupling is proposed in this paper for inclined well bore stability analysis of water-sensitive and laminated rock, developed basing on the recently established FEA technology for transversely isotropic rock with hydraulic-mechanical-damage (HMD) coupling. The chemical activity of the drilling fluid is considered as phenomenological hydration behavior, the moisture content and parameters of rock considering hydration could be determined with time. The finite element (FE) solutions of numerical well bore model considering the chemical activity of drilling fluid, damage tensor calculation and weak plane strength criterion for transversely isotropic rock are developed for researching the well bore failure characteristics and computing the time-dependent collapse and fracture pressure of laminated rock as shale reservoirs. A three-dimensional FE model and elastic solid deformation and seepage flow coupled equations are developed, and the damage tensor calculation technology for transversely isotropic rock are realized by introducing effect of the hydration and the stress state under the current load. The proposed method utilizing weak plane strength criterion fully reflects the strength parameters in rock matrix and weak plane. To the end, an effective and reliable numerically three-step FEA strategy is well established for well bore stability analysis. Numerical examples are given to show that the proposed method can establish efficient and applicable FE model and be suitable for analyzing the time-dependent solutions of pore pressure and stresses, and the evolution region considering the hydration surrounding well bore, furthermore to compute the collapse cycling time and the safe mud weight for collapse and fracture pressure of transversely isotropic rock.

    An unsymmetric 8-node hexahedral element with high distortion tolerance
    International Journal for Numerical Methods in Engineering, 109 (8) (2017) 1130-1158.
    P.L. Zhou, J.B. Huang, C.F. Li, Q. Zhang and S. Cen

    Among all 3D 8-node hexahedral solid elements in current finite element library, the 'best' one can produce good results for bending problems using coarse regular meshes. However, once the mesh is distorted, the accuracy will drop dramatically. And how to solve this problem is still a challenge that remains outstanding. This paper develops an 8-node, 24-DOF (three conventional DOFs per node) hexahedral element based on the virtual work principle, in which two different sets of displacement fields are employed simultaneously to formulate an unsymmetrical element stiffness matrix. The first set simply utilizes the formulations of the traditional 8-node trilinear isoparametric element, while the second set mainly employs the analytical trial functions in terms of 3D oblique coordinates (R, S, T). The resulting element, denoted by US-ATFH8, contains no adjustable factor and can be used for both isotropic and anisotropic cases. Numerical examples show it can strictly pass both the first-order (constant stress/strain) patch test and the second-order patch test for pure bending, remove the volume locking, and provide the invariance for coordinate rotation. Especially, it is insensitive to various severe mesh distortions.


    Multiphase SPH simulation for interactive fluids and solids
    ACM Transactions on Graphics (SIGGRAPH 2016), 35 (4) (2016), Article No. 79.
    X. Yan, Y.T. Jiang, C.F. Li, R.R. Martin and S.M. Hu

    This work extends existing multiphase-fluid SPH frameworks to cover solid phases, including deformable bodies and granular materials. In our extended multiphase SPH framework, the distribution and shapes of all phases, both fluids and solids, are uniformly represented by their volume fraction functions. The dynamics of the multiphase system is governed by conservation of mass and momentum within different phases. The behavior of individual phases and the interactions between them are represented by corresponding constitutive laws, which are functions of the volume fraction fields and the velocity fields. Our generalized multiphase SPH framework does not require separate equations for specific phases or tedious interface tracking. As the distribution, shape and motion of each phase is represented and resolved in the same way, the proposed approach is robust, efficient and easy to implement. Various simulation results are presented to demonstrate the capabilities of our new multiphase SPH framework, including deformable bodies, granular materials, interaction between multiple fluids and deformable solids, flow in porous media, and dissolution of deformable solids.

    Statistical reconstruction and Karhunen-Loève expansion for multiphase random media
    International Journal for Numerical Methods in Engineering, 105 (1) (2016) 3-32.
    J.W. Feng, S. Cen, C.F. Li and D.R.J. Owen

    Termed as random media, rocks, composites, alloys and many other heterogeneous materials consist of multiple material phases that are randomly distributed through the medium. This paper presents a robust and efficient algorithm for reconstructing random media, which can then be fed into stochastic finite element solvers for statistical response analysis. The new method is based on nonlinear transformation of Gaussian random fields, and the reconstructed media can meet the discrete-valued marginal probability distribution function and the two-point correlation function of the reference medium. The new method, which avoids iterative root-finding computation, is highly efficient and particularly suitable for reconstructing large-size random media or a large number of samples. Also, benefiting from the high efficiency of the proposed reconstruction scheme, a Karhunen–Loève (KL) representation of the target random medium can be efficiently estimated by projecting the reconstructed samples onto the KL basis. The resulting uncorrelated KL coefficients can be further expressed as functions of independent Gaussian random variables to obtain an approximate Gaussian representation, which is often required in stochastic finite element analysis.

    A 4-node quadrilateral flat shell element formulated by the shape-free HDF plate and HSF membrane elements
    Engineering Computations, 33 (3) (2016) 713-741.
    Y. Shang, S. Cen and C.F. Li

    Purpose – The purpose of this paper is to propose an efficient low-order quadrilateral flat shell element that possesses all outstanding advantages of novel shape-free plate bending and plane membrane elements proposed recently. Design/methodology/approach – By assembling a shape-free quadrilateral hybrid displacement-function (HDF) plate bending element HDF-P4-11ß (Cen et al., 2014) and a shape-free quadrilateral hybrid stress-function (HSF) plane membrane element HSF-Q4?-7ß (Cen et al., 2011b) with drilling degrees of freedom (DOF), a new 4-node, 24-DOF (six DOFs per node) quadrilateral flat shell element is successfully constructed. The trial functions for resultant/stress fields within the element are derived from the analytical solutions of displacement and stress functions for plate bending and plane problems, respectively, so that they can a priori satisfy the related governing equations. Furthermore, in order to take the influences of moderately warping geometry into consideration, the rigid link correction strategy (Taylor, 1987) is also employed. Findings – The element stiffness matrix of a new simple 4-node, 24-DOF quadrilateral flat shell element is obtained. The resulting models, denoted as HDF-SH4, not only possesses all advantages of original HDF plate and HSF plane elements when analyzing plate and plane structures, but also exhibits good performances for analyses of complicated spatial shell structures. Especially, it is quite insensitive to mesh distortions. Originality/value – This work presents a new scheme, which can take the advantages of both analytical and discrete methods, to develop low-order mesh-distortion resistant flat shell elements.

    Fast SPH simulation for gaseous fluids
    The Visual Computer, 32 (4) (2016) 523-534.
    B. Ren, X. Yan, T. Yang, C.F. Li, M.C. Lin and S.M. Hu

    This paper presents a fast smoothed particle hydro-dynamics (SPH) simulation approach for gaseous fluids. Unlike previous SPH gas simulators, which solve the transparent air flow in a fixed simulation domain, the proposed approach directly solves the visible gas without involving the transparent air. By compensating the density and force calculation for the visible gas particles, we completely avoid the need of computational cost on ambient air particles in previous approaches. This allows the computational resources to be exclusively focused on the visible gas, leading to significant performance improvement of SPH gas simulation. The proposed approach is at least ten times faster than the standard SPH gas simulation strategy and is able to reduce the total particle number by 25–400 times in large open scenes. The proposed approach also enables fast SPH simulation of complex scenes involving liquid–gas transition, such as boiling and evaporation. A particle splitting and merging scheme is proposed to handle the degraded resolution in liquid–gas phase transition. Various examples are provided to demonstrate the effectiveness and efficiency of the proposed approach.

    Quasi-static crack propagation modeling using shape-free hybrid stress-function elements with drilling degrees of freedom
    International Journal of Computational Methods, 13 (3) (2016) 1650014.
    S. Cen, Y. Bao and C.F. Li

    A new plane shape-free multi-node singular hybrid stress-function (HS-F) element with drilling degrees of freedom, which can accurately capture the stress intensity factors at the crack tips, is developed. Then, a quasi-static 2D crack propagation modeling strategy is established by combination of the new singular element and a shape-free 4-node HS-F plane element with drilling degrees of freedom proposed recently. Only simple re-meshing with an unstructured mesh is needed for each simulation step. Numerical results show that the proposed scheme is an effective and robust technique for dealing with the crack propagation problems.


    A simple approach for bubble modelling from multiphase fluid simulation
    Computational Visual Media, 1 (2) (2015) 171-181.
    B. Ren, Y.T. Jiang, C.F. Li, T. Yang and Ming C. Lin

    This article presents a novel and flexible bubble modeling technique for multi-fluid simulations using a volume fraction representation. By combining the volume fraction data obtained from a primary multi-fluid simulation with simple and efficient secondary bubble simulation, a range of real-world bubble phenomena are captured with a high degree of physical realism, including large bubble deformation, sub-cell bubble motion, bubble stacking over the liquid surface, bubble volume change, dissolving of bubbles, etc. Without any change in the primary multi-fluid simulator, our bubble modeling approach is applicable to any multi-fluid simulator based on the volume fraction representation.

    An unsymmetric 4-node, 8-DOF plane membrane element perfectly breaking through MacNeal's theorem
    International Journal for Numerical Methods in Engineering, 103 (7) (2015) 469-500.
    S. Cen, P.L. Zhou, C.F. Li and C.J. Wu

    Among numerous finite element techniques, few models can perfectly (without any numerical problems) break through MacNeal’s theorem: any 4-node, 8-DOF membrane element will either lock in in-plane bending or fail to pass a C0 patch test when the element’s shape is an isosceles trapezoid. In this paper, a 4-node plane quadrilateral membrane element is developed following the unsymmetrical formulation concept, which means two different sets of interpolation functions for displacement fields are simultaneously used. The first set employs the shape functions of the traditional 4-node bilinear isoparametric element, while the second set adopts a novel composite coordinate interpolation scheme with analytical trail function method, in which the Cartesian coordinates .x; y/ and the second form of quadrilateral area coordinates (QACM-II) (S,T) are applied together. The resulting element US-ATFQ4 exhibits amazing performance in rigorous numerical tests. It is insensitive to various serious mesh distortions, free of trapezoidal locking, and can satisfy both the classical first-order patch test and the second-order patch test for pure bending. Furthermore, because of usage of the second form of quadrilateral area coordinates (QACM-II), the new element provides the invariance for the coordinate rotation. It seems that the behaviors of the present model are beyond the well-known contradiction defined by MacNeal’s theorem.

    Automatic surface repairing, defeaturing and meshing algorithms based on an extended B-Rep
    Advances in Engineering Software, 86 (2015) 55-69.
    J.J. Chen, B.W. Cao, Y. Zheng, L.J. Xie, C.F. Li and Z.F. Xiao

    This paper presents an extended surface boundary representation (B-rep), where each topology entity can have dual geometric representations to accommodate various defects (e.g., gaps and overlaps) commonly present in CAD models. Keeping a uniform B-rep and the unsuppressed geometry data enables the use of various existing repairing, defeaturing and meshing algorithms to process CAD models with small gaps and overlaps on surface boundaries. The continuous geometry of the input model remains untouched in the repairing, defeaturing and meshing process, and the output mesh is loyal to this geometry. Such feature is often desirable in numerical simulations that require meshes with high geometry

    Two generalized conforming quadrilateral Mindlin–Reissner plate elements based on the displacement function
    Finite Element in Analysis and Design, 99 (2015) 24-38.
    Y. Shang, S. Cen, C.F. Li and X.R. Fu

    This work presents two 4-node, 12-DOF quadrilateral displacement-based finite elements for analysis of the Mindlin-Reissner plate. Derived from the fundamental analytical solutions of the displacement function F, the deflection and rotation fields of the proposed elements satisfy a priori all related governing equations. The unknown coefficients are determined through the generalized conforming element method, a relaxed and rational conforming approach. The resulting elements perform like nonconforming elements on a coarse mesh, and with mesh refinement they converge as conforming elements. Numerical benchmarks demonstrate that the new elements are insensitive to mesh distortion and free of shear locking, and can provide satisfactory results for most cases.

    An effective hybrid displacement function element method for solving the edge effect of Mindlin-Reissner plate
    International Journal for Numerical Methods in Engineering, 102 (8) (2015) 1449-1487.
    Y. Shang, S. Cen, C.F. Li and J.B. Huang

    In bending problems of Mindlin-Reissner plate, the resultant forces often vary dramatically within a narrow range near free and soft simply-supported (SS1) boundaries. This is so-called the edge effect or the boundary layer effect, a challenging problem for conventional finite element method. In this paper, an effective finite element method for analysis of such edge effect is developed. The construction procedure is based on the hybrid displacement function (HDF) element method [1], a simple hybrid-Trefftz stress element method proposed recently. What is different is that an additional displacement function f related to the edge effect is considered, and its analytical solutions are employed as the additional trial functions for the first time. Furthermore, the free and the SS1 boundary conditions are also applied to modify the element assumed resultant fields. Then, two new special elements, HDF-P4-Free and HDF-P4-SS1, are successfully constructed. These new elements are allocated along the corresponding boundaries of the plate, while the other region is modeled by the usual HDF plate element HDF-P4-11 \beta [1]. Numerical tests demonstrate that the present method can effectively capture the edge effects and exactly satisfy the corresponding boundary conditions by only using relatively coarse meshes.

    A priori error estimation for the stochastic perturbation method
    Computer Methods in Applied Mechanics and Engineering, 286 (2015) 1-21.
    X.Y. Wang, S. Cen, C.F. Li and D.R.J. Owen

    The perturbation method has been among the most popular stochastic finite element methods due to its simplicity and efficiency. The error estimation for the perturbation method is well established for deterministic problems, but until now there has not been an error estimation developed in the probabilistic context. This paper presents a priori error estimation for the perturbation method in solving stochastic partial differential equations. The physical problems investigated here come from linear elasticity of heterogeneous materials, where the material parameters are represented by stochastic fields. After applying the finite element discretization to the physical problem, a stochastic linear algebraic equation system is formed with a random matrix on the left hand side. Such systems have been efficiently solved by using the stochastic perturbation approach, without knowing how accurate/inaccurate the perturbation solution is. In this paper, we propose a priori error estimation to directly link the error of the solution vector with the variation of the source stochastic field. A group of examples are presented to demonstrate the effectiveness of the proposed error estimation.


    The mechanical behavior and reliability prediction of the HTR graphite component at various temperature and neutron dose ranges
    Nuclear Engineering and Design 276 (2014) 9-18.
    X. Fang, S.Y. Yu, H.T. Wang and C.F. Li

    In a pebble-bed high temperature gas-cooled reactor (HTR), nuclear graphite serves as the main structuralmaterial of the side reflectors. The reactor core is made up of a large number of graphite bricks. Inthe normal operation case of the reactor, the maximum temperature of the helium coolant commonlyreaches about 750◦C. After around 30 years' full power operation, the peak value of in-core fast neutron-cumulative dose reaches to 1 × 1022n cm−2(EDN). Such high temperature and neutron irradiation strongly impact the behavior of graphite component, causing obvious deformation. The temperature and neutron-dose are unevenly distributed inside a graphite brick, resulting in stress concentrations. The deformation and stress concentration can both greatly affect safety and reliability of the graphite component. In addition, most of the graphite properties (such as Young's modulus and coefficient of thermal expansion)change remarkably under high temperature and neutron irradiations. The irradiation-induced creep also plays a very important role during the whole process, and provides a significant impact on the stress accumulation. In order to simulate the behavior of graphite component under various in-core conditions,all of the above factors must be considered carefully. In this paper, the deformation, stress distribution and failure probability of a side graphite component are studied at various temperature points and neutron-dose levels. 400◦C, 500◦C, 600◦C and 750◦C are selected as the core-side temperature of the graphite component, while the range of neutron dose is 0–1 × 1022n cm−2(EDN). It is shown in the numerical-analysis that temperature strongly affects the histories of both stresses and failure probabilities as fast-neutron dose increases. In addition, the behavior of graphite brick at 750◦C is observed obviously different from that at other temperatures.

    A quasi-static crack propagation simulation based on shape-free hybrid stress-function finite elements with simple remeshing
    Computer Methods in Applied Mechanics and Engineering 275 (2014) 159-188.
    M.J. Zhou, S. Cen, Y. Bao and C.F. Li

    In this paper, a new shape-free multi-node singular hybrid stress-function (HSF) element and a shape-free 8-node plane HSF element proposed recently are employed to simulate the quasi-static 2D crack propagation problem. Compared with other well-known methods, such new scheme exhibits four advantages: (i) for the singular element, the shape and the number of nodes can be flexibly adjusted as required; (ii) high precision for stress intensity factors (SIF) can be obtained due to the advantages of the HSF method; (iii) only simple remeshing with a very coarse mesh is needed for each simulation step; (iv) unstructured mesh containing extremely distorted elements can be used without losing precision. It demonstrates that the proposed scheme is an effective technique for dealing with crack propagation problems.

    Multiple-fluid SPH Simulation Using a Mixture Model
    ACM Transactions on Graphics / SIGGRAPH, 33 (5) (2014) 171.
    R. Bo, C.F. Li, X. Yan, M.C. Lin, J. Bonet and S.M. Hu

    This paper presents a versatile and robust SPH simulation approach for multiple-fluid flows. The spatial distribution of different phases or components is modeled using the volume fraction representation, the dynamics of multiple-fluid flows is captured by using an improved mixture model, and a stable and accurate SPH formulation is rigorously derived to resolve the complex transport and transformation processes encountered in multiple-fluid flows. The new approach can capture a wide range of realworld multiple-fluid phenomena, including mixing/unmixing of miscible and immiscible fluids, diffusion effect and chemical reaction etc. Moreover, the new multiple-fluid SPH scheme can be readily integrated into existing state-of-the-art SPH simulators, and the multiple-fluid simulation is easy to set up. Various examples are presented to demonstrate the effectiveness of our approach.

    Statistical reconstruction of two-phase random media
    Computers and Structures, 137 (2014) 78-92.
    J.W. Feng, C.F. Li, S. Cen and D.R.J. Owen

    A robust and efficient algorithm is proposed to reconstruct two-phase composite materials with random morphology, according to given samples or given statistical characteristics. The new method is based on nonlinear transformation of Gaussian random fields, where the correlation of the underlying Gaussian field is determined explicitly rather than through iterative methods. The reconstructed media can meet the binary-valued marginal probability distribution function and the two point correlation function of the reference media. The new method, whose main computation is completed using fast Fourier transform (FFT), is highly efficient and particularly suitable for reconstructing large size random media or a large number of samples. Its feasibility and performance are examined through a series of practical examples with comparisons to other state-of-the-art methods in random media reconstruction.

    Hybrid displacement function element method: a simple hybrid-Trefftz stress element method for analysis of Mindlin–Reissner plate
    International Journal for Numerical Methods in Engineering, 98 (2014) 203-234.
    S. Cen, Y. Shang, C.F. Li and H.G. Li

    In order to develop robust finite element models for analysis of thin and moderately thick plates, a simple hybrid displacement function element method is presented. First, the variational functional of complementary energy for Mindlin-Reissner plates is modified to be expressed by a displacement function F, which can be used to derive displacement components satisfying all governing equations. Second, the assumed element resultant force fields, which can satisfy all related governing equations, are derived from the fundamental analytical solutions of F . Third, the displacements and shear strains along each element boundary are determined by the locking-free formulae based on the Timoshenko's beam theory. Finally, by applying the principle of minimum complementary energy, the element stiffness matrix related to the conventional nodal displacement DOFs is obtained. Because the trial functions of the domain stress approximations a priori satisfy governing equations, this method is consistent with the hybrid-Trefftz stress element method. As an example, a 4-node, 12-DOF quadrilateral plate bending element, HDF-P4-11?, is formulated. Numerical benchmark examples have proved that the new model possesses excellent precision. It is also a shape-free element that performs very well even when a severely distorted mesh containing concave quadrilateral and degenerated triangular elements is employed.

    Artistic preprocessing for painterly rendering and image stylization
    The Visual Computer, 30 (2014) 969-979.
    Y. Zang, H. Huang and C.F. Li

    A practical image enhancing technique is presented as a preprocessing step for painterly rendering and image stylization, which transforms the input image mimicking the vision of artists. The new method contains mainly two parts dealing with artistic enhancement and color adjustment, respectively. First, through feature extraction and simplification, an abstract shadow map is constructed for the input image, which is then taken as a guide for emphasizing the light–shadow contrast and the important shadow lines. Next, to simulate the intense color emotion often subjectively added by the artist, a color adjustment technique is proposed to generate lively colors with sharp contrast imitating the artistic vision. The preprocessing operation is compatible with existing stylization and stroke-based painterly rendering techniques, and it can also produce different types of stylization results independently.


    Generalized Neumann expansion and its application in stochastic finite element methods
    Mathematical Problems in Engineering, Volume 2013 (2013), Article ID 325025,
    X.Y. Wang, S. Cen and C.F. Li

    An acceleration technique, termed generalized Neumann expansion (GNE), is presented for evaluating the responses of uncertain systems. The GNE method, which solves stochastic linear algebraic equations arising in stochastic finite element analysis, is easy to implement and is of high efficiency. The convergence condition of the new method is studied, and a rigorous error estimator is proposed to evaluate the upper bound of the relative error of a given GNE solution. It is found that the third-order GNE solution is sufficient to achieve a good accuracy even when the variation of the source stochastic field is relatively high. The relationship between the GNEmethod, the perturbationmethod, and the standard Neumann expansion method is also discussed. Based on the links between these three methods, quantitative error estimations for the perturbation method and the standard Neumannmethod are obtained for the first time in the probability context.

    Stroke style analysis for painterly rendering
    Journal of Computer Science and Technology, 28 (5) (2013) 762-775.
    Y. Zang, H. Huang and C.F. Li

    We propose a novel method that automatically analyzes stroke-related artistic styles of paintings. A set of adaptive interfaces are also developed to connect the style analysis with existing painterly rendering systems, so that the specific artistic style of a template painting can be effectively transferred to the input photo with minimal effort. Different from conventional texture-synthesis based rendering techniques that focus mainly on texture features, this work extracts, analyzes and simulates high-level style features expressed by artists' brush stroke techniques. Through experiments, user studies and comparisons with ground truth, we demonstrate that the proposed style-orientated painting framework can significantly reduce tedious parameter adjustment, and it allows amateur users to efficiently create desired artistic styles simply by specifying a template painting.

    Flow field modulation
    IEEE Transaction on Visualization and Computer Graphics, 19 (10) (2013) 1708-1719.
    B. Ren, C.F. Li, Ming C. Lin, T. Kim and S.M. Hu

    The nonlinear and non-stationary nature of Navier-Stokes equations produces fluid flows that can be noticeably different in appearance with subtle changes. In this paper we introduce a method that can analyze the intrinsic multiscale features of flow fields from a decomposition point of view, by using the Hilbert-Huang transform method on 3D fluid simulation. We show how this method can provide insights to flow styles and help modulate the fluid simulation with its internal physical information.We provide easy-to-implement algorithms that can be integrated with standard grid-based fluid simulation methods, and demonstrate how this approach can modulate the flow field and guide the simulation with different flow styles. The modulation is straightforward and relates directly to the flow's visual effect, with moderate computational overhead.

    View-dependent multiscale fluid simulation
    IEEE Transaction on Visualization and Computer Graphics, 19 (2) (2013) 178-188.
    Y. Gao, C.F. Li, B. Ren and S.M. Hu

    Fluid flows are highly nonlinear and non-stationary, with turbulence occurring and developing at different length and time scales. In real-life observations, the multiscale flow generates different visual impacts depending on the distance to the viewer. We propose a new fluid simulation framework that adaptively allocates computational resources according to the viewer's position. First, a 3D empirical mode decomposition scheme is developed to obtain the velocity spectrum of the turbulent flow. Then, depending on the distance to the viewer, the fluid domain is divided into a sequence of nested simulation partitions. Finally, the multiscale fluid motions revealed in the velocity spectrum are distributed non-uniformly to these view-dependent partitions, and the mixed velocity fields defined on different partitions are solved separately using different grid sizes and time steps. The fluid flow is solved at different spatial-temporal resolutions, such that higher-frequency motions closer to the viewer are solved at higher resolutions and vice versa. The new simulator better utilizes the computing power, producing visually plausible results with realistic fine-scale details in a more efficient way. It is particularly suitable for large scenes with the viewer inside the fluid domain. Also, as high-frequency fluid motions are distinguished from low-frequency motions in the simulation, the numerical dissipation is effectively reduced.


    Failure probability study of HTR graphite component using microstructure-based model
    Nuclear Engineering and Design, 253 (2012) 192-199.
    S.Y. Yu, X. Fang, H.T. Wang and C.F. Li

    Nuclear graphite is widely used as in-core structural material in the high temperature gas-cooled reactors (HTRs). As mechanical properties of graphite change with neutron irradiation and temperature, the reliability evaluation of graphite internals of a HTR is commonly accomplished by the finite element analysis using various constitutive creep models and the following structural failure prediction with introduction of failure models that have essential impact on the final evaluation result. In this paper, a microstructurebased failure model proposed by Burchell in 1996 in conjunction with the finite element code INET-GRA3D of INET is used to study failure probability of a graphite side reflector design in the pebble-bed HTR during its entire service life from the microstructural view. The corresponding irradiation-induced creep stress analysis is carried out using both the UKAEA model and the Kennedy model. H-451 graphite is selected as material input since its irradiation material data from both macrostructure tests and microstructure observation is available. The results are compared with those using the macrostructure-based Weibull model commonly accepted in the engineering field and consistent trends are observed. Moreover, The Burchell model leads to the failure probability more sensitive to stress levels than the Weibull model.

    An evaluation on fatigue crack growth in a fine-grained isotropic graphite, Nuclear Engineering and Design
    Nuclear Engineering and Design, 250 (2012) 197-206.
    H.T. Wang, L.B. Sun, C.F. Li, H.T. Wang and L. Shi

    The aim of this paper is to investigate the mechanism of fatigue crack propagation in IG-11 graphite, and determine the crack growth rate in relation to the stress level. Experimental studies were performed at both micro and macro scales. For fatigue microcrack propagation, single-edge-notch specimens were chosen for testing and the fatigue crack growth was measured in situ with a scanning electron microscope. For fatigue macrocrack propagation, CT specimens were used and the fatigue crack growth was measured with a high-accuracy optic microscope. Combining the two groups of experimental results, the following conclusions are derived: (1) The heterogeneous microstructures of the graphite material have significant impact on the fatigue microcrack growth, while their influence on fatigue macrocrack growth is very limited. (2) The relationship between the fatigue crack growth rate and the crack-tip stress intensity factor range can be expressed in the form of Paris formulae, which contains three stages: an initial rising part with a small slope, an abrupt rise with a very large acceleration, and a short final part with a small slope. (3) The fatigue fracture surface of the graphite material contains considerable sliding of leaf-shape graphite flakes combined with small cotton-shape plastic deformations. These sliding traces are approximately aligned at 45?, showing the main cause of the fatigue fracture is the shear stress. There are also a large amount of secondary cracks inside unit cells and on cell walls, which indicates the fracture mechanism of the IG-11 graphite is mainly cleavage fracture.

    A novel joint diagonalization approach for linear stochastic systems and reliability analysis
    Engineering Computations, 29 (2) (2012) 221-244.
    F. Wang, C.F. Li, J.W. Feng, S. Cen and D.R.J. Owen

    We present a novel gradient based iterative algorithm for the joint diagonalization of a set of real symmetric matrices. The approximate joint diagonalization of a set of matrices is an important tool for solving stochastic linear equations. As an application, reliability analysis of structures by using the stochastic finite element analysis based on the joint diagonalization approach is also introduced in this paper, and it provides useful references to practical engineers. By starting with a least squares (LS) criterion, we obtain a classical nonlinear cost-function and transfer the joint diagonalization problem into a least squares like minimization problem. A gradient method for minimizing such a cost function is derived and tested against other techniques in engineering applications. A novel approach is presented for joint diagonalization for a set of real symmetric matrices. The new algorithm works on the numerical gradient base, and solves the problem with iterations. Demonstrated by examples, the new algorithm shows the merits of simplicity, effectiveness, and computational efficiency.

    The various creep models for irradiation behavior of nuclear graphite
    Nuclear Engineering and Design, 242 (2012) 19-25.
    X. Fang, H.T. Wang, S.Y. Yu and C.F. Li

    Graphite is a widely used material in nuclear reactors, especially in high temperature gascooled Q2 reactors (HTRs), in which it plays three main roles: moderator, reflector and structure material. Irradiationinduced creep has a significant impact on the behavior of nuclear graphite as graphite is used in high temperature and neutron irradiation environments. Thus the creep coefficient becomes a key factor in stress analysis and lifetime prediction of nuclear graphite. Numerous creep models have been established, including the visco-elastic model, UK model, and Kennedy model. A Fortran code based on user subroutines of MSC.MARC was developed in INET in order to perform three-dimensional finite element analysis of irradiation behavior of the graphite components for HTRs in 2008, and the creep model used is for the visco-elastic model only. Recently the code has been updated and can be applied to two other models—the UK model and the Kennedy model. In the present study, all three models were used for calculations in the temperature range of 280-450 ◦C and the results are contrasted. The associated constitutive law for the simulation of irradiated graphite covering properties, dimensional changes, and creep is also briefly reviewed in this paper. It is shown that the trends of stresses and life prediction of the three models are the same, but in most cases the Kennedy model gives the most conservative results while the UK model gives the least conservative results. Additionally, the influence of the creep strain ratio is limited, while the absence of primary creep strain leads to a great increase of failure probability.

    Fisheye video correction
    IEEE Transaction on Visualization and Computer Graphics, 18 (10) (2012) 1771-1783.
    J. Wei, C.F. Li, S.M. Hu, R.R. Martin and C.L. Tai

    Various types of video can be captured with fisheye lenses, particularly surveillance video, due to their ability to capture a wide field of view. However, fisheye lenses introduce distortion, which changes as objects in the scene move. This makes fisheye video difficult to interpret. Current still fisheye image correction methods are either limited to small angles of view, or are strongly content-dependent, and therefore unsuitable for processing video streams. We present an efficient and robust scheme for fisheye video correction, which minimizes time-varying distortions and preserves salient content in a coherent manner. Our optimization process is controlled by user annotation, and takes into account a wide set of measures addressing different aspects of natural scene appearance. Each is represented as a quadratic term in an energy minimization problem, leading to a closed form solution via a sparse linear system. We illustrate our method with a range of examples, demonstrating coherent natural-looking video output; the visual quality of individual frames is comparable to those produced by state-of-the-art methods for fisheye still photograph correction.


    Painterly rendering with content-dependent natural paint strokes
    The Visual Computer, 27 (9) (2011) 861-871
    H. Huang, T.N. Fu and C.F. Li

    We present a new painterly rendering method that simulates artists' content-dependent painting process and the natural variation of hand-painted strokes. First, a new stroke layout strategy is proposed to enhance the contrast between large and small paint strokes, which is an important characteristic of hand-painted paintings. Specifically, the input image is partitioned into nonuniform grids according to its importance map, and determined by the grid size, an individually constructed paint stroke is applied in each grid. Second, an anisotropic digital brush is designed to simulate a real paint brush. In particular, each bristle of the digital brush has an individual color, so that strokes rendered by the new brush can have multiple colors and naturally varied textures. Finally, we present a novel method to add lighting effects to the canvas. This lighting imitation method is robust and very easy to implement, and it can significantly improve the quality of rendering. Comparing with traditional painterly rendering approaches, the new method simulates more closely the real painting procedure, and our experimental results show that it can produce vivid paintings with fewer artifacts.

    Shape-free finite element method: The plane hybrid stress-function (HS-F) element method for anisotropic materials
    Science China - Physics, Mechanics & Astronomy, 54 (4) (2011) 653-665.
    S. Cen, X.R. Fu, G.H. Zhou M.J. Zhou and C.F. Li

    The sensitivity problem to mesh distortion and the low accuracy problem of the stress solutions are two inherent difficulties in the finite element method. By applying the fundamental analytical solutions (in global Cartesian coordinates) to the Airy stress function φ of the anisotropic materials, 8- and 12- node plane quadrilateral hybrid stress-function (HS-F) elements are successfully developed based on the principle of the minimum complementary energy. Numerical results show that the present new elements exhibit much better and more robust performance in both displacement and stress solutions than those obtained from other models. They can still perform very well even when the element shapes degenerate into a triangle and a concave quadrangle. It is also demonstrated that the proposed construction procedure is an effective way for developing shape-free finite element models which can completely overcome the sensitivity problem to mesh distortion and can produce highly accurate stress solutions.


    A joint diagonalisation approach for linear stochastic systems
    Computers and Structures, 88 (19-20) (2010) 1137-1148
    C.F. Li, S. Adhikari, S. Cen, Y.T. Feng and D.R.J. Owen

    In stochastic finite element problems the solution of a system of coupled linear random algebraic equations is needed. This problem in turn requires the calculation of the inverse of a random matrix. Over the past four decades several approximate analytical methods and simulation methods have been proposed for the solution of this problem in the context of probabilistic structural mechanics. In this paper we present a new solution method for stochastic linear equations. The proposed method is based on Neumann expansion and the recently developed joint diagonalisation solution strategy. Unlike the classical Neumann expansion, here only the inversion of a diagonal matrix is needed. Numerical examples are given to illustrate the use of the expressions derived in the paper.

    A hybrid-stress element based on Hamilton principle
    Acta Mechanica Sinica, 26 (4) (2010) 625-634  
    S. Cen, T. Zhang, C.F. Li, X.R. Fu and Y.Q. Long

    A novel hybrid-stress finite element method is proposed for constructing simple 4-node quadrilateral plane elements, and the new element is denoted as HH4-3ß here. Firstly, the theoretical basis of the traditional hybrid-stress elements, i.e., the Hellinger-CReissner variational principle, is replaced by the Hamilton variational principle, in which the number of the stress variables is reduced from 3 to 2. Secondly, three stress parameters and corresponding trial functions are introduced into the system equations. Thirdly, the displacement fields of the conventional bilinear isoparametric element are employed in the new models. Finally, from the stationary condition, the stress parameters can be expressed in terms of the displacement parameters, and thus the new element stiffness matrices can be obtained. Since the required number of stress variables in the Hamilton variational principle is less than that in the Hellinger-CReissner variational principle, and no additional incompatible displacement modes are considered, the new hybrid-stress element is simpler than the traditional ones. Furthermore, in order to improve the accuracy of the stress solutions, two enhanced post-processing schemes are also proposed for element HH4-3ß. Numerical examples show that the proposed model exhibits great improvements in both displacement and stress solutions, implying that the proposed technique is an effective way for developing simple finite element models with high performance.

    Example-based painting guided by color features
    Visual Computer, 26 (6-8) (2010) 933-942  
    H. Huang, Y. Zang and C.F. Li

    In this paper, by analyzing and learning the color features of the reference painting with a novel set of measures, an example-based approach is developed to transfer some key color features from the template to the source image. First, color features of a given template painting is analyzed in terms of hue distribution and the overall color tone. These features are then extracted and learned by the algorithm through an optimization scheme. Next, to ensure the spatial coherence of the final result, a segmentation based post processing is performed. Finally, a new color blending model, which avoids the dependence of edge detection and adjustment of inconvenient tune parameters, is developed to provide a flexible control for the accuracy of painting. Experimental results show that the new example-based painting system can produce paintings with specific color features of the template, and it can also be applied to changing color themes of art pieces, designing color styles of paintings/ real images, and specific color harmonization.

    The elastic T-stress for slightly curved or kinked cracks
    International Journal of Solids and Structures, 47 (2010) 1753-1763
    D.F. Li, C.F. Li, H. Qing and J. Lu

    This work presents a solution for the elastic T-stress at the tip of a slightly curved or kinked crack based on a perturbation approach. Compared to other exact or numerical solutions the present solution is accurate for considerable deviations from straightness. The T-stress variation as crack extends along a curved trajectory is subsequently examined. It is predicted that T-stress always keeps negative during crack extension when the crack has an initial negative T-stress. In the case of a positive T-stress and non-zero first and second stress intensity factors initially accompanying the crack, the T-stress is not positive with increasing the extension length until a threshold is exceeded. Based on directional stability criterion with respect to the sign of the T-stress, this result implies that for a straight crack with a positive T-stress, the crack extension path will not turn immediately and instead keep a stable growth until a critical length is reached. This prediction is consistent with experimental observations.

    Axisymmetric vibration of single-walled carbon nanotubes in water
    Physics Letters A, 373 (2010) 2467-2474
    C. Y. Wang, C. F. Li and S. Adhikari

    A double shell-Stokes flow model is developed to study the axisymmetric vibration of single-walled carbon nanotubes (SWCNTs) immerged in water. In contrast to macroscopic solid-liquid system, a submerged SWCNT is coupled with surrounding water via the van der Waals interaction. It is shown that this unique feature substantially reduces viscous damping of the axisymmetric radial, longitudinal and torsional vibrations and significantly up-shifts the frequency of the radial vibration of an SWCNT. The study offers a theoretical explanation for the experimental observation and molecular dynamics simulations available in particular cases, and provides an efficient modeling tool and useful guidance for the study of the general dynamic behavior of SWCNTs in a fluid.

    A directed Monte Carlo solution of linear stochastic algebraic system of equations
    Finite Elements in Analysis and Design, 46 (2010) 462-473
    Y.T. Feng, C.F. Li and D.R.J. Owen

    This paper proposes a modified Monte Carlo simulation method for the solution of a linear stochastic algebraic system of equations arising from the stochastic finite element modelling of linear elastic problems. The basic idea is to direct Monte Carlo samples along straight lines and then utilize their spatial proximity or order to provide high quality initial approximations in order to significantly accelerate the convergence of iterative solvers at each sample. The method, termed the directed Monte Carlo (DMC) simulation, is developed first for one random variable using the preconditioned conjugate gradient equipped with an initial approximation prediction scheme, and then extended to multiple Gaussian random variable cases by the adoption of a general hyper-spherical transformation. The eigen-poperties of the linear system are also briefly discussed to reveal the suitability of several preconditioning schemes for iterative solvers. Two numerical examples with up to around 6000 DOFs are provided to assess the performance of the proposed solution strategy and associated numerical techniques in terms of computational costs and solution accuracy.

    Analytical trial function method for development of new 8-node plane element based on the variational principle containing Airy stress function
    Engineering Computations, 24 (4) (2010) 442-463
    X.R. Fu, S. Cen, C.F. Li and X.M. Chen

    The purpose of this paper is to propose a novel and simple strategy for construction of hybrid-"stress function" plane element. First, a complementary energy functional, in which the Airy stress function is taken as the functional variable, is established within an element for analysis of plane problems. Second, 15 basic analytical solutions (in global Cartesian coordinates) of the stress function are taken as the trial functions for an 8-node element, and meanwhile, 15 unknown constants are then introduced. Third, according to the principle of minimum complementary energy, the unknown constants can be expressed in terms of the displacements along element edges, which are interpolated by element nodal displacements. Finally, the whole system can be rewritten in terms of element nodal displacement vector. A new hybrid element stiffness matrix is obtained. The resulting 8-node plane element, denoted as analytical trial function (ATF-Q8), possesses excellent performance in numerical examples. Furthermore, some numerical defects, such as direction dependence and interpolation failure, are not found in present model. This paper presents a new strategy for developing finite element models exhibits advantages of both analytical and discrete method.


    Simulating Gaseous Fluids with Low and High Speeds
    Computer Graphics Forum, 28 (7) (2009) 1845-1852   
    G. Yue, C.F. Li, S.M. Hu and B.A. Barsky

    Gaseous fluids may move slowly, as smoke does, or at high speed, such as occurs with explosions. High-speed gas flow is always accompanied by low-speed gas flow, which produces rich visual details in the fluid motion. Realistic visualization involves a complex dynamic flow field with both low and high speed fluid behavior. In computer graphics, algorithms to simulate gaseous fluids address either the low speed case or the high speed case, but no algorithm handles both efficiently. With the aim of providing visually pleasing results, we present a hybrid algorithm that efficiently captures the essential physics of both low- and high-speed gaseous fluids. We model the low speed gaseous fluids by a grid approach and use a particle approach for the high speed gaseous fluids. In addition, we propose a physically sound method to connect the particle model to the grid model. By exploiting complementary strengths and avoiding weaknesses of the grid and particle approaches, we produce some animation examples and analyze their computational performance to demonstrate the effectiveness of the new hybrid method.

    Dynamic behaviors of microtubules in cytosol
    Journal of Biomechanics, 42 (9) (2009) 1270-1274.    
    C.Y. Wang, C.F. Li and S. Adhikari

    Highly anisotropic micro-tubules (MTs) immersed in cytosol are a central part of the cytoskeleton in eukaryotic cells. The dynamic behaviors of an MT-cytosol system are of major interest in biomechanics community. Such a solid-fluid system is characterized by a Reynolds number of the order 10-3 and a slip ionic layer formed at the MT-cytosol interface. In view of these unique features, an orthotropic shell-Stokes flow model with a slip boundary condition has been developed to explore the distinctive dynamic behaviors of MTs in cytosol. Three types of motions have been identified, i.e., (a) undamped and damped torsional vibration, (b) damped longitudinal vibration, and (c) overdamped bending and radial motions. The exponentially decaying bending motion given by the present model is found to be in qualitative agreement with the existing experimental observation [Felgner et al., 1996. Flexural rigidity of microtubules measured with the use of optical tweezers, Journal of Cell Science 109, 509-516].

    Quadrilateral membrane elements with analytical element stiffness matrices formulated by the new quadrilateral area coordinate method (QACM-II)
    International Journal for Numerical Methods in Engineering, 77 (8) (2009) 1172-1200    
    S. Cen, X.M. Chen, C.F. Li and X.R. Fu

    A novel strategy for developing low-order membrane elements with analytical element stiffness matrices is proposed. First, some complete low-order basic analytical solutions for plane stress problems are given in terms of the new quadrilateral area coordinates method (QACM-II). Then, these solutions are taken as the trial functions for developing new membrane elements. Thus, the interpolation formulae for displacement fields naturally possess second-order completeness in physical space (Cartesian coordinates). Finally, by introducing nodal conforming conditions, new 4-node and 5-node membrane elements with analytical element stiffness matrices are successfully constructed. The resulting models, denoted as QAC-ATF4 and QAC-ATF5, have high computational efficiency since the element stiffness matrices are formulated explicitly and no internal parameter is added. These two elements exhibit excellent performance in various bending problems with mesh distortion. It is demonstrated that the proposed strategy possesses advantages of both the analytical and the discrete method, and the QACM-II is a powerful tool for constructing high-performance quadrilateral finite element models.

    A Fourier–Karhunen–Loève discretization scheme for stationary random material properties in SFEM
    International Journal for Numerical Methods in Engineering, 73 (13) (2008) 1942-1965    
    C.F. Li, Y.T. Feng, D.R.J. Owen, D.F. Li and I.M. Davies

    In order to overcome the computational difficulties in Karhunen-Loeve (K-L) expansions of stationary random material properties in stochastic finite element method (SFEM) analysis, a Fourier-Karhunen- Loeve (F-K-L) discretization scheme is developed in this paper, by following the harmonic essence of stationary random material properties and solving a series of specific technical challenges encountered in its development. Three numerical examples are employed to investigate the overall performance of the new discretization scheme and to demonstrate its use in practical SFEM simulations. The proposed F-K-L discretization scheme exhibits a number of advantages over the widely used K-L expansion scheme based on FE meshes, including better computational efficiency in terms of memory and CPU time, convenient a priori error-control mechanism, better approximation accuracy of random material properties, explicit methods for predicting the associated eigenvalue decay speed and geometrical compatibility for random medium bodies of different shapes.


    A fast and accurate analysis of the interacting cracks in linear elastic solids
    Internatioinal Journal of Fracture, 151 (2008) 169-185   
    D.F. Li, C.F. Li, S.Q. Shu, Z.X. Wang and J. Lu

    This paper presents a fast and accurate solution for crack interaction problems in infinite- and half- plane solids. The newsolution is based on themethod of complex potentials developed by Muskhelishvili for the analysis of plane linear elasticity, and it is formulated through three steps. First, the problem is decomposed into a set of basic problems, and for each sub-problem, there is only one crack in the solid. Next, after a crack-dependent conformal mapping, the modified complex potentials associated with the sub-problems are expanded into Laurent's series with unknown coefficients, which in turn provides a mechanism to exactly implement in the form of Fourier series the boundary condition in each sub-problem. Finally, taking into account the crack interaction via a perturbation approach, an iterative algorithm based on fast Fourier transforms (FFT) is developed to solve the unknown Fourier coefficients, and the solution of the whole problem is readily obtained with the superposition of the complex potentials in each sub-problem. The performance of the proposed method is fully investigated by comparing with benchmark results in the literatures, and superb accuracy and efficiency is observed in all situations including patterns where cracks are closely spaced. Also, the new method is able to cope with interactions among a large number of cracks, and this capability is demonstrated by a calculation of effective moduli of an elastic solid with thousands of randomly-spaced cracks.

    Discrete thermal element modelling of heat conduction in particle systems: Basic formulations
    Journal of Computational Physics, 227 (2008) 5072-5089
    Y.T. Feng, K. Han, C.F. Li and D.R.J. Owen

    This paper proposes a novel numerical methodology, termed the discrete thermal element method (DTEM), for the effective modelling of heat conduction in systems comprising a large number of circular particles in 2D cases. Based on an existing analytical integral solution for the temperature distribution over a circular domain subjected to the Neumann boundary condition, a linear algebraic system of thermal conductivity equations for each particle is derived in terms of the average temperatures and the resultant fluxes at the contact zones with its neighboring particles. Thus, each particle is treated as an individual element with the number of (temperature) unknowns equal to the number of particles that it is in contact with. The element thermal conductivity matrix can be very effectively evaluated and is entirely dependent on the characteristics of the contact zones, including the contact positions and contact angles. This new element shares the same form and properties with its conventional thermal finite element counterpart. In particular, the whole solution procedure can follow exactly the same steps as those involved in the finite element analysis. Unlike finite elements or other modern numerical techniques, however, no discretization errors are involved in the discrete thermal elements. The modelling error mainly stems from the assumption made about the heat flux distribution within the contact zones. Based on some theoretical work, an enhanced version is suggested to improve the approximation. The numerical assessment against finite element results indicates that the basic version of DTEM can achieve a very reasonable solution accuracy, while the enhanced version further improves the accuracy to a high level. In addition, thermal resistance phenomena between the contact zones can be readily incorporated into the current modelling framework.


    Fourier representation of random media fields in stochastic finite element modelling
    Engineering Computations, 23 (7) (2006) 794-817
    C.F. Li, Y.T. Feng, D.R.J. Owen and I.M. Davies

    To provide an explicit representation for wide-sense stationary stochastic fields which can be used in stochastic finite element modelling to describe random material properties. Design/methodology/approach - This method represents wide-sense stationary stochastic fields in terms of multiple Fourier series and a vector of mutually uncorrelated random variables, which are obtained by minimizing the mean-squared error of a characteristic equation and solving a standard algebraic eigenvalue problem. The result can be treated as a semi-analytic solution of the Karhunen-Loeve expansion. According to the Karhunen-Loe`ve theorem, a second-order stochastic field can be decomposed into a random part and a deterministic part. Owing to the harmonic essence of wide-sense stationary stochastic fields, the decomposition can be effectively obtained with the assistance of multiple Fourier series. The proposed explicit representation of wide-sense stationary stochastic fields is accurate, efficient and independent of the real shape of the random structure in consideration. Therefore, it can be readily applied in a variety of stochastic finite element formulations to describe random material properties. This paper discloses the connection between the spectral representation theory of wide-sense stationary stochastic fields and the Karhunen-Loe`ve theorem of general second-order stochastic fields, and obtains a Fourier-Karhunen-Loe`ve representation for the former stochastic fields.

    Explicit solution to the stochastic system of linear algebraic equations (a1A1 + a2A2 + ? + amAm)x = b
    Computer Methods in Applied Mechanics and Engineering, 195 (44-47) (2006) 6560-6576.
    C.F. Li, Y.T. Feng and D.R.J. Owen

    This paper presents a novel solution strategy to the stochastic system of linear algebraic equations (a1A1 + a2A2 + ... + amAm)x = b arising from stochastic finite element modelling in computational mechanics, in which ai (i = 1,. . . ,m) denote random variables, Ai (i = 1,. . . ,m) real symmetric deterministic matrices, b a deterministic/random vector and x the unknown random vector to be solved. The system is first decoupled by simultaneously diagonalizing all the matrices Ai via a similarity transformation, and then it is trivial to invert the sum of the diagonalized stochastic matrices to obtain the explicit solution of the stochastic equation system. Unless all the matrices Ai share exactly the same eigen-structure, the joint diagonalization can only be approximately achieved. Hence, the solution is approximate and corresponds to a particular average eigen-structure of the matrix family. Specifically, the classical Jacobi algorithm for the computation of eigenvalues of a single matrix is modified to accommodate multiple matrices and the resulting Jacobi-like joint diagonalization algorithm preserves the fundamental properties of the original version including its convergence and an explicit solution for the optimal Givens rotation angle. Three numerical examples are provided to illustrate the performance of the method.

    SMB: Collision detection based on temporal coherence
    Computer Methods in Applied Mechanics and Engineering, 195 (19-22) (2006) 2252-2269
    C.F. Li, Y.T. Feng and D.R.J. Owen

    The paper presents a novel collision detection algorithm, termed the sort moving boxes (SMB) for large number of moving 2D/3D objects which are represented by their axis-aligned bounding boxes (AABBs). The main feature of the algorithm is the full exploitation of the temporal coherence of the objects exhibited in a dynamic environment. In the algorithm, the AABBs are first projected to each Cartesian axis. The projected intervals on the axes are separately sorted by the diminishing increment sort (DIS) and further divided into subsections. By processing all the intervals within the subsections to check if they overlap, a complete contact list can be built. The SMB is a fast and robust collision detection algorithm, particularly for systems involving a large number of moving AABBs, and also supports for the dynamic insertion and deletion of objects. Its performance in terms of both expected total detection time and memory requirements is proportional to the total number of AABBs, N, and is not influenced by size differences of AABBs, the space size and packing density over a large range up to ten times difference. The only assumption made is that the sorted list at one time step will remain an almost sorted list at the next time step, which is valid for most applications whose movement and deformation of each AABB and the dynamic change of the total number N are approximately continuous.


    Actual morphing: a physics-based approach to blending
    ACM Symposium on Solid Modelling and Applications, 309-314, Genoa, Italy, June 9-11, 2004
    S.M. Hu, C.F. Li and H. Zhang

    When two topologically identical shapes are blended, various possible transformation paths exist from the source shape to the target shape. Which one is the most plausible? Here we propose that the transformation process should obey a quasi-physical law. This paper combines morphing with deformation theory from continuum mechanics. By using strain energy, which reflects the magnitude of deformation, as an objective function, we convert the problem of path interpolation into an unconstrained optimization problem. To reduce the number of variables in the optimization we adopt shape functions, as used in the finite element method (FEM). A point-to-point correspondence between the source and target shapes is naturally established using these polynomial functions plus a distance map.


    Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates. Part II. Verification and numerical investigation
    Composites: Part A Applied Science and Manufacturing, 33 (8) (2002) 1063-1072.
    C.F. Li, N. Hu, J.G. Cheng, H. Fukunaga and H. Sekine

    In Part I of the current work [this issue], we have developed a numerical model for simulating the process of low-velocity impact damage in composite laminates using the finite element method (FEM). This FEM model based on the Mindlin plate element can describe various impact-induced damages and their mutual effects. Some new and effective techniques have also been put forward in that paper, which can significantly increase the computational efficiency. In the current paper, i.e. Part II of the two-part series on the study of impact of composites, we focus on the following two aspects: (a) verification of our numerical model through the comparison with other researchers' results; (b) investigation of the impact-induced damage in the laminated plates using the present numerical model. For the first aspect, some previous experimental data have been adopted for comparison to validate the present numerical model. For the second, we have mainly studied the effects on the impact damage in detail in such aspects as the size of target plate, the boundary conditions of target plate, impact velocity, impactor mass, etc. From these computations, the understanding of the low-velocity impact damage in laminates can be improved. q 2002 Elsevier Science Ltd. All rights reserved.

    Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates. Part I. An FEM numerical model
    Composites: Part A Applied Science and Manufacturing, 33 (8) (2002) 1055-1062
    C.F. Li, N. Hu, Y.J. Yin, H. Sekine and H Fukunaga

    A numerical model for simulating the process of low-velocity impact damage in composite laminates using the finite element method is presented in this paper, i.e. Part I of this two part series on the study of impact. In this model, the 9-node Lagrangian element of the Mindlin plate with consideration of large deformation analysis is employed. To analyze the transient response of the laminated plates, a modified Newmark time integration algorithm previously proposed by the authors is adopted here. We also proved that the impact process between a rigid ball and laminated plates is a stiff system, therefore a kind of A(alpha) stable method has been advocated here to solve the motion equation of the rigid ball. Furthermore, various types of damages including delamination, matrix cracking and fiber breakage, etc. and their mutual influences are modeled and investigated in detail. To overcome the difficulty of numerical oscillation or instability in the analysis of the dynamic contact problem between delaminated layers using the traditional penalty methods, we have employed dynamic spring constraints to simulate the contact effect, which are added to the numerical model by a kind of continuous penalty function. Moreover, an effective technique to calculate the strain energy release rate based on the Mindlin plate model is proposed, which can attain high precision. Finally, some techniques of adaptive analyses have been realized for improving the computational efficiency. Based on this model, a program has been developed for numerically simulating the damage process of cross-ply fiber-reinforced carbon/epoxy composite laminates under lowvelocity impact load. In Part II, this numerical model will be verified by comparing with the experimental results. Also the impact damage will be investigated in detail using this numerical approach.