Three dimensional hybrid-Trefftz stress finite elements for plates and shells are proposed. Two independent fields are approximated: stresses within the element and displacement on their boundary. The required stress field derived from the Papkovitch-Neuber solution of the Navier equation, which a priori satisfies the Trefftz constraint, is generated using homogeneous harmonic polynomials. Restriction on the polynomial degree in the coordinate measured along the thickness direction is imposed to reduce the number of independent terms. Explicit expressions of the corresponding independent polynomials are listed up to the fifth order. Illustrative applications to evaluate displacements and stresses are carried out by hexahedral hybrid-Trefftz stress element models. The hierarchical *p*-refinement and *h*-refinement strategy are exploited in the numerical tests.

Greedy algorithm has been widely adopted for the point selection procedure of radial basis function based mesh deformation. However, in large deformation simulations with thousands of points selected, the greedy point selection will be too expensive and thus become a performance bottleneck. To improve the efficiency of the point selection procedure, a parallel multi-selection greedy method has been developed in this paper. Multiple points are selected at each step to accelerate the convergence speed of the greedy algorithm. In addition, two strategies are presented to determine the specific selecting number. The parallelization of the greedy point selection is realized based on a master-slave model, and a hybrid decomposition algorithm is proposed to address the load imbalance problem. Numerical benchmarks show that both our multi-selection method and the parallelization could obviously improve the point selection efficiency. Specifically, total speedups of 20 and 55 are separately obtained for the 3*D* undulating fish with 10^{6} cell mesh and the 3*D* rotating hydrofoil with 11 million cell mesh. This article is protected by copyright. All rights reserved.

Dynamic relaxation is an iterative method to solve nonlinear systems of equations, which is frequently used for form finding and analysis of structures that undergo large displacements. It is based on the solution of a fictitious dynamic problem where the vibrations of the structure are traced by means of a time integration scheme until a static equilibrium is reached. Fictitious values are used for the mass and damping parameters. Heuristic rules exist to determine these values in such a way that the time integration procedure converges rapidly without becoming unstable. blackCentral to these heuristic rules is the assumption that the highest convergence rate is achieved when the ratio of the highest and lowest eigenfrequency of the structure is minimal.

blackThis short communication shows that all eigenfrequencies become identical when a fictitious mass matrix proportional to the stiffness matrix is used. If, in addition, specific values are used for the fictitious damping parameters and the time integration step, the dynamic relaxation method becomes completely equivalent to the Newton-Raphson method. The Newton-Raphson method can therefore be regarded as a specific form of dynamic relaxation. This insight may help to interpret and improve nonlinear solvers based on dynamic relaxation and/or the Newton-Raphson method. This article is protected by copyright. All rights reserved.

With the fast development of Additive Manufacturing (AM) technology, topology optimization involving multiple materials has received ever increasing attention. Traditionally, this kind of optimization problem is solved within the implicit solution framework using the SIMP (Solid Isotropic Material with Penalization) or level set method. This treatment, however, will inevitably lead to a large number of design variables especially when many types of materials are involved and three-dimensional (3D) problems are considered. This is due to the fact that for each type of material, a corresponding density field/level function defined on the entire design domain must be introduced to describe its distribution. In the present paper, a novel approach for topology optimization with multiple materials is established based on the Moving Morphable Component (MMC) framework. With use of this approach, topology optimization problems with multiple materials can be solved with much less numbers of design variables and degrees of freedom. Numerical examples provided demonstrate the effectiveness of the proposed approach.

The unsymmetric FEM is a promising technique to produce distortion-immune finite elements. In this work, a simple but robust 4-node 12-DOF unsymmetric quadrilateral membrane element is formulated. The test function of this new element is determined by a concise isoparametric-based displacement field which is enriched by the Allman-type drilling degrees of freedom. Meanwhile, a rational stress field, instead of the one in the original unsymmetric formulation, is directly adopted to be the element's trial function. This stress field is obtained based on the analytical solutions of the plane stress/strain problem and the quasi-conforming technique. Thus it can *a priori* satisfy related governing equations. Numerical tests show that the presented new unsymmetric element, named as US-Q4θ, exhibits excellent capabilities in predicting results of both displacement and stress, in most cases, superior to other existing 4-node element models. In particular, it can still work very well in severely distorted meshes even when the element shape deteriorates into concave quadrangle or degenerated triangle.

This paper presents a new concurrent simulation approach to couple isogeometric analysis (IGA) with the meshfree method for studying of crack problems. In the present method, the overall physical domain is divided into two sub-domains which are formulated with the IGA and meshfree method, respectively. In the meshfree sub-domain, the moving least squares (MLS) shape function is adopted for the discretization of the area around crack tips, and the IGA sub-domain is adopted in the remaining area. Meanwhile, the interface region between the two sub-domains is represented by coupled shape functions. The resulting shape function, which comprises both IGA and meshfree shape functions, satisfies the consistency condition, thus ensuring convergence of the method. Moreover, the meshfree shape functions augmented with the enriched basis functions to predict the singular stress fields near a crack tip are presented. The proposed approach is also applied to simulate the crack propagation under a mixed-mode condition. Several numerical examples are studied to demonstrate the utility and robustness of the proposed method. This article is protected by copyright. All rights reserved.

The stability and reflection-transmission properties of the bipenalty method are studied in application to explicit finite element analysis of one-dimensional contact-impact problems. It is known that the standard penalty method, where an additional stiffness term corresponding to contact boundary conditions is applied, attacks the stability limit of finite element model. Generally, the critical time step size rapidly decreases with increasing penalty stiffness. Recent comprehensive studies have shown that the so-called bipenalty technique, using mass penalty together with standard stiffness penalty, preserves the critical time step size associated to contact-free bodies. In this paper, the influence of the penalty ratio (ratio of stiffness and mass penalty parameters) on stability and reflection-transmission properties in one-dimensional contact-impact problems using the same material and mesh size for both domains is studied. The paper closes with numerical examples which demonstrate the stability and reflection-transmission behaviour of the bipenalty method in one-dimensional contact-impact and wave propagation problems of homogeneous materials. This article is protected by copyright. All rights reserved.

The Koiter-Newton (KN) method is a combination of local multi-mode polynomial approximations inspired by Koiter's initial post-buckling theory and global corrections using the standard Newton's method. In the original formulation, the local polynomial approximation, called a reduced order model, is used to make significantly more accurate predictions compared to the standard linear prediction used in conjunction with Newton's method. The correction to the exact equilibrium path relied exclusively on Newton-Raphson method using the full model.

In this paper, we proposed a modified Newton-type Koiter-Newton method to trace the geometrically nonlinear response of structures . The developed predictor-corrector strategy is applied to each predicted solution of the reduced order model. The reduced order model can be used also in the correction phase, and the exact full nonlinear model is applied only to calculate force residuals. Remainder terms in both the displacement expansion and the reduced order model are well considered and constantly updated during correction. The same augmented FEM system is used for both the construction of the reduced order model and the iterations for correction. Hence, the developed method can be seen as a particular modified Newton method with a constant iteration matrix over the single Koiter-Newton step. This significantly reduces the computational cost of the method. As a side product, the method has better error control, leading to more robust step size adaptation strategies. Numerical results demonstrate the effectiveness of the method in treating nonlinear buckling problems. This article is protected by copyright. All rights reserved.

We introduce a coupled finite and boundary element formulation for acoustic scattering analysis over thin shell structures. A triangular Loop subdivision surface discretisation is used for both geometry and analysis fields. The Kirchhoff-Love shell equation is discretised with the finite element method and the Helmholtz equation for the acoustic field with the boundary element method. The use of the boundary element formulation allows the elegant handling of infinite domains and precludes the need for volumetric meshing. In the present work the subdivision control meshes for the shell displacements and the acoustic pressures have the same resolution. The corresponding smooth subdivision basis functions have the *C*^{1} continuity property required for the Kirchhoff-Love formulation and are highly efficient for the acoustic field computations. We redverify the proposed isogeometric formulation through a closed-form solution of acoustic scattering over a thin shell sphere. Furthermore, we demonstrate the ability of the proposed approach to handle complex geometries with arbitrary topology that provides an integrated isogeometric design and analysis workflow for coupled structural-acoustic analysis of shells. This article is protected by copyright. All rights reserved.

In this paper, we present a new method for inserting several triangulated surfaces into an existing tetrahedral mesh generated by the meccano method. The result is a conformal mesh where each inserted surface is approximated by a set of faces of the final tetrahedral mesh. First, the tetrahedral mesh is refined around the inserted surfaces to capture their geometric features. Second, each immersed surface is approximated by a set of faces from the tetrahedral mesh. Third, following a novel approach the nodes of the approximated surfaces are mapped to the corresponding immersed surface. Fourth, we untangle and smooth the mesh by optimizing a regularized shape distortion measure for tetrahedral elements in which we move all the nodes of the mesh, restricting the movement of the edge and surface nodes along the corresponding entity they belong to. The refining process allows approximating the immersed surface for any initial meccano tetrahedral mesh. Moreover, the proposed projection method avoids computational expensive geometric projections. Finally, the applied simultaneous untangling and smoothing process delivers a high-quality mesh and ensures that the immersed surfaces are interpolated. Several examples are presented to assess the properties of the proposed method. This article is protected by copyright. All rights reserved.

In this paper, the coupling of the improved interpolating element-free Galerkin (IIEFG) method and the variable-order infinite acoustic wave envelope element (WEE) method is studied. A coupled IIEFG-WEE method for computing sound radiation is proposed to make use of their advantages while evading their disadvantages. The coupling is achieved by constructing the hybrid shape function of continuity and compatibility on the interface between the IIEFG and WEE domains. In the IIEFG domain, the improved interpolating moving least-squares (IIMLS) method is employed to form the shape functions satisfying the Kronecker delta condition while nonsingular weight functions can be used. The impacts of the size of the influence domain and the shape parameter on the performance of this coupled method are investigated. The numerical results show that the coupled IIEFG-WEE method can take full advantage of both the IIEFG and WEE methods, and that it not only can achieve higher accuracy but also has a faster convergence speed than the conventional method of the finite element coupled with the WEE. The experimental results show that the method is very flexible for acoustic radiation prediction in the infinite domain. This article is protected by copyright. All rights reserved.

