Dr. Chenfeng Li 1, 2

1 Zienkiewicz Centre for Computational Engineering Post: College of Engineering, Swansea University Bay Campus, Swansea SA1 8EN, UK.
2 Energy Safety Research Institute Tel: +44-01792-602256
College of Engineering Email:


  • Computational Engineering
  • Computer Graphics
  • Civil and Structural Engineering
  • Geo-mechanics, Oil & Gas Reservoir
  • Concrete, Composite, Heterogeneous Materials

  • Computer simulation of deformation, fracture, energy and mass transport
  • Physically-based modelling with visual applications
  • Uncertainty quantification, reliability and risk assessment

    JOURNAL Papers


    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, 2018,
    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 well bore stability of transversely isotropic 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 unsymmetrical 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 modeling 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 unsymmetrical 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 re-meshing
    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 modified hybrid-stress / hybrid-Trefftz 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 image painterly rendering and 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 singlewall 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 micro-tubules 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-Loeve 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.