This study proposes a framework for uncertainty analysis by incorporating explicit numerical solutions of governing equations for flood wave propagation with the expectation operator. It aims at effectively evaluating the effect of variations in initial and boundary conditions on the estimation of flood waves. Spatiotemporal semi-variogram models are employed to quantify the correlation of the variables in time and space. The 1-D nonlinear kinematic wave equation for the overland flow (named *EVO_NS_KWE*) is applied in the model development. Model validation is made by comparison with the Monte Carlo simulation (*MCS*) model in the calculation of statistical properties of model outputs (i.e. flow depths), that is, the mean, standard deviation and coefficient of variation (*CV*). The results from the model validation show that the *EVO_NS_KWE* model can produce excellent approximations of the mean and less satisfactory approximations of the standard deviation and *CV* compared with those obtained using the *MCS* model. It concludes that the uncertainties of flow depths in the domain are significantly affected by variations in the boundary condition. Future, application of the *EVO_NS_KWE* model enable the evaluation of uncertainty in model outputs induced by the initial and boundary condition subject to uncertainty, and will also provide corresponding probabilistic information for risk quantification method.

A methodology for the calculation of gradients with respect to design parameters in general Fluid-Structure Interaction problems is presented. It is based on fixed-point iterations on the adjoint variables of the coupled system using Algorithmic Differentiation. This removes the need for the construction of the analytic Jacobian for the coupled physical problem, which is the usual limitation for the computation of adjoints in most realistic applications. The formulation is shown to be amenable to partitioned solution methods for the adjoint equations. It also poses no restrictions to the nonlinear physics in either the fluid or structural field, other than the existence of a converged solution to the primal problem from which to compute the adjoints. We demonstrate the applicability of this procedure and the accuracy of the computed gradients on coupled problems involving viscous flows with geometrical and material non-linearities in the structural domain. This article is protected by copyright. All rights reserved.

Two stable approximation space configurations are treated for the mixed finite element method for elliptic problems based on curved meshes. Their choices are guided by the property that, in the master element, the image of the flux space by the divergence operator coincides with the potential space. By using static condensation, the sizes of global condensed matrices, which are proportional to the dimension of border fluxes, are the same in both configurations. The meshes are composed of different topologies (tetrahedra, hexahedra or prisms). Simulations using asymptotically affine uniform meshes, exactly fitting a spherical like region, and constant polynomial degree distribution k, show L2-errors of order k+1 or k+2 for the potential variable, while keeping order k+1 for the flux in both configurations. The first case corresponds to RT(k) and BDFM(k+1) spaces for hexahedral and tetrahedral meshes, respectively, but holding for prismatic elements as well. The second case, further incrementing the order of approximation of the potential variable, holds for the three element topologies. The case of hp-adaptive meshes is considered for a problem modelling a porous media flow around a cylindrical horizontal well with elliptical drainage area. The effect of parallelism and static condensation in CPU time reduction is illustrated. This article is protected by copyright. All rights reserved.

Beam-to-beam contact can be dually described depending on the scale as contact between curves or contact between surfaces. The surface-to-surface algorithm requires a computationally expensive continuum description of the beams. The curve-to-curve (CTC) contact algorithm is generated from the corresponding closest point projection (CPP) procedure determining a minimal distance between curves and involves less expensive beam finite elements.

The existence and uniqueness criterion of the CPP procedure formulated already in earlier publications is analyzed in detail, especially in the case of multiple solutions for “parallel curves“. It is shown that the CTC contact algorithm cannot be applied for general parallel curves, which can be generated in an arbitrary range of angles between the tangent lines. A new Curve-To-Solid Beam contact algorithm is here developed. This algorithm requires a special Solid Beam finite element and is combining the advantages of surface contact together with the flexibility of the application to beam-to-beam contact.

Further, a special case with ”parallel tangent“ vectors is considered from the classical Hertzian contact problem. It is found numerically that the corresponding length of a contact zone is well correlated with the Hertz theory. A possible application of the developed strategy for wire ropes is analyzed numerically. This article is protected by copyright. All rights reserved.

In this paper, we propose a novel unfitted finite element method for the simulation of multiple body contact. The computational mesh is generated independently of the geometry of the interacting solids, which can be arbitrarily complex. The key novelty of the approach is the combination of elements of the CutFEM technology, namely the enrichment of the solution field via the definition of overlapping fictitious domains with a dedicated penalty-type regularisation of discrete operators, and the LaTIn hybrid-mixed formulation of complex interface conditions. Furthermore, the novel P1-P1 discretisation scheme that we propose for the unfitted LaTIn solver is shown to be stable, robust and optimally convergent with mesh refinement. Finally, the paper introduces a high-performance 3D level-set/CutFEM framework for the versatile and robust solution of contact problems involving multiple bodies of complex geometries, with more than two bodies interacting at a single point. This article is protected by copyright. All rights reserved.

This work concerns an application of the Tsallis entropy to homogenization problem of the fiber-reinforced and also of the particle-filled composites with random material and geometrical characteristics. Calculation of the effective material parameters is done with two alternative homogenization methods – the first is based upon the deformation energy of the Representative Volume Element (RVE) subjected to the few specific deformations while the second uses explicitly the so-called homogenization functions determined under periodic boundary conditions imposed on this RVE. Probabilistic homogenization is made with the use of three concurrent non-deterministic methods, namely Monte-Carlo simulation, iterative generalized stochastic perturbation technique as well as the semi-analytical approach. The last two approaches are based on the Least Squares Method with polynomial basis of the statistically optimized order – this basis serves for further differentiation in the tenth order stochastic perturbation technique, while semi-analytical method uses it in probabilistic integrals. These three approaches are implemented all as the extensions of the traditional Finite Element Method (FEM) with contrastively different mesh sizes and they serve in computations of Tsallis entropies of the homogenized tensor components as the functions of input coefficient of variation.

Fractures tend to propagate along the least resistance paths, and homogeneous-based models may not be able to reliably predict the true crack paths, as they are not capable of capturing nonlinearities and local damage induced by local inhomogeneity. This paper presents a stochastic numerical modelling framework for simulating fracturing in natural heterogeneous materials. Fracture propagation is modelled using Francfort and Marigo's variational theory, and randomness in the material properties is introduced by random field principle. A computational strategy on the basis of nonlinear dimensionality reduction framework is developed that maps domain of spatially variable properties of the materials to a low-dimensional space. This strategy allows us to predict the most probable fracture patterns leading to failure by an optimisation algorithm. The reliability and performance of the developed methodology are examined through simulation of experimental case studies and comparison of predictions with measured data.

The present work addresses a multiscale framework for Fast-Fourier-Transform based computational homogenization. The framework considers the scale-bridging between microscopic and macroscopic scale. While the macroscopic problem is discretized with finite elements, the microscopic problems are solved by means of Fast Fourier Transforms on periodic representative volume elements (RVE). In such multiscale scenario, the computation of the *effective properties* of the microstructure is crucial. While effective quantities in terms of stresses and deformations can be computed from surface integrals along the boundary of the RVE, the computation of the associated moduli is not straight-forward. The key contribution of the present paper is thus the derivation and implementation of an *algorithmically consistent macroscopic tangent operator* which directly resembles the effective moduli of the microstructure. The macroscopic tangent is derived by means of the classical Lippmann-Schwinger equation and can be computed from a simple system of linear equations. This is performed through an efficient Fast-Fourier-Transform based approach along with a conjugate gradient solver. The viability and efficiency of the method is demonstrated for a number of two- and three-dimensional boundary value problems incorporating linear and nonlinear elasticity as well as viscoelastic material response. This article is protected by copyright. All rights reserved.

This paper deals with model-order reduction of parametric partial differential equations (PPDE). More specifically, we consider the problem of finding a good approximation subspace of the solution manifold of the PPDE when only partial information on the latter is available. We assume that two sources of information are available: *i)* a “rough” prior knowledge, taking the form of a manifold containing the target solution manifold; *ii)* partial linear measurements of the solutions of the PPDE (the term partial refers to the fact that observation operator cannot be inverted). We provide and study several tools to derive good approximation subspaces from these two sources of information. We first identify the best worst-case performance achievable in this setup and propose simple procedures to approximate the corresponding optimal approximation subspace. We then provide, in a simplified setup, a theoretical analysis relating the achievable reduction performance to the choice of the observation operator and the prior knowledge available on the solution manifold. This article is protected by copyright. All rights reserved.

We present 2 adaptive refinement techniques, namely, adaptive local refinement and adaptive hierarchical refinement, for isogeometric analysis. An element-wise point of view is adopted, exploiting Bézier extraction, which facilitates the implementation of adaptive refinement in isogeometric analysis. Locally refined and hierarchical T-splines are used for the description of the geometry as well as for the approximation of the solution space in the analysis. The refinement is conducted with the aid of a subdivision operator, which is computed by again exploiting the Bézier extraction operator. The concept and algorithm of an element-based adaptive isogeometric analysis are illustrated. Numerical examples are given to examine the accuracy, the convergence, and the condition number.

This paper presents a framework for r-adaptive quasi-static configurational force (CF) brittle crack propagation, cast within a discontinuous Galerkin (DG) symmetric interior penalty (SIPG) finite element scheme. Cracks are propagated in discrete steps, with a staggered algorithm, along element interfaces, which align themselves with the predicted crack propagation direction. The key novelty of the work is the exploitation of the DG face stiffness terms existing along element interfaces to propagate a crack in a mesh-independent r-adaptive quasi-static fashion, driven by the CF at the crack tip. This adds no new degrees of freedom to the data structure. Additionally, as DG methods have element-specific degrees of freedom, a geometry-driven p-adaptive algorithm is also easily included allowing for more accurate solutions of the CF on a moving crack front. Further, for nondeterminant systems, we introduce an average boundary condition that restrains rigid body motion leading to a determinant system. To the authors' knowledge, this is the first time that such a boundary condition has been described. The proposed formulation is validated against single and multiple crack problems with single- and mixed-mode cracks, demonstrating the predictive capabilities of the method.

The extended finite element method was introduced in 1999 to treat problems involving discontinuities with no or minimal remeshing through appropriate enrichment functions. This enables elements to be split by a discontinuity, strong or weak, and hence requires the integration of discontinuous functions or functions with discontinuous derivatives over elementary volumes. A variety of approaches have been proposed to facilitate these special types of numerical integration, which have been shown to have a large impact on the accuracy and the convergence of the numerical solution. The smoothed extended finite element method (XFEM), for example, makes numerical integration elegant and simple by transforming volume integrals into surface integrals. However, it was reported in the literature that the strain smoothing is inaccurate when non-polynomial functions are in the basis. In this paper, we investigate the benefits of a recently developed Linear smoothing procedure which provides better approximation to higher-order polynomial fields in the basis. Some benchmark problems in the context of linear elastic fracture mechanics are solved and the results are compared with existing approaches. We observe that the stress intensity factors computed through the proposed linear smoothed XFEM is more accurate than that obtained through smoothed XFEM.

In this paper, a level-set-based method is presented to deal with the multi-material topology optimization of compliant mechanisms with stress constraints. A novel stress-based multi-material topology optimization model of compliant mechanisms is proposed. In this model, the multi-material level set topology description model and the separable stress interpolation scheme are adopted. The weighted sum method is used to deal with the multi-objective optimization of the output displacement and compliance of compliant mechanisms. The penalty of stresses is also considered in the objective function to control the local stress level in different materials. To solve the optimization problem, the parametric level set method is employed, and the sensitivity analysis is conducted. Application of the method is demonstrated by 2 numerical examples. Results show that the multi-material structures without undesirable de facto hinges can be obtained. The output displacement and compliance of the compliant mechanisms are optimized, and stress constraints in different materials are simultaneously satisfied.

This paper introduces multivariate input-output models to predict the errors and bases dimensions of local parametric Proper Orthogonal Decomposition reduced-order models. We refer to these mappings as the multivariate predictions of local reduced-order model characteristics (MP-LROM) models. We use Gaussian processes and artificial neural networks to construct approximations of these multivariate mappings. Numerical results with a viscous Burgers model illustrate the performance and potential of the machine learning-based regression MP-LROM models to approximate the characteristics of parametric local reduced-order models. The predicted reduced-order models errors are compared against the multifidelity correction and reduced-order model error surrogates methods predictions, whereas the predicted reduced-order dimensions are tested against the standard method based on the spectrum of snapshots matrix. Since the MP-LROM models incorporate more features and elements to construct the probabilistic mappings, they achieve more accurate results. However, for high-dimensional parametric spaces, the MP-LROM models might suffer from the curse of dimensionality. Scalability challenges of MP-LROM models and the feasible ways of addressing them are also discussed in this study.

The paper deals with the use of model order reduction within a posteriori error estimation procedures in the context of the finite element method. More specifically, it focuses on the constitutive relation error concept, which has been widely used over the last 40 years for FEM verification of computational mechanics models. A technical key-point when using constitutive relation error is the construction of admissible fields, and we propose here to use the proper generalized decomposition to facilitate this task. In addition to making the implementation into commercial FE software easier, it is shown that the use of proper generalized decomposition enables to optimize the verification procedure and to get both accurate and reasonably expensive upper bounds on the discretization error. Numerical illustrations are presented to assess the performance of the proposed approach.

A detailed and systematic analysis is performed on the local and global properties of the recently developed harmonic polynomial cell (HPC) method, a very accurate and efficient field solver for problems governed by the Laplace equation. At the local cell level, a simple rule is identified for the proper choice of harmonic polynomials in the local representation of the velocity potential in cells with symmetry properties. The local solution error, its convergence rate, its dependence on the cell topology, its distribution inside the cell, and its features across cells with different dimensions are carefully examined with relevant findings for HPC numerical implementations. At the global level, the error convergence rate is analytically estimated in terms of error contributions from the boundary conditions and from inside the liquid domain. In most cases, the error associated with boundary conditions dominates the global error. In order to minimize it, Quadtree grid strategies or high-order local expressions of the velocity potential are proposed for cells near critical boundary portions. To model accurately the boundary conditions on rigid or deformable surfaces with generic geometries, 3 different grid strategies are proposed by adopting concepts of immersed boundary method and overlapping grids. They are comparatively studied for a circular rigid cylinder in infinite fluid and for the propagation of a free-surface wave. Then, an immersed boundary strategy, using numerical choices suggested in this paper, is successfully compared against a fully nonlinear boundary element method for the case of a surface-piercing circular cylinder heaving in otherwise calm water.

The paper presents a novel model order reduction technique for large-scale linear parameter-varying (LPV) systems. The approach is based on decoupling the original dynamics into smaller dimensional LPV subsystems that can be independently reduced by parameter-varying reduction methods. The decomposition starts with the construction of a modal transformation that separates the modal subsystems. Hierarchical clustering is applied then to collect the dynamically similar modal subsystems into larger groups. The resulting parameter-varying subsystems are then independently reduced. This approach substantially differs from most of the previously proposed LPV model reduction techniques, since it performs manipulations on the LPV model itself, instead of on a set of linear time-invariant models defined at fixed scheduling parameter values. Therefore, the interpolation, which is often a challenging part in reduction techniques, is inherently solved. The applicability of the developed algorithm is thoroughly investigated and demonstrated by numerical case studies.

This paper focuses on topology optimization for linear transient thermo-mechanical problems. The latter are, for example, encountered for extreme precision tools and instrumentation. Due to the transient nature, a standard adjoint sensitivity analysis will result in a backward transient analysis for the adjoint variables, leading to both storage and computational inefficiencies. A method is proposed that rigorously eliminates the backward transient integration for the adjoint sensitivity analysis. At the basis is a model-order reduction technique, which relies on a reduced thermal modal basis combined with static correction. The modal amplitudes can be readily obtained semi-analytically using simple convolutions. This accurate but reduced-order model is the starting point for an adjoint sensitivity analysis. Via a tactic selection of adjoint variables, the backward transient analysis for the adjoints is completely eliminated, whereas computational efficiency and consistency are maintained. The effectiveness of the resulting adjoint sensitivities and their application in topology optimization are demonstrated on the basis of several test examples.

We present a three-dimensional vector level set method coupled to a recently developed stable extended finite element method (XFEM). We further investigate a new enrichment approach for XFEM adopting discontinuous linear enrichment functions in place of the asymptotic near-tip functions. Through the vector level set method, level set values for propagating cracks are obtained via simple geometrical operations, eliminating the need for solution of differential evolution equations. The first XFEM variant ensures optimal convergence rates by means of geometrical enrichment, ie, the use of enriched elements in a fixed volume around the crack front, without giving rise to conditioning problems. The linear enrichment approach, significantly simplifies implementation and reduces the computational cost associated with numerical integration, while providing nonoptimal convergence rates similar to standard finite elements. The 2 dicretization schemes are tested for different benchmark problems, and their combination to the vector level set method is verified for nonplanar crack propagation problems.

A novel Lagrangian gradient smoothing method (L-GSM) is developed to solve “solid-flow” (flow media with material strength) problems governed by Lagrangian form of Navier-Stokes equations. It is a particle-like method, similar to the smoothed particle hydrodynamics (SPH) method but without the so-called tensile instability that exists in the SPH since its birth. The L-GSM uses gradient smoothing technique to approximate the gradient of the field variables, based on the standard GSM that was found working well with Euler grids for general fluids. The Delaunay triangulation algorithm is adopted to update the connectivity of the particles, so that supporting neighboring particles can be determined for accurate gradient approximations. Special techniques are also devised for treatments of 3 types of boundaries: no-slip solid boundary, free-surface boundary, and periodical boundary. An advanced GSM operation for better consistency condition is then developed. Tensile stability condition of L-GSM is investigated through the von Neumann stability analysis as well as numerical tests. The proposed L-GSM is validated by using benchmarking examples of incompressible flows, including the Couette flow, Poiseuille flow, and 2D shear-driven cavity. It is then applied to solve a practical problem of solid flows: the natural failure process of soil and the resultant soil flows. The numerical results are compared with theoretical solutions, experimental data, and other numerical results by SPH and FDM to evaluate further L-GSM performance. It shows that the L-GSM scheme can give a very accurate result for all these examples. Both the theoretical analysis and the numerical testing results demonstrate that the proposed L-GSM approach restores first-order accuracy unconditionally and does not suffer from the tensile instability. It is also shown that the L-GSM is much more computational efficient compared with SPH, especially when a large number of particles are employed in simulation.

This paper presents a (higher-order) finite element approach for the simulation of heat diffusion and thermoelastic deformations in NC-milling processes. The inherent continuous material removal in the process of the simulation is taken into account via continuous removal-dependent refinements of a paraxial hexahedron base-mesh covering a given workpiece. These refinements rely on isotropic bisections of these hexahedrons along with subdivisions of the latter into tetrahedrons and pyramids in correspondence to a milling surface triangulation obtained from the application of the marching cubes algorithm. The resulting mesh is used for an element-wise defined characteristic function for the milling-dependent workpiece within that paraxial hexahedron base-mesh. Using this characteristic function, a (higher-order) fictitious domain method is used to compute the heat diffusion and thermoelastic deformations, where the corresponding ansatz spaces are defined for some hexahedron-based refinement of the base-mesh. Numerical experiments compared to real physical experiments exhibit the applicability of the proposed approach to predict deviations of the milled workpiece from its designed shape because of thermoelastic deformations in the process.

We study time integration schemes for discontinuous Galerkin discretizations of problems with moving immersed interfaces. Two approaches have been discussed in literature so far: splitting schemes and space-time methods. Splitting schemes are cheap and easy to implement, but are non-conservative, inherently limited to low orders of accuracy, and require extremely small time steps. Space-time methods, on the other hand, are conservative, allow for large time steps, and generalize to arbitrary orders of accuracy. However, these advantages come at the expense of a severe growth the systems to be solved. Within this work, we present a generic strategy that combines the advantages of both concepts by evaluating numerical fluxes in a moving reference frame and by making use of a conservative cell-agglomeration strategy. We study the performance of this strategy in combination with backward-difference-formulas and explicit Runge-Kutta schemes in the context of the scalar transport equation, the Burgers equation, and the heat equation. Our results indicate that higher order spatial and temporal convergence rates can be achieved without increasing the size of the systems to be solved.

The paper describes the main features of an automatic and three-dimensional Cartesian mesher specifically designed for compressible inviscid/viscous flow solvers based on an immersed boundary technique. The development of a meshing tool able of conducting non-isotropic cell refinements is a very tiresome task. The major difficulty is to imagine, at the pre-design phase, a light but flexible data management, which minimizes the memory and CPUs' resources. In particular, the embedded geometry has to be detected by means of a fast and robust tagging procedure. Cells in proximity of the wall have to be refined in a proper way to adequately solve large flow gradients. Smooth variation of mesh density between differently refined zones has to be guaranteed to increase the flow solver robustness. A procedure to obtain accurate data on the geometry surfaces should be foresee. Here, a robust algorithm is developed to reconstruct a surface triangulation starting from the intersection points among volume cells and the geometry surfaces. The paper attempts to address all the above issues to help the readers in designing their own tools and suggesting them a way forward.

In this paper, conserving time-stepping algorithms for frictionless and full stick friction dynamic contact problems are presented. Time integration algorithms for frictionless and full stick friction dynamic contact problems have been designed to preserve the conservation of key discrete properties satisfied at the continuum level. Energy and energy-momentum–preserving algorithms for frictionless and full stick friction dynamic contact problems, respectively, have been designed and implemented within the framework of the direct elimination method, avoiding the drawbacks linked to the use of penalty-based or Lagrange multipliers methods. An assessment of the performance of the resulting formulation is shown in a number of selected and representative numerical examples, under full stick friction and slip frictionless contact conditions. Conservation of key discrete properties exhibited by the time-stepping algorithm is shown.

We introduce a new methodology for modeling problems with both weak and strong discontinuities independently of the finite element discretization. At variance with the eXtended/Generalized Finite Element Method (X/GFEM), the new method, named the Discontinuity-Enriched Finite Element Method (DE-FEM), adds enriched degrees of freedom only to nodes created at the intersection between a discontinuity and edges of elements in the mesh. Although general, the method is demonstrated in the context of fracture mechanics, and its versatility is illustrated with a set of traction-free and cohesive crack examples. We show that DE-FEM recovers the same rate of convergence as the standard FEM with matching meshes, and we also compare the new approach to X/GFEM.

This paper presents a high fidelity homogenization method for periodically layered composite structures that accounts for plasticity in the matrix material and quasi-brittle damage in the reinforcing layers, combined with strong geometrical nonlinearities. A set of deliberately chosen internal kinematic variables results in a rigorous representation of the kinematics of the 2 constituents, which in turn allows for complex constitutive laws per constituent to be used directly in the formulation. The model accounts for hyper-elastoplastic behavior in the matrix phase and hyper-elastic behavior in the reinforcement as well as for the bending stiffness of the reinforcement layers. Additionally to previously proposed models, the present method includes Lemaitre-type damage for the reinforcement, making it applicable to a wider range of engineering applications. The capability of the proposed method in representing the combined effect of plasticity, damage, and buckling at microlevel within a homogenized setting is demonstrated by means of direct comparisons to a reference discrete model.

A methodology is proposed for the efficient solution of probabilistic nonconvex constrained optimization problems with uncertain. Statistical properties of the underlying stochastic generator are characterized from an initial statistical sample of function evaluations. A diffusion manifold over the initial set of data points is first identified and an associated basis computed. The joint probability density function of this initial set is estimated using a kernel density model and an Itô stochastic differential equation (ISDE) constructed with this model as its invariant measure. This ISDE is adapted to fluctuate around the manifold yielding additional joint realizations of the uncertain parameters, design variables, and function values, which are obtained as solutions of the ISDE. The expectations in the objective function and constraints are then accurately evaluated without performing additional function evaluations. The methodology brings together novel ideas from manifold learning and stochastic Hamiltonian dynamics to tackle an outstanding challenge in stochastic optimization. Three examples are presented to highlight different aspects of the proposed methodology.

Laminated piezocomposite energy harvesters (LAPEHs) are multilayer arrangements of piezoelectric and nonpiezoelectric materials. Multiple materials and physics, and dynamic analysis need to be considered in their design. Usually these devices are designed for harmonic excitation; however, they are subjected to other types of excitations. Thus, a novel topology optimization formulation is developed for designing LAPEHs that considers a combination of harmonic and transient optimization problems with the aim of designing the so-called “multi-entry" devices in which the power generated is the same for different types of excitation. LAPEHs are modeled by the finite element method, and the material model used for the piezoelectric layer is based on penalization and polarization model who controls material distribution and corresponding polarization. To optimize the RLC circuit, a novel linear interpolation model of coupled electrical impedance is also introduced to consider different magnitudes of the coupled impedance. The topology optimization problem seeks to maximize the active power generated by the LAPEH at its RLC circuit, to minimize its response time measured as the slope of the power versus time curve, and to maximize its stiffness. Numerical examples are shown to illustrate the potential of the method.

The conventional SPH method has been frustrated with the low accuracy originated from the particle inconsistency. The finite particle method (FPM) is an improved SPH method with a higher-order accuracy. However, the numerical accuracy and stability of FPM depend on the uniformity of particles. Facing the disorder particles, the conventional FPM may suffer from an ill-conditioned corrective matrix and is difficult to conduct longtime simulations. A popular and robust particle regularization scheme is the so-called particle shifting technology (PST), which can effectively enhance the accuracy and stability of particle methods. In the context of FPM, PST is analyzed and discussed, and a modified PST (MPST) is proposed. Modified PST saves great amount of computational cost with respect to the conventional PST and acquires better features of accuracy and stability. Finally, the coupled FPM method by combining MPST and δ-SPH is developed to simulate a series of viscous flows. The numerical results indicate that MPST is effective in improving accuracy, stability, and efficiency of PST, and the coupled FPM is shown to be robust for simulating viscous flows and has a higher accuracy and stability.

A novel approach, referred to as the homotopy stochastic finite element method, is proposed to solve the eigenvalue problem of a structure associated with some amount of uncertainty based on the homotopy analysis method. For this approach, an infinite multivariate series of the involved random variables is proposed to express the random eigenvalue or even a random eigenvector. The coefficients of the multivariate series are determined using the homotopy analysis method. The convergence domain of the derived series is greatly expanded compared with the Taylor series due to the use of an approach function of the parameter *h*. Therefore, the proposed method is not limited to random parameters with small fluctuation. However, in practice, only single-variable and double-variable approximations are employed to simplify the calculation. The numerical examples show that with a suitable choice of the auxiliary parameter *h*, the suggested approximations can produce very accurate results and require reduced or similar computational efforts compared with the existing methods.

Modelling the entire ductile fracture process remains a challenge. On the one hand, continuous damage models succeed in capturing the initial diffuse damage stage but are not able to represent discontinuities or cracks. On the other hand, discontinuous methods, as the cohesive zones, which model the crack propagation behaviour, are suited to represent the localised damaging process. However, they are unable to represent diffuse damage. Moreover, most of the cohesive models do not capture triaxiality effect. In this paper, the advantages of the two approaches are combined in a single damage to crack transition framework. In a small deformation setting, a nonlocal elastic damage model is associated with a cohesive model in a discontinuous Galerkin finite element framework. A cohesive band model is used to naturally introduce a triaxiality-dependent behaviour inside the cohesive law. Practically, a numerical thickness is introduced to recover a 3D state, mandatory to incorporate the in-plane stretch effects. This thickness is evaluated to ensure the energy consistency of the method and is not a new numerical parameter. The traction-separation law is then built from the underlying damage model. The method is numerically shown to capture the stress triaxiality effect on the crack initiation and propagation.

Fictitious domain methods are attractive for shape optimization applications, since they do not require deformed or regenerated meshes. A recently developed such method is the CutFEM approach, which allows crisp boundary representations and for which uniformly well-conditioned system matrices can be guaranteed. Here, we investigate the use of the CutFEM approach for acoustic shape optimization, using as test problem the design of an acoustic horn for favorable impedance-matching properties. The CutFEM approach is used to solve the Helmholtz equation, and the geometry of the horn is implicitly described by a level-set function. To promote smooth algorithmic updates of the geometry, we propose to use the nodal values of the Laplacian of the level-set function as design variables. This strategy also improves the algorithm's convergence rate, counteracts mesh dependence, and, in combination with Tikhonov regularization, controls small details in the optimized designs. An advantage with the proposed method is that the exact derivatives of the discrete objective function can be expressed as boundary integrals, as opposed to when using a traditional method that uses mesh deformations. The resulting horns possess excellent impedance-matching properties and exhibit surprising subwavelength structures, not previously seen, which are possible to capture due to the fixed mesh approach.

Several cases of nonlinear wave propagation are studied numerically in two dimensions within the framework of potential flow. The Laplace equation is solved with the harmonic polynomial cell (HPC) method, which is a field method with high-order accuracy. In the HPC method, the computational domain is divided into overlapping cells. Within each cell, the velocity potential is represented by a sum of harmonic polynomials. Two different methods denoted as immersed boundary (IB) and multigrid (MG) are used to track the free surface. The former treats the free surface as an IB in a fixed Cartesian background grid, while the latter uses a free-surface fitted grid that overlaps with a Cartesian background grid. The simulated cases include several nonlinear wave mechanisms, such as high steepness and shallow-water effects. For one of the cases, a numerical scheme to suppress local wave breaking is introduced. Such scheme can serve as a practical mean to ensure numerical stability in simulations where local breaking is not significant for the result. For all the considered cases, both the IB and MG method generally give satisfactory agreement with known reference results. Although the two free-surface tracking methods mostly have similar performance, some differences between them are pointed out. These include aspects related to modeling of particular physical problems as well as their computational efficiency when combined with the HPC method.

In goal-oriented adaptivity, the error in the quantity of interest is represented using the error functions of the direct and adjoint problems. This error representation is subsequently bounded above by element-wise error indicators that are used to drive optimal refinements. In this work, we propose to replace, in the error representation, the adjoint problem by an alternative operator. The main advantage of the proposed approach is that, when judiciously selecting such alternative operator, the corresponding upper bound of the error representation becomes sharper, leading to a more efficient goal-oriented adaptivity.

While the method can be applied to a variety of problems, we focus here on two- and three-dimensional (2-D and 3-D) Helmholtz problems. We show via extensive numerical experimentation that the upper bounds provided by the alternative error representations are sharper than the classical ones and lead to a more robust *p*-adaptive process. We also provide guidelines for finding operators delivering sharp error representation upper bounds. We further extend the results to a convection-dominated diffusion problem as well as to problems with discontinuous material coefficients. Finally, we consider a sonic logging-while-drilling problem to illustrate the applicability of the proposed method.

To increase the robustness of a Padé-based approximation of parametric solutions to finite element problems, an a priori estimate of the poles is proposed. The resulting original approach is shown to allow for a straightforward, efficient, subsequent Padé-based expansion of the solution vector components, overcoming some of the current convergence and robustness limitations. In particular, this enables for the intervals of approximation to be chosen a priori in direct connection with a given choice of Padé approximants. The choice of these approximants, as shown in the present work, is theoretically supported by the Montessus de Ballore theorem, concerning the convergence of a series of approximants with fixed denominator degrees. Key features and originality of the proposed approach are (1) a component-wise expansion which allows to specifically target subsets of the solution field and (2) the a priori, simultaneous choice of the Padé approximants and their associated interval of convergence for an effective and more robust approximation. An academic acoustic case study, a structural-acoustic application, and a larger acoustic problem are presented to demonstrate the potential of the approach proposed.

The finite analytic method (FAM) is developed to solve the 2D steady fluid flows in heterogeneous porous media with full tensor permeability on unstructured grids. The proposed FAM is constructed based upon the power-law analytic nodal solution in the angular domain with arbitrary shape. When approaching the grid node joining the subdomains, 3 different flow patterns may exist: power-law flow, linear flow, or the stagnant flow. Based on the nodal analytic solution, the triangle-based FAM is proposed. Numerical examples show that the proposed numerical scheme makes the convergences much quickly than the traditional methods, typically the weighted harmonic mean method under the cell refinement. In practical applications, the grid refinement parameter *n* = 2 or *n* = 3 is recommended, and the relative error of the calculated equivalent permeability will below 4% independent of the heterogeneity. In contrast, when using the traditional numerical scheme the refinement ratio for the grid cell needs to increase dramatically to get an accurate result, especially for strong heterogeneous porous medium.

Proper generalized decomposition (PGD) is often used for multiquery and fast-response simulations. It is a powerful tool alleviating the curse of dimensionality affecting multiparametric partial differential equations. Most implementations of PGD are intrusive extensions based on in-house developed FE solvers. In this work, we propose a nonintrusive PGD scheme using off-the-shelf FE codes (such as certified commercial software) as an external solver. The scheme is implemented and monitored by in-house flow-control codes. A typical implementation is provided with downloadable codes. Moreover, a novel parametric separation strategy for the PGD resolution is presented. The parametric space is split into two- or three-dimensional subspaces, to allow PGD technique solving problems with constrained parametric spaces, achieving higher convergence ratio. Numerical examples are provided. In particular, a practical example in biomechanics is included, with potential application to patient-specific simulation.

This paper presents an adaption of periodic boundary conditions (BC), which is termed tessellation BC. While periodic BC restrict strain localization zones to obey the periodicity of the microstructure, the proposed tessellation BC adjust the periodicity frame to meet the localization zone. Thereby, arbitrary developing localization zones are permitted. Still the formulation is intrinsically unsusceptible against spurious localization. Additionally, a modification of the Hough transformation is derived, which constitutes an unbiased criterion for the detection of the localization zone. The behavior of the derived formulation is demonstrated by various examples and compared with other BC. It is thereby shown that tessellation BC lead to a reasonable dependence of the effective stress on the localization direction. Furthermore, good convergence of stiffness values with increasing size of the representative volume element is shown as well as beneficial characteristics in use with strain softening material.

In this paper, a new numerical method, element differential method (EDM), is proposed for solving general thermal-mechanical problems. The key point of the method is the direct differentiation of the shape functions of Lagrange isoparametric elements used to characterize the geometry and physical variables. A set of analytical expressions for computing the first- and second-order partial derivatives of the shape functions with respect to global coordinates are derived. Based on these expressions, a new collocation method is proposed for establishing the system of equations, in which the equilibrium equations are collocated at nodes inside elements, and the traction equilibrium equations are collocated at interface nodes between elements and outer surface nodes of the problem. Attributed to the use of the Lagrange elements that can guarantee the variation of physical variables consistent through all elemental nodes, EDM has higher stability than the traditional collocation method. The other main features of EDM are that no mathematical or mechanical principles are required to set up the system of equations and no integrals are involved to form the coefficients of the system. A number of numerical examples of 2- and 3-dimensional problems are given to demonstrate the correctness and efficiency of the proposed method.

There is increasing interest in the material point method (MPM) as a means of modelling solid mechanics problems in which very large deformations occur, e.g. in the study of landslides and metal forming; however, some aspects vital to wider use of the method have to date been ignored, in particular methods for imposing essential boundary conditions in the case where the problem domain boundary does not coincide with the background grid element edges. In this paper, we develop a simple procedure originally devised for standard finite elements for the imposition of essential boundary conditions, for the MPM, expanding its capabilities to model boundaries of any inclination. To the authors' knowledge, this is the first time that a method has been proposed that allows arbitrary Dirichlet boundary conditions (zero and nonzero values at any inclination) to be imposed in the MPM. The method presented in this paper is different from other MPM boundary approximation approaches, in that (1) the boundaries are independent of the background mesh, (2) artificially stiff regions of material points are avoided, and (3) the method does not rely on *mirroring* of the problem domain to impose symmetry. The main contribution of this work is equally applicable to standard finite elements and the MPM.

An efficient method for generating the mass matrix inverse of structural dynamic problems is presented, which can be tailored to improve the accuracy of target frequency ranges and/or wave contents. The present method bypasses the use of biorthogonal construction of a kernel inverse mass matrix that requires special procedures for boundary conditions and free edges or surfaces and constructs the free-free inverse mass matrix using the standard FEM procedure. The various boundary conditions are realized by the the method of localized Lagrange multipliers. In particular, the present paper constructs the kernel inverse matrix by using the standard FEM elemental mass matrices. It is shown that the accuracy of the present inverse mass matrix is almost identical to that of a conventional consistent mass matrix or a combination of lumped and consistent mass matrices. Numerical experiments with the proposed inverse mass matrix are conducted to validate its effectiveness when applied to vibration analysis of bars, beams, and plain stress problems.

Gradient-dependent plasticity can be used to achieve mesh-objective results upon loss of well-posedness of the initial/boundary value problem because of the introduction of strain softening, non-associated flow, and geometric nonlinearity. A prominent class of gradient plasticity models considers a dependence of the yield strength on the Laplacian of the hardening parameter, usually an invariant of the plastic strain tensor. This inclusion causes the consistency condition to become a partial differential equation, in addition to the momentum balance. At the internal moving boundary, one has to impose appropriate boundary conditions on the hardening parameter or, equivalently, on the plastic multiplier. This internal boundary condition can be enforced without tracking the elastic-plastic boundary by requiring -continuity with respect to the plastic multiplier. In this contribution, this continuity has been achieved by using nonuniform rational B-splines as shape functions both for the plastic multiplier and for the displacements. One advantage of this isogeometric analysis approach is that the displacements can be interpolated one order higher, making it consistent with the interpolation of the plastic multiplier. This is different from previous approaches, which have been exploited. The regularising effect of gradient plasticity is shown for 1- and 2-dimensional boundary value problems.

This paper presents a level set–based topology optimization method for the design of an acoustic metamaterial with negative bulk modulus in an acoustic-elastic coupled system. To represent the coupled system, we propose a two-phase material model in which the solid and acoustic phases are mixed. This model allows the coupled system to be represented uniformly and avoids the need to impose the coupling conditions typically required in analyses of acoustic-elastic coupled systems. The effective bulk modulus of the acoustic metamaterial is represented using an S-parameter–based method and the optimization problem is formulated as a level set–based topology optimization. The topological derivatives for the two-phase material model are obtained using the relationship between the topological and shape derivatives. Numerical examples demonstrate the validity of the proposed two-phase material model and topological derivatives, and an optimal unit cell design for an acoustic metamaterial with negative bulk modulus is obtained at a targeted frequency.

A new method of topology optimization is introduced in which a continuous material field is combined with adaptive mesh refinement. Using a continuous material field with different analysis and design meshes allows the method to produce optimal designs that are free of numerical artifacts like checkerboard patterns and material islands. Adaptive mesh refinement is then applied to both meshes to precisely locate the optimal boundary of the final structure. A Helmholtz-type density filter is used to prevent the appearance of small topological features as the mesh refinement proceeds. Results are presented for several test problems, including problems with geometrically complex domain boundaries.

We explore diffuse formulations of Nitsche's method for consistently imposing Dirichlet boundary conditions on phase-field approximations of sharp domains. Leveraging the properties of the phase-field gradient, we derive the variational formulation of the diffuse Nitsche method by transferring all integrals associated with the Dirichlet boundary from a geometrically sharp surface format in the standard Nitsche method to a geometrically diffuse volumetric format. We also derive conditions for the stability of the discrete system and formulate a diffuse local eigenvalue problem, from which the stabilization parameter can be estimated automatically in each element. We advertise metastable phase-field solutions of the Allen-Cahn problem for transferring complex imaging data into diffuse geometric models. In particular, we discuss the use of mixed meshes, that is, an adaptively refined mesh for the phase-field in the diffuse boundary region and a uniform mesh for the representation of the physics-based solution fields. We illustrate accuracy and convergence properties of the diffuse Nitsche method and demonstrate its advantages over diffuse penalty-type methods. In the context of imaging-based analysis, we show that the diffuse Nitsche method achieves the same accuracy as the standard Nitsche method with sharp surfaces, if the inherent length scales, ie, the interface width of the phase-field, the voxel spacing, and the mesh size, are properly related. We demonstrate the flexibility of the new method by analyzing stresses in a human vertebral body.

This article presents a nonlinear solver combining regression analysis and a multiscale simulation scheme. First, the proposed method repeats microscopic analysis of a local simulation domain, which is extracted from the entire global domain, to statistically estimate the relation(s) between the value of a dependent variable at a point and values at surrounding points. The relation is called regression function. Subsequent global analysis reveals the behavior of the global domain with only coarse-grained points using the regression function quickly at low computational cost, which can be accomplished using a multiscale numerical solver, called the seamless-domain method. The objective of the study is to solve a nonlinear problem accurately and at low cost by combining the 2 techniques. We present an example problem of a nonlinear steady-state heat conduction analysis of a heterogeneous material. The proposed model using fewer than 1000 points generates a solution with precision similar to that of a standard finite-element solution using hundreds of thousands of nodes. To investigate the relationship between the accuracy and computational time, we apply the seamless-domain method under varying conditions such as the number of iterations of the prior analysis for statistical data learning.

We consider an efficient preconditioner for a boundary integral equation (BIE) formulation of the two-dimensional Stokes equations in porous media. While BIEs are well-suited for resolving the complex porous geometry, they lead to a dense linear system of equations that is computationally expensive to solve for large problems. This expense is further amplified when a significant number of iterations is required in an iterative Krylov solver such as generalized minimial residual method (GMRES). In this paper, we apply a fast inexact direct solver, the inverse fast multipole method, as an efficient preconditioner for GMRES. This solver is based on the framework of
-matrices and uses low-rank compressions to approximate certain matrix blocks. It has a tunable accuracy *ε* and a computational cost that scales as
. We discuss various numerical benchmarks that validate the accuracy and confirm the efficiency of the proposed method. We demonstrate with several types of boundary conditions that the preconditioner is capable of significantly accelerating the convergence of GMRES when compared to a simple block-diagonal preconditioner, especially for pipe flow problems involving many pores.

The Koiter-Newton method had recently demonstrated a superior performance for nonlinear analyses of structures, compared to traditional path-following strategies. The method follows a predictor-corrector scheme to trace the entire equilibrium path. During a predictor step, a reduced-order model is constructed based on Koiter's asymptotic postbuckling theory that is followed by a Newton iteration in the corrector phase to regain the equilibrium of forces. In this manuscript, we introduce a robust mixed solid-shell formulation to further enhance the efficiency of stability analyses in various aspects. We show that a Hellinger-Reissner variational formulation facilitates the reduced-order model construction omitting an expensive evaluation of the inherent fourth-order derivatives of the strain energy. We demonstrate that extremely large step sizes with a reasonable out-of-balance residual can be obtained with substantial impact on the total number of steps needed to trace the complete equilibrium path. More importantly, the numerical effort of the corrector phase involving a Newton iteration of the full-order model is drastically reduced thus revealing the true strength of the proposed formulation. We study a number of problems from engineering and compare the results to the conventional approach in order to highlight the gain in numerical efficiency for stability problems.

This paper presents a projection method to obtain high-resolution, manufacturable structures from efficient and coarse-scale homogenization-based topology optimization results. The presented approach bridges coarse and fine scale, such that the complex periodic microstructures can be represented by a smooth and continuous lattice on the fine mesh. A heuristic methodology allows control of the projected topology, such that a minimum length scale on both solid and void features is ensured in the final result. Numerical examples show excellent behavior of the method, where performances of the projected designs are almost equal to the homogenization-based solutions. A significant reduction in computational cost is observed compared to conventional topology optimization approaches.

Topology optimization using stress constraints and considering uncertainties is a serious challenge, since a reliability problem has to be solved for each stress constraint, for each element in the mesh. In this paper, an alternative way of solving this problem is used, where uncertainty quantification is performed through the first-order perturbation approach, with proper validation by Monte Carlo simulation. Uncertainties are considered in the loading magnitude and direction. The minimum volume problem subjected to local stress constraints is formulated as a robust problem, where the stress constraints are written as a weighted average between their expected value and standard deviation. The augmented Lagrangian method is used for handling the large set of local stress constraints, whereas a gradient-based algorithm is used for handling the bounding constraints. It is shown that even in the presence of small uncertainties in loading direction, different topologies are obtained when compared to a deterministic approach. The effect of correlation between uncertainties in loading magnitude and direction on optimal topologies is also studied, where the main observed result is loss of symmetry in optimal topologies.

In this work, a fully explicit partitioned method for the simulation of Fluid Structure Interaction (FSI) problems is presented. The fluid domain is modelled with an explicit Particle Finite Element Method (PFEM) based on the hypothesis of weak compressibility. The Lagrangian description of the fluid is particularly effective in the simulation of FSI problems with free surface flows and large structural displacements, since the fluid boundaries are automatically defined by the position of the mesh nodes. A distinctive feature of the proposed FSI strategy is that the solid domain is modelled using the explicit integration FEM in an off-the-shelf commercial software (Abaqus/Explicit). This allows to perform simulations with a complete and advanced description on the structural domain, including advanced structural material models and contact. The structure-to-fluid coupling algorithm is based on a technique derived from the Domain Decomposition Methods, namely, the Gravouil and Combescure algorithm. The method allows for arbitrarily large interface displacements using different time incrementation and nonconforming meshes in the different domains, which is an essential feature for the efficiency of an explicit solver involving different materials. The resulting fully explicit and fully lagrangian finite element approach is particularly appealing for the possibility of its efficient application in a large variety of highly non-linear engineering problems.

The conventional approach to construct quadratic elements for an *n*-sided polygon will yield *n*(*n*+1)/2 shape functions, which increases the computational effort. It is well known that the serendipity elements based on isoparametric formulation suffers from mesh distortion. Floater and Lai proposed a systematic way to construct higher-order approximations over arbitrary polygons using the generalized barycentric and triangular coordinates. This approach ensures 2*n* shape functions with nodes only on the boundary of the polygon. In this paper, we extend the polygonal splines approach to 3 dimensions and construct serendipity shape functions over hexahedra and convex polyhedra. This is done by expressing the shape functions using the barycentric coordinates and the local tetrahedral coordinates. The quadratic shape functions possess Kronecker delta property and satisfy constant, linear, and quadratic precision. The accuracy and the convergence properties of the quadratic serendipity shape elements are demonstrated with a series of standard patch tests. The numerical results show that the quadratic serendipity elements pass the patch test, yield optimal convergence rate, and can tolerate *extreme* mesh distortion.

The most common technique for the numerical implementation of peridynamic theory is the uniform discretization together with constant horizon. However, unlike the nonuniform discretization and varying horizons, it is not a natural and intrinsic component of the adaptive refinement analysis and multiscale modeling. Besides, it encounters discretization difficulty in analyzing irregular structures. Therefore, to analyze problems with nonuniform discretization and varying horizons and solve the resulting problems of ghost forces and spurious wave reflection, the dual-horizon peridynamics based on uniform discretization is extended to the nonuniform discretization based on Voronoi diagrams, for which we call it Voronoi-based peridynamics. We redefine the damage definition as well. Next, an adaptive refinement analysis method based on the proposed Voronoi-based peridynamics and its numerical implementation is introduced. Finally, the traditional bond-based peridynamics and the Voronoi-based peridynamics with or without refinement are used to simulate 4 benchmark problems. The examples of 2-D quasi-static elastic deformation, 2-D wave propagation, 2-D dynamic crack growth, and 3-D simulation of the Kalthoff-Winkler experiment demonstrate the efficiency and effectivity of the proposed Voronoi-based peridynamics. Further, the adaptive refinement analysis can be used to obtain reasonable crack path and crack propagation speed with reduced computational burden.

The material point method for the analysis of deformable bodies is revisited and originally upgraded to simulate crack propagation in brittle media. In this setting, phase-field modelling is introduced to resolve the crack path geometry. Following a particle in cell approach, the coupled continuum/phase-field governing equations are defined at a set of material points and interpolated at the nodal points of an Eulerian, ie, non-evolving, mesh. The accuracy of the simulated crack path is thus decoupled from the quality of the underlying finite element mesh and relieved from corresponding mesh-distortion errors. A staggered incremental procedure is implemented for the solution of the discrete coupled governing equations of the phase-field brittle fracture problem. The proposed method is verified through a series of benchmark tests while comparisons are made between the proposed scheme, the corresponding finite element implementation, and experimental results.

Unsaturated soils are solid-water-air systems that include a solid skeleton, pore water, and pore air. Heterogeneities in porosity or degree of saturation are salient features of unsaturated soils. These heterogeneities may trigger localized deformation (eg, shear banding) in such materials as demonstrated by numerical simulations via a pseudo three-phase model. In this article, we formulate a true three-phase mathematical framework implemented via stabilized low-order mixed finite elements. With this mathematical framework, we study the evolution of pore air pressure and its role in the inception of strain localization triggered by initial heterogeneity either in porosity or suction. The numerical simulations show that pore air pressure is nonzero and nonuniform in the process of progressive failure in unsaturated soils. The heterogeneity of pore air pressure may also play a significant role in the onset of localized deformation of unsaturated soils. Therefore, a three-phase model considering the pore air phase is physically more appropriate for modeling strain localization in unsaturated soils.

Topology optimization is a methodology for assigning material or void to each point in a design domain in a way that extremizes some objective function, such as the compliance of a structure under given loads, subject to various imposed constraints, such as an upper bound on the mass of the structure. Geometry projection is a means to parameterize the topology optimization problem, by describing the design in a way that is independent of the mesh used for analysis of the design's performance; it results in many fewer design parameters, necessarily resolves the ill-posed nature of the topology optimization problem, and provides sharp descriptions of the material interfaces. We extend previous geometric projection work to 3 dimensions and design unit cells for lattice materials using inverse homogenization. We perform a sensitivity analysis of the geometric projection and show it has smooth derivatives, making it suitable for use with gradient-based optimization algorithms. The technique is demonstrated by designing unit cells comprised of a single constituent material plus void space to obtain light, stiff materials with cubic and isotropic material symmetry. We also design a single-constituent isotropic material with negative Poisson's ratio and a light, stiff material comprised of 2 constituent solids plus void space.

This paper presents a method for the optimization of multicomponent structures comprised of 2 and 3 materials considering large motion sliding contact and separation along interfaces. The structural geometry is defined by an explicit level set method, which allows for both shape and topology changes. The mechanical model assumes finite strains, a nonlinear elastic material behavior, and a quasi-static response. Identification of overlapping surface position is handled by a coupled parametric representation of contact surfaces. A stabilized Lagrange method and an active set strategy are used to model frictionless contact and separation. The mechanical model is discretized by the extended FEM, which maintains a clear definition of geometry. Face-oriented ghost penalization and dynamic relaxation are implemented to improve the stability of the physical response prediction. A nonlinear programming scheme is used to solve the optimization problem, which is regularized by introducing a perimeter penalty into the objective function. Design sensitivities are determined by the adjoint method. The main characteristics of the proposed method are studied by numerical examples in 2 dimensions. The numerical results demonstrate improved design performance when compared to models optimized with a small strain assumption. Additionally, examples with load path dependent objectives display nonintuitive designs.

The material point method (MPM) enhanced with B-spline basis functions, referred to as B-spline MPM (BSMPM), is developed and demonstrated using representative quasi-static and dynamic example problems. Smooth B-spline basis functions could significantly reduce the cell-crossing error as known for the original MPM. A Gauss quadrature scheme is designed and shown to be able to diminish the quadrature error in the BSMPM analysis of large-deformation problems for the improved accuracy and convergence, especially with the quadratic B-splines. Moreover, the increase in the order of the B-spline basis function is also found to be an effective way to reduce the quadrature error and to improve accuracy and convergence. For plate impact examples, it is demonstrated that the BSMPM outperforms the generalized interpolation material point (GIMP) and convected particle domain interpolation (CPDI) methods in term of the accuracy of representing stress waves. Thus, the BSMPM could become a promising alternative to the MPM, GIMP, and CPDI in solving certain types of transient problems.

]]>An arbitrary Lagrangian-Eulerian framework, which combines the advantages of both Lagrangian and Eulerian methods, is presented to solve incompressible multiphase flow problems. The incompressible Navier-Stokes equations are discretized using the side-centered unstructured finite volume method, where the velocity vector components are defined at the midpoint of each cell face, while the pressure term is defined at element centroids. The pressure field is treated to be discontinuous across the interface with the discontinuous treatment of density and viscosity. The surface tension term at the interface is treated as a force tangent to the interface and computed with several different approaches including the use of Legendre polynomials. In addition, the several different discretizations of interface kinematic boundary conditions are investigated. For the application of the interface kinematic boundary condition, a special attention is given to satisfy both local and global discrete geometric conservation law to conserve the total mass of both species at machine precision. The mesh vertices are deformed by solving the linear elasticity equations due to the normal displacement of interface. The resulting algebraic equations are solved in a fully coupled manner, and a one-level restricted additive Schwarz preconditioner with a block-incomplete factorization within each partitioned subdomain is used for the resulting fully coupled system. The method is validated by simulating the classical benchmark problem of a single rising bubble in a viscous fluid due to buoyancy. The results of numerical simulations are found out to be in an excellent agreement with the earlier results in the literature. The mass of the bubble is conserved, and discontinuous pressure field is obtained to avoid errors due to the incompressibility condition in the vicinity of the interface, where the density and viscosity jumps occur.

Adaptive hierarchical refinement in isogeometric analysis is developed to model cohesive crack propagation along a prescribed interface. In the analysis, the crack is introduced by knot insertion in the NURBS basis, which yields continuous basis functions. To capture the stress state smoothly ahead of the crack tip, the hierarchical refinement of the spline basis functions is used starting from a coarse initial mesh. A multilevel mesh is constructed, with a fine mesh used for quantifying the stresses ahead of the crack tip, knot insertion to insert the crack, and coarsening in the wake of the crack tip, since a lower resolution suffices there. This technique can be interpreted as a moving mesh around the crack tip. To ensure compatibility with existing finite element programs, an element-wise point of view is adopted using Bézier extraction. A detailed description is given how the approach can be implemented in a finite element data structure. The accuracy of the approach to cohesive fracture modelling is demonstrated by several numerical examples, including a double cantilever beam, an L-shaped specimen, and a fibre embedded in an epoxy matrix.

In this article, we develop a dynamic version of the variational multiscale (D-VMS) stabilization for nearly/fully incompressible solid dynamics simulations of viscoelastic materials. The constitutive models considered here are based on Prony series expansions, which are rather common in the practice of finite element simulations, especially in industrial/commercial applications. Our method is based on a mixed formulation, in which the momentum equation is complemented by a pressure equation in rate form. The unknown pressure, displacement, and velocity are approximated with piecewise linear, continuous finite element functions. To prevent spurious oscillations, the pressure equation is augmented with a stabilization operator specifically designed for viscoelastic problems, in that it depends on the viscoelastic dissipation. We demonstrate the robustness, stability, and accuracy properties of the proposed method with extensive numerical tests in the case of linear and finite deformations.

A machine learning–based framework for modeling the error introduced by surrogate models of parameterized dynamical systems is proposed. The framework entails the use of high-dimensional regression techniques (eg, random forests, and LASSO) to map a large set of inexpensively computed “error indicators” (ie, features) produced by the surrogate model at a given time instance to a prediction of the surrogate-model error in a quantity of interest (QoI). This eliminates the need for the user to hand-select a small number of informative features. The methodology requires a training set of parameter instances at which the time-dependent surrogate-model error is computed by simulating both the high-fidelity and surrogate models. Using these training data, the method first determines regression-model locality (via classification or clustering) and subsequently constructs a “local” regression model to predict the time-instantaneous error within each identified region of feature space. We consider 2 uses for the resulting error model: (1) as a correction to the surrogate-model QoI prediction at each time instance and (2) as a way to statistically model arbitrary functions of the time-dependent surrogate-model error (eg, time-integrated errors). We apply the proposed framework to model errors in reduced-order models of nonlinear oil-water subsurface flow simulations, with time-varying well-control (bottom-hole pressure) parameters. The reduced-order models used in this work entail application of trajectory piecewise linearization in conjunction with proper orthogonal decomposition. When the first use of the method is considered, numerical experiments demonstrate consistent improvement in accuracy in the time-instantaneous QoI prediction relative to the original surrogate model, across a large number of test cases. When the second use is considered, results show that the proposed method provides accurate statistical predictions of the time- and well-averaged errors.

We propose a new computational framework for the treatment of acousto-magneto-mechanical coupling that arises in low-frequency electro-magneto-mechanical systems such as magnetic resonance imaging scanners. Our transient Newton–Raphson strategy involves the solution of a monolithic system obtained from the linearisation of the coupled system of equations. Moreover, this framework, in the case of excitation from static and harmonic current sources, allows us to propose a simple linearised system and rigorously motivate a single-step strategy for understanding the response of systems under different frequencies of excitation. Motivated by the need to solve industrial problems rapidly, we restrict ourselves to solving problems consisting of axisymmetric geometries and current sources. Our treatment also discusses in detail the computational requirements for the solution of these coupled problems on unbounded domains and the accurate discretisation of the fields using *h**p*–finite elements. We include a set of academic and industrially relevant examples to benchmark and illustrate our approach. Copyright © 2017 The Authors. *International Journal for Numerical Methods in Engineering* Published by John Wiley & Sons, Ltd.

Numerical stability by using certain time integration scheme is a critical issue for accurate simulation of discontinuous deformations of solids. To investigate the effects of the time integration schemes on the numerical stability of the numerical manifold method, the implicit time integration schemes, ie, the Newmark, the HHT-α, and the WBZ-α methods, and the explicit time integration algorithms, ie, the central difference, the Zhai's, and Chung-Lee methods, are implemented. Their performance is examined by conducting transient response analysis of an elastic strip subjected to constant loading, impact analysis of an elastic rod with an initial velocity, and excavation analysis of jointed rock masses, respectively. Parametric studies using different time steps are conducted for different time integration algorithms, and the convergence efficiency of the open-close iterations for the contact problems is also investigated. It is proved that the Hilber-Hughes-Taylor-α (HHT-α), Wood-Bossak-Zienkiewicz-α (WBZ-α), Zhai's, and Chung-Lee methods are more attractive in solving discontinuous deformation problems involving nonlinear contacts. It is also found that the examined explicit algorithms showed higher computational efficiency compared to those implicit algorithms within acceptable computational accuracy.

We present and analyze a validation procedure for a given state estimate *u*^{⋆} of the true field *u*^{true} based on Monte Carlo sampling of experimental observation functionals. Our method provides, given a set of *N* possibly noisy local experimental observation functionals over the spatial domain Ω, confidence intervals for the *L*^{2}(Ω) error in state and the error in *L*^{2}(Ω) outputs. For *L*^{2}(Ω) outputs, our approach also provides a confidence interval for the output itself, which can be used to improve the initial output estimate based on *u*^{⋆}. Our approach implicitly takes advantage of variance reduction, through the proximity of *u*^{⋆} to *u*^{true}, to provide tight confidence intervals even for modest values of *N*. We present results for a synthetic model problem to illustrate the elements of the methodology and confirm the numerical properties suggested by the theory. Finally, we consider an experimental thermal patch configuration to demonstrate the applicability of our approach to real physical systems.

The present study proposes a method of micro-macro concurrent topology optimization for a two-phase nonlinear solid to minimize the end compliance of its macrostructure undergoing large deformation. To reduce the computational costs to solve a 2-scale boundary value problem under geometrically nonlinear setting, we use the so-called method of *decoupling multiscale structural analysis*, in which the microscopic and macroscopic boundary value problems are decoupled in the homogenization process. An isotropic hyperelasticity model is used for the constitutive model for microstructures, while an orthotropic one is assumed to represent the macroscopic material behavior. Owing to this decoupling framework, the micro-macro concurrent optimization problem can be split into 2 individual problems at the microscale and macroscale for the sake of algorithmic simplicity. Also, a 2-scale adjoint sensitivity analysis can be performed within the framework of computational homogenization. It is verified from a series numerical examples that the proposed method is capable of computing the optimal structures at both microscale and macroscale, according to the level of applied load.

In this paper, new spherical Hankel shape functions are used to reformulate boundary element method for 2-dimensional elastostatic and elastodynamic problems. To this end, the dual reciprocity boundary element method is reconsidered by using new spherical Hankel shape functions to approximate the state variables (displacements and tractions) of Navier's differential equation. Using enrichment of a class of radial basis functions (RBFs), called spherical Hankel RBFs hereafter, the interpolation functions of a Hankel boundary element framework has been derived. For this purpose, polynomial terms are added to the functional expansion that only uses spherical Hankel RBF in the approximation. In addition to polynomial function fields, the participation of spherical Bessel function fields has also increased robustness and efficiency in the interpolation. It is very interesting that there is no Runge phenomenon in equispaced Hankel macroelements, unlike equispaced classic Lagrange ones. Several numerical examples are provided to demonstrate the effectiveness, robustness and accuracy of the proposed Hankel shape functions and in comparison with the classic Lagrange ones, they show much more accurate and stable results.

To coarsen a mesh, we usually remove a set of selected nodes one by one. Currently, the basic operation used to remove a node is edge collapsing, which does not perform well when applied to handling narrow regions in a tetrahedron mesh and could produce low-quality elements or even fail to give valid results. To overcome the drawbacks of edge collapsing, we present a new node-removal operator created by revising a topological transformation called small polyhedron reconnection. This new operator can guarantee success if the cavity that forms after a node is removed is meshable, and it produces higher-quality results and keeps the nodes unmoved, which is preferred for applications such as multigrid hierarchies.

In addition, 2 other aspects of mesh coarsening that determine whether a node should be removed and the sequence in which to remove the selected nodes are also studied. Our strategy consists of constructing a coarse node set using the sphere-packing method and removing the nodes in a reversed kd-tree sequence. The excellent performance of the new method is demonstrated by applying it to examples of adaptive meshing and multigrid hierarchy creation and comparing the results with those of the edge collapsing method.

The paper introduces a novel multiresolution scheme to topology optimization in the framework of the isogeometric analysis. A new variable parameter space is added to implement multiresolution topology optimization based on the Solid Isotropic Material with Penalization approach. Design density variables defined in the variable space are used to approximate the element analysis density by the bivariate B-spline basis functions, which are easily obtained using *k*-refinement strategy in the isogeometric analysis. While the nonuniform rational B-spline basis functions are used to exactly describe geometric domains and approximate unknown solutions in finite element analysis. By applying a refined sensitivity filter, optimized designs include highly discrete solutions in terms of solid and void materials without using any black and white projection filters. The Method of Moving Asymptotes is used to solve the optimization problem. Various benchmark test problems including plane stress, compliant mechanism inverter, and 2-dimensional heat conduction are examined to demonstrate the effectiveness and robustness of the present method.

No abstract is available for this article.

]]>We propose mixed hybrid Discontinuous Galerkin Finite Element (DGFEM) formulations for the Stokes problem characterized by the introduction of Lagrange multipliers associated with the traces of the velocity and pressure fields on the edges of the elements to weakly impose the transmission conditions. Both velocity and pressure multipliers are stabilized and, as a consequence of these stabilizations, we prove existence and uniqueness of solution for the local problems. All velocity and pressure degrees-of-freedom can be eliminated at the element level by static condensation leading to a global problem in the multipliers only. The proposed methodology is able to recover stability of very convenient choices of finite element spaces, such as those adopting equal order polynomial approximations for all fields. Numerical experiments illustrate the flexibility and robustness of the proposed formulations and show optimal rates of convergence. Copyright © 2017 John Wiley & Sons, Ltd.

This paper presents a comprehensive study on the use of Irwin's crack closure integral for direct evaluation of mixed-mode stress intensity factors (SIFs) in curved crack problems, within the extended finite element method. The approach employs high-order enrichment functions derived from the standard Williams asymptotic solution, and SIFs are computed in closed form without any special post-processing requirements. Linear triangular elements are used to discretize the domain, and the crack curvature within an element is represented explicitly. An improved quadrature scheme using high-order isoparametric mapping together with a generalized Duffy transformation is proposed to integrate singular fields in tip elements with curved cracks. Furthermore, because the Williams asymptotic solution is derived for straight cracks, an appropriate definition of the angle in the enrichment functions is presented and discussed. This contribution is an important extension of our previous work on straight cracks and illustrates the applicability of the SIF extraction method to curved cracks.

The performance of the method is studied on several circular and parabolic arc crack benchmark examples. With two layers of elements enriched in the vicinity of the crack tip, striking accuracy, even on relatively coarse meshes, is obtained, and the method converges to the reference SIFs for the circular arc crack problem with mesh refinement. Furthermore, while the popular interaction integral (a variant of the J-integral method) requires special auxiliary fields for curved cracks and also needs cracks to be sufficiently apart from each other in multicracks systems, the proposed approach shows none of those limitations. Copyright © 2017 John Wiley & Sons, Ltd.

Over the past two decades, meshfree methods have undergone significant development as a numerical tool to solve partial differential equations (PDEs). In contrast to finite elements, the basis functions in meshfree methods are smooth (nonpolynomial functions), and they do not rely on an underlying mesh structure for their construction. These features render meshfree methods to be particularly appealing for higher-order PDEs and for large deformation simulations of solid continua. However, a deficiency that still persists in meshfree Galerkin methods is the inaccuracies in numerical integration, which affects the consistency and stability of the method. Several previous contributions have tackled the issue of integration errors with an eye on consistency, but without explicitly ensuring stability. In this paper, we draw on the recently proposed virtual element method, to present a formulation that guarantees both the consistency and stability of the approximate bilinear form. We adopt maximum-entropy meshfree basis functions, but other meshfree basis functions can also be used within this framework. Numerical results for several two-dimensional and three-dimensional elliptic (Poisson and linear elastostatic) boundary-value problems that demonstrate the effectiveness of the proposed formulation are presented. Copyright © 2017 John Wiley & Sons, Ltd.

The boundary condition represented by polygons in the moving particle semi-implicit method can accurately represent geometries and treat complex geometry with high efficiency. However, inaccurate wall contribution to the Poisson's equation leads to drastic numerical oscillation. To address this issue, in this research, we analyzed the problems of the Poisson's equation used in the boundary condition represented by polygons. The new Poisson's equation is proposed based on the improved source term (Tanaka and Masunaga, Trans Jpn Soc Comput Eng Sci, 2008). The asymmetric gradient model (Khayyer and Gotoh, Coastal Engineering Journal, 2008) is also adopted to further suppress the numerical oscillation of fluid particles. The proposed method can dramatically improve the pressure distribution to arbitrary geometry in three dimensions and keep the efficiency. Four examples including the hydrostatic simulation, dam break simulation, and two complex geometries are verified to show the general applicability of the proposed method. Copyright © 2017 John Wiley & Sons, Ltd.

This work gives new statement of the vertex solution theorem for exact bounds of the solution to linear interval equations and its novel proof by virtue of the convex set theory. The core idea of the theorem is to transform linear interval equations into a series of equivalent deterministic linear equations. Then, the important theorem is extended to find the upper and lower bounds of static displacements of structures with interval parameters. Following discussions about the computational efforts, a coupled framework based on vertex method (VM) is established, which allows us to solve many large-scale engineering problems with uncertainties using deterministic finite element software. Compared with the previous works, the contribution of this work is not only to obtain the exact bounds of static displacements but also lay the foundation for development of an easy-to-use interval finite element software. Numerical examples demonstrate the good accuracy of VM. Meanwhile, the implementation of VM and availability of the coupled framework are demonstrated by engineering example. Copyright © 2017 John Wiley & Sons, Ltd.

A novel density-based topology optimization framework for plastic energy absorbing structural designs with maximum damage constraint is proposed. This framework enables topologies to absorb large amount of energy via plastic work before failure occurs. To account for the plasticity and damage during the energy absorption, a coupled elastoplastic ductile damage model is incorporated with topology optimization. Appropriate material interpolation schemes are proposed to relax the damage in the low-density regions while still ensuring the convergence of Newton-Raphson solution process in the nonlinear finite element analyses. An effective method for obtaining path-dependent sensitivities of the plastic work and maximum damage via adjoint method is presented, and the sensitivities are verified by the central difference method. The effectiveness of the proposed methodology is demonstrated through a series of numerical examples. Copyright © 2017 John Wiley & Sons, Ltd.

We present the theory of novel time-stepping algorithms for general nonlinear, non-smooth, coupled, and thermomechanical problems. The proposed methods are thermodynamically consistent in the sense that their solutions rigorously comply with the two laws of thermodynamics: for isolated systems, they preserve the total energy and the entropy never decreases. Extending previous works on the subject, the newly proposed integrators are applicable to coupled mechanical systems with non-smooth kinetics and can be formulated in arbitrary variables. The ideas are illustrated with a simple non-smooth problem: a rheological model for a thermo-elasto-plastic material with hardening. Numerical simulations verify the qualitative features of the proposed methods and illustrate their excellent numerical stability, which stems precisely from their ability to preserve the structure of the evolution equations they discretize. Copyright © 2017 John Wiley & Sons, Ltd.

The numerical manifold method (NMM) builds up a unified framework that is used to describe continuous and discontinuous problems; it is an attractive method for simulating a cracking phenomenon. Taking into account the differences between the generalized degrees of freedom of the physical patch and nodal displacement of the element in the NMM, a decomposition technique of generalized degrees of freedom is deduced for mixed mode crack problems. An analytic expression of the energy release rate, which is caused by a virtual crack extension technique, is proposed. The necessity of using a symmetric mesh is demonstrated in detail by analysing an additional error that had previously been overlooked. Because of this necessity, the local mathematical cover refinement is further applied. Finally, four comparison tests are given to illustrate the validity and practicality of the proposed method. The aforementioned aspects are all implemented in the high-order NMM, so this study can be regarded as the development of the virtual crack extension technique and can also be seen as a prelude to an h-version high-order NMM. Copyright © 2017 John Wiley & Sons, Ltd.

The scaled boundary radial point interpolation method (SBRPIM), a new semi-analytical technique, is introduced and applied to the analysis of the stress–strain problems. The proposed method combines the advantages of the scaled boundary finite element method and the boundary radial point interpolation method. In this method, no mesh is required, nodes are scattered only on the boundary of the domain, no fundamental solution is required, and as the shape functions constructed based on the radial point interpolation method possess the Kronecker delta function property, the boundary conditions of problems can be imposed accurately without additional efforts. The main ideas of the SBRPIM are introducing a new method based on boundary scattered nodes without the need to element connectivity information, satisfying Kronecker delta function property, and being capable to handle singular problems. The equations of the SBRPIM in stress–strain fields are outlined in this paper. Several benchmark examples of 2-D elastostatic are analyzed to validate the accuracy and efficiency of the proposed method. It is found that the SBRPIM is very easy to implement and the obtained results of the proposed method show a very good agreement with the analytical solution. Copyright © 2017 John Wiley & Sons, Ltd.