Search

Partial differential equations having diffusive, convective, and reactive terms appear in the modeling of a large variety of processes in several branches of science. Often, several species or components interact with each other, rendering strongly coupled systems of convection-diffusion-reaction equations. Exact solutions are available in extremely few cases lacking practical interest due to the simplifications made to render such equations amenable by analytical tools. Then, numerical approximation remains the best strategy for solving these problems. The properties of these systems of equations, particularly the lack of sufficient physical diffusion, cause most traditional numerical methods to fail, with the appearance of violent and nonphysical oscillations, even for the single equation case. For systems of equations, the situation is even harder due to the lack of fundamental principles guiding numerical discretization. Therefore, strategies must be developed in order to obtain physically meaningful and numerically stable approximations. Such stabilization techniques have been extensively developed for the single equation case in contrast to the multiple equations case. This paper presents a perturbation-based stabilization technique for coupled systems of one-dimensional convection-diffusion-reaction equations. Its characteristics are discussed, providing evidence of its versatility and effectiveness through a thorough assessment.

SUMMARY 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 which 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 two- and three-dimensional problems are given to demonstrate the correctness and efficiency of the proposed method.

We propose a formulation for linear elastic fracture mechanics (LEFM) in which the stress intensity factors (SIF) are found directly from the solution vector of an extended boundary element method (XBEM) formulation. The enrichment is embedded in the BEM formulation, rather than adding new degrees of freedom for each enriched node. Therefore, a very limited number of new degrees of freedom is added to the problem, which contributes to preserving the conditioning of the linear system of equations. The Stroh formalism is used to provide BEM fundamental solutions for any degree of anisotropy, and these are used for both conventional and enriched degrees of freedom. Several numerical examples are shown with benchmark solutions to validate the proposed method.

This manuscript presents the formulation and application of the Green’s Discrete Transformation Method (GDTM) for the meshfree simulation of transient diffusion problems, including those with moving boundaries. The GDTM implements a linear combination of time-dependent Green’s basis functions defined on a set of source points to approximate the field in the form of a solution series. A discrete transformation is implemented to evaluate unknown coefficients of this series, which eliminates the need to use time integration schemes. We will study the optimal number and location of the GDTM source points that yield the highest level of accuracy, while maintaining a manageable condition number for the resulting linear system of equations. The optimal values of these parameters, which are inherently independent of the domain geometry, are determined such that the basis functions have appropriate features for approximating the field. A comprehensive convergence study is presented to show the precision and convergence rate of the GDTM for modeling various diffusion problems. We also demonstrate the application of this method for simulating three diffusion problems with complex and evolving morphologies: heat transfer in a turbine blade, thermal response of a porous material, and localized (pitting) corrosion in stainless steel. Copyright c 2015 John Wiley & Sons, Ltd.

Summary We propose a coupled boundary element method (BEM) and a finite element method (FEM) for modelling localized damage growth in structures. BEM offers the flexibility of modelling large domains efficiently, while the non-linear damage growth is accurately accounted by a local FEM mesh. An integral-type nonlocal continuum damage mechanics with adapting FEM mesh is used to model multiple damage zones and follow their propagation in the structure. Strong form coupling, BEM hosted, is achieved using Lagrange multipliers. Because the non-linearity is isolated in the FEM part of the system of equations, the system size is reduced using Schur complement approach, then the solution is obtained by a monolithic Newton method that is used to solve both domains simultaneously. The coupled BEM/FEM approach is verified by a set of convergence studies, where the reference solution is obtained by a fine FEM. In addition, the method is applied to multiple fractures growth benchmark problems and shows good agreement with the literature. Copyright © 2015 John Wiley & Sons, Ltd.

A methodology is proposed in this paper to construct an adaptive sparse polynomial chaos (PC) expansion of the response of stochastic systems whose input parameters are independent random variables modeled as random fields. The proposed methodology utilizes the concept of variability response function (VRF) in order to compute an a priori low cost estimate of the spatial distribution of the second-order error of the response, as a function of the number of terms used in the truncated Karhunen-Lo eve (KL) expansion. This way the influence of the response variance to the spectral content (correlation structure) of the random input is taken into account through a spatial variation of the truncated KL terms. The criterion for selecting the number of KL terms at different parts of the structure is the uniformity of the spatial distribution of the second order error. This way a significantly reduced number of PC coefficients, with respect to classical PC expansion, is required in order to reach a uniformly distributed target second-order error. This results in an increase of sparsity of the coefficient matrix of the corresponding linear system of equations leading to an enhancement of the computational efficiency of the spectral stochastic finite element method (SSFEM). Copyright c © 2014 John Wiley & Sons, Ltd.

Summary An algorithm is derived for the computation of eigenpair derivatives of asymmetric quadratic eigenvalue problem with distinct and repeated eigenvalues. In the proposed method, the eigenvector derivatives of the damped systems are divided into a particular solution and a homogeneous solution. By introducing an additional normalization condition, we construct two extended systems of linear equations with nonsingular coefficient matrices to calculate the particular solution. The method is numerically stable, and the homogeneous solutions are computed by the second-order derivatives of the eigenequations. Two numerical examples are used to illustrate the validity of the proposed method. Copyright © 2015 John Wiley & Sons, Ltd.

SUMMARY In this work, an explicit family of time marching procedures with adaptive dissipation control is introduced. The proposed technique is conditionally stable, second-order accurate, and has controllable algorithm dissipation, which adapts according to the properties of the governing system of equations. Thus, spurious modes can be more effectively dissipated and accuracy is improved. Because this is an explicit time integration technique, the new family is quite efficient, requiring no system of equations to be dealt with at each time step. Moreover, the technique is simple and very easy to implement. Numerical results are presented along the paper, illustrating the good performance of the proposed method, as well as its potentialities. Copyright © 2014 John Wiley & Sons, Ltd.

In the present work, a monolithic solution approach for thermo-structure interaction problems motivated by the challenging application of the behaviour of rocket nozzles is proposed. Structural and thermal fields are independently discretised via finite elements. The resulting system of equations is solved via a monolithic thermo-structure interaction scheme, which is constructed by a block Gauss-Seidel preconditioner in combination with algebraic multigrid methods. The proposed method is tested for three numerical examples, the second Danilovskaya problem, a simplified rocket nozzle configuration and an internal loaded hollow sphere. Good agreement of the numerical results with results from the literature is observed. Furthermore, it is shown that the monolithic solution algorithm can handle the complete range of the parameter spectrum, whereas partitioned algorithms are limited to a certain parameter range only. Moreover, the monolithic algorithm exhibits improved efficiency and robustness compared to partitioned algorithms.

SUMMARY The present work deals with a new approach to frictional large deformation contact problems. In particular, a new formulation of the frictional kinematics is introduced that is based on a specific augmentation technique used for the introduction of additional variables. This augmentation technique substantially simplifies the formulation of the whole system. A size reduction of the resulting system of algebraic equations is proposed. Consequently, the augmentation technique does not lead to an increase in size of the algebraic system of equations to be ultimately solved. The size reduction retains the simplicity of the formulation and preserves important conservation laws such as conservation of angular momentum. Copyright © 2013 John Wiley & Sons, Ltd.

SUMMARY A unified framework of dual-primal finite element tearing and interconnecting (FETI-DP) algorithms is proposed for solving the system of linear equations arising from the mixed finite element approximation of incompressible Stokes equations. A distinctive feature of this framework is that it allows using both continuous and discontinuous pressures in the algorithm, whereas previous FETI-DP methods only apply to discontinuous pressures. A preconditioned conjugate gradient method is used in the algorithm with either a lumped or a Dirichlet preconditioner, and scalable convergence rates are proved. This framework is also used to describe several previously developed FETI-DP algorithms and greatly simplifies their analysis. Numerical experiments of solving a two-dimensional incompressible Stokes problem demonstrate the performances of the discussed FETI-DP algorithms represented under the same framework.Copyright © 2012 John Wiley & Sons, Ltd.

An enriched partition of unity FEM is developed to solve time-dependent diffusion problems. In the present formulation, multiple exponential functions describing the spatial and temporal diffusion decay are embedded in the finite element approximation space. The resulting enrichment is in the form of a local asymptotic expansion. Unlike previous works in this area where the enrichment must be updated at each time step, here, the temporal decay in the solution is embedded in the asymptotic expansion. Thus, the system matrix that is evaluated at the first time step may be decomposed and retained for every next time step by just updating the right-hand side of the linear system of equations. The advantage is a significant saving in the computational effort where, previously, the linear system must be reevaluated and resolved at every time step. In comparison with the traditional finite element analysis with p-version refinements, the present approach is much simpler, more efficient, and yields more accurate solutions for a prescribed number of DoFs. Numerical results are presented for a transient diffusion equation with known analytical solution. The performance of the method is analyzed on two applications: the transient heat equation with a single source and the transient heat equation with multiple sources. The aim of such a method compared with the classical FEM is to solve time-dependent diffusion applications efficiently and with an appropriate level of accuracy.

SUMMARY This paper analyses the buckling behaviour of single-walled and double-walled carbon nanotubes. The total potential of the atomic structure consists of the bonded energy and the non-bonded energy, both resulting from interatomic potentials, as well as the energy of external contributions. In particular, the influence of the in-layer and inter-layer non-bonded interactions is investigated. These non-bonded interactions are important to avoid the nanotubes from self-intersection or penetration and govern the morphology of the buckled tubes. The simulation model is based on a molecular statics approach embedded in the finite element framework. The search for equilibrium configurations of the structure leads to a non-linear system of equations, which is linearised and solved iteratively. Therefore, the model relies on the fully non-linear description of the interatomic potential and the atomic kinematics. Stability points of the loaded carbon nanotubes are detected by an accompanying eigenvalue analysis in combination with a bisection algorithm. By means of branch switching, the continuation of the non-linear load-deformation path in the postbuckling regime is enabled. With this framework, the buckling characteristics of carbon nanotubes in consequence of different loading conditions is studied. Results of numerical simulations are given for carbon nanotubes under torsion, axial compression and bending. Copyright © 2012 John Wiley & Sons, Ltd.

SUMMARY In this paper, an exact reanalysis algorithm based on an incremental Cholesky factorization is presented, which can solve a linear system of equations when a small portion of the coefficient matrix is modified. The efficiency of the proposed algorithm is demonstrated by modeling quasi-static crack growth in the extended finite element method. The example presented shows that a 60% to 70% reduction in computational time is achievable by using the reanalysis approach for solving crack growth problems. It is shown that the reanalysis approach has increasing benefits as the mesh density increases. Copyright © 2012 John Wiley & Sons, Ltd.

SUMMARY Frequency sweeps in structural dynamics, acoustics, and vibro-acoustics require evaluating frequency response functions for a large number of frequencies. The brute force approach for performing these sweeps leads to the solution of a large number of large-scale systems of equations. Several methods have been developed for alleviating this computational burden by approximating the frequency response functions. Among these, interpolatory model order reduction methods are perhaps the most successful. This paper reviews this family of approximation methods with particular attention to their applicability to specific classes of frequency response problems and their performance. It also includes novel aspects pertaining to the iterative solution of large-scale systems of equations in the context of model order reduction and frequency sweeps. All reviewed computational methods are illustrated with realistic, large-scale structural dynamic, acoustic, and vibro-acoustic analyses in wide frequency bands. These highlight both the potential of these methods for reducing CPU time and their limitations. Copyright © 2012 John Wiley & Sons, Ltd.

In single point incremental forming (SPIF), the sheet is incrementally deformed by a small spherical tool following a lengthy tool path. The simulation by the finite element method of SPIF requires extremely long computing times that limit the application to simple academic cases. The main challenge is to perform thousands of load increments modelling the lengthy tool path with elements that are small enough to model the small contact area. Because of the localised deformation in the process, a strong nonlinearity is observed in the vicinity of the tool. The rest of the sheet experiences an elastic deformation that introduces only a weak nonlinearity because of the change of shape. The standard use of the implicit time integration scheme is inefficient because it applies an iterative update (Newton–Raphson) strategy for the entire system of equations. The iterative update is recommended for the strong nonlinearity that is active in a small domain but is not required for the large part with only weak nonlinearities. It is proposed in this paper to split the finite element mesh into two domains. The first domain models the plastically deforming zone that experiences the strong nonlinearity. It applies a full nonlinear update for the internal force vector and the stiffness matrix every iteration. The second domain models the large elastically deforming zone of the sheet. It applies a pseudolinear update strategy based on a linearization at the beginning of each increment.Within the increment, it reuses the stiffness matrix and linearly updates the internal force vector. The partly linearized update strategy is cheaper than the full nonlinear update strategy, resulting in a reduction of the overall computing. Furthermore, in this paper, adaptive refinement is combined with the two domain method. It results in accelerating the standard SPIF implicit simulation of 3200 shell elements by a factor of 3.6.

SUMMARY The dual-primal finite element tearing and interconnecting method (FETI-DP) is extended to systems of linear equations arising from a finite element discretization for a class of fluid–structure interaction problems in the frequency domain. A preconditioned generalized minimal residual method is used to solve the linear equations for the Lagrange multipliers introduced on the subdomain boundaries to enforce continuity of the solution. The coupling between the fluid and the structure on the fluid–structure interface requires an appropriate choice of coarse level degrees of freedom in the FETI-DP algorithm to achieve fast convergence. Several choices are proposed and tested by numerical experiments on three-dimensional fluid–structure interaction problems in the mid-frequency regime that demonstrate the greatly improved performance of the proposed algorithm over the standard FETI-DP method. Copyright © 2011 John Wiley & Sons, Ltd.

SUMMARY Three preconditioners proposed by Eriksson, Choi and Merkel, and Turkel are implemented in a 2D upwind Euler flow solver on unstructured meshes. The mathematical formulations of these preconditioning schemes for different sets of primitive variables are drawn, and their eigenvalues and eigenvectors are compared with each other. For this purpose, these preconditioning schemes are expressed in a unified formulation. A cell-centered finite volume Roe's method is used for the discretization of the preconditioned Euler equations. The accuracy and performance of these preconditioning schemes are examined by computing steady low Mach number flows over a NACA0012 airfoil and a two-element NACA4412–4415 airfoil for different conditions. The study shows that these preconditioning schemes greatly enhance the accuracy and convergence rate of the solution of low Mach number flows. The study indicates that the preconditioning methods implemented provide nearly the same results in accuracy; however, they give different performances in convergence rate. It is demonstrated that although the convergence rate of steady solutions is almost independent of the choice of primitive variables and the structure of eigenvectors and their orthogonality, the condition number of the system of equations plays an important role, and it determines the convergence characteristics of solutions. Copyright © 2011 John Wiley & Sons, Ltd.

SUMMARY A high-order local transmitting boundary to model the propagation of acoustic or elastic, scalar or vectorvalued waves in unbounded domains of arbitrary geometry is proposed. It is based on an improved continuedfraction solution of the dynamic stiffness matrix of an unbounded medium. The coefficient matrices of the continued-fraction expansion are determined recursively from the scaled boundary finite element equation in dynamic stiffness. They are normalised using a matrix-valued scaling factor, which is chosen such that the robustness of the numerical procedure is improved. The resulting continued-fraction solution is suitable for systems with many DOFs. It converges over the whole frequency range with increasing order of expansion and leads to numerically more robust formulations in the frequency domain and time domain for arbitrarily high orders of approximation and large-scale systems. Introducing auxiliary variables, the continued-fraction solution is expressed as a system of linear equations in i! in the frequency domain. In the time domain, this corresponds to an equation of motion with symmetric, banded and frequency-independent coefficient matrices. It can be coupled seamlessly with finite elements. Standard procedures in structural dynamics are directly applicable in the frequency and time domains. Analytical and numerical examples demonstrate the superiority of the proposed method to an existing approach and its suitability for time-domain simulations of large-scale systems. Copyright © 2011 John Wiley & Sons, Ltd. Received 6 January 2011; Revised 29 March 2011; Accepted 20 April 2011

Extensions of deflation techniques previously developed for the Poisson equation to static elasticity are presented. Compared to the (scalar) Poisson equation (J. Comput. Phys. 2008; 227(24):10196–10208; Int. J. Numer. Meth. Engng 2010; DOI: 10.1002/nme.2932; Int. J. Numer. Meth. Biomed. Engng 2010; 26(1):73–85), the elasticity equations represent a system of equations, giving rise to more complex low-frequency modes (Multigrid. Elsevier: Amsterdam, 2000). In particular, the straightforward extension from the scalar case does not provide generally satisfactory convergence. However, a simple modification allows to recover the remarkable acceleration in convergence and CPU time reached in the scalar case. Numerous examples and timings are provided in a serial and a parallel context and show the dramatic improvements of up to two orders of magnitude in CPU time for grids with moderate graph depths compared to the non-deflated version. Furthermore, a monotonic decrease of iterations with increasing subdomains, as well as a remarkable acceleration for very few subdomains are also observed if all the rigid body modes are included. Copyright © 2011 John Wiley & Sons, Ltd.

A variationally consistent alpha finite element method (VC αFEM) for quadrilateral isoparametric elements is presented by constructing an assumed strain field in which the gradient of the compatible strain field is scaled with a free parameter α. The assumed strain field satisfies the orthogonal condition and the Hellinger–Reissner variational principle is used to establish the discretized system of equations. It is shown that the strain energy is a second-order continuous function of α, and the VC αFEM can produce both lower and upper bounds to the exact solution in the strain energy for all elasticity problems by choosing a proper α∈[0, αupper]. Based on this bound property, an exact-α approach has been devised to give an ultra-accurate solution that is very close to the exact one in the strain energy. Furthermore, the exact-α approach also works well for volumetric locking problems, by simply replacing the strain gradient matrix by a stabilization matrix. In addition, a regularization-α approach has also been suggested to overcome possible hourglass instability. Intensive numerical studies have been conducted to confirm the properties of the present VC αFEM, and a very good performance has been found in comparing to a large number of existing FEM models. Copyright © 2010 John Wiley & Sons, Ltd.

In this work a fast solver for large-scale three-dimensional elastodynamic crack problems is presented, implemented, and tested. The dual boundary element method in the Laplace transform domain is used for the accurate dynamic analysis of cracked bodies. The fast solution procedure is based on the use of hierarchical matrices for the representation of the collocation matrix for each computed value of the Laplace parameter. An ACA (adaptive cross approximation) algorithm is used for the population of the low rank blocks and its performance at varying Laplace parameters is investigated. A preconditioned GMRES is used for the solution of the resulting algebraic system of equations. The preconditioners are built exploiting the hierarchical arithmetic and taking full advantage of the hierarchical format. An original strategy, based on the computation of some local preconditioners only, is presented and tested to further speed up the overall analysis. The reported numerical results demonstrate the effectiveness of the technique for both uncracked and cracked solids and show significant reductions in terms of both memory storage and computational time.

The coupling of flexible structures to incompressible fluids draws a lot attention during the last decade. Many different solution schemes have been proposed. In this contribution we concentrate on strong coupling fluid-structure interaction by means of monolithic solution schemes. Therein, a Newton-Krylov method is applied to the monolithic set of nonlinear equations. Such schemes require good preconditioning to be efficient. We propose two preconditioners that apply algebraic multigrid techniques to the entire fluid-structure interaction system of equations. The first is based on a standard block Gauss-Seidel approach where approximate inverses of the individual field blocks are based on a algebraic multigrid hierarchy tailored for the type of the underlying physical problem. The second is based on a monolithic coarsening scheme for the coupled system that makes use of prolongation and restriction projections constructed for the individual fields. The resulting nonsymmetric monolithic algebraic multigrid method therefore involves coupling of the fields on coarse approximations to the problem yielding significantly enhanced performance.

We present stochastic projection schemes for approximating the solution of a class of deterministic linear elliptic partial differential equations defined on random domains. The key idea is to carry out spatial discretization using a combination of finite element methods and stochastic mesh representations. We prove a result to establish the conditions that the input uncertainty model must satisfy to ensure the validity of the stochastic mesh representation and hence the well posedness of the problem. Finite element spatial discretization of the governing equations using a stochastic mesh representation results in a linear random algebraic system of equations in a polynomial chaos basis whose coefficients of expansion can be non-intrusively computed either at the element or the global level. The resulting randomly parametrized algebraic equations are solved using stochastic projection schemes to approximate the response statistics. The proposed approach is demonstrated for modeling diffusion in a square domain with a rough wall and heat transfer analysis of a three-dimensional gas turbine blade model with uncertainty in the cooling core geometry. The numerical results are compared against Monte–Carlo simulations, and it is shown that the proposed approach provides high-quality approximations for the first two statistical moments at modest computational effort.

SUMMARY Population balance systems are models for processes in nature and industry that lead to a coupled system of equations (Navier–Stokes equations, transport equations, etc.) where the equations are defined in domains with different dimensions. This paper will study the impact of using different schemes for solving the three-dimensional (3D) equation of a precipitation process in a two-dimensional flow domain. The numerical schemes for the 3D equation are assessed with respect to the median of the volume fraction of the particle size distribution and the computational costs. It turns out that in the case of a structured flow field with small variations in time all schemes give qualitatively the same results. For a highly timedependent flow field, the evolution of the median of the volume fraction differs considerably between first order and higher order schemes. Copyright q 2010 John Wiley & Sons, Ltd.

SUMMARY In the computational electromagnetism community it is known how the differential formulation of an eddy-current problem, can be translated into a finite dimensional system of equations involving circulations and fluxes, by means of the so-called Discrete Geometric Approach. This is done by exploiting the geometric structure behind Maxwell’s equations. In this paper, we will show how the same Discrete Geometric Approach can be profitably used also to discretize an eddy-current problem formulated in an integral way. We rely on a purely geometric definition of a novel set of face vector basis functions that we use to construct the discrete counterparts—matrices—of both the Ohm’s constitutive relation and of the integral relation between the magnetic vector potential and the eddy-current density vector. The symmetry and positive-definiteness of such matrices will be demonstrated and their geometric structure will be apparent. Copyright q 2010 John Wiley & Sons, Ltd.

A number of recent papers (see, e.g. (Int. J. Mech. Sci. 2007; 49:454–465; Eur. J. Mech. A/Solids 2008; 27:859–881; Eng. Struct. 2008; 30:664–674; Int. J. Mech. Sci. 2009; 51:179–191)) have shown that classical limit analysis can be extended to incorporate such important features as geometric non-linearity, softening and various so-called ductility constraints. The generic formulation takes the form of a challenging (nonconvex and nonsmooth) optimization problem referred to in the mathematical programming literature as a mathematical program with equilibrium constraints (MPEC). Similar to a classical limit analysis, the aim is to compute in a single step a bound (upper bound, in the case of the extended problem) to the maximum load. The solution algorithm so far proposed to solve the MPEC is to convert it into an iterative non-linear programming problem and attempts to process this using a standard non-linear optimizer. Motivated by the fact that no method is guaranteed to solve such MPECs and by the need to avoid the use of an optimization approach, which is unfamiliar to most practising engineers, we propose, in the present paper, a novel numerical scheme to solve the MPEC as a constrained non-linear system of equations. We illustrate the application of this approach using the simple class of elastoplastic softening skeletal structures for which certain ductility conditions are prescribed. Copyright © 2009 John Wiley & Sons, Ltd.

The Riemann problem for an irreversibly compressible medium has been solved. This solution introduces the maximum medium density that is attained in the process of active loading. To close the system of equations an extra equation is required to describe the transfer of this parameter through the Eulerian mesh. The possible wave configurations have been analyzed and the corresponding equations for the evaluation of the contact pressure and velocity have been obtained. The existence and uniqueness of the solution have been proven. The technique of the Riemann problem's solution for the arbitrary Lagrange Euler mesh was developed. Examples of the Riemann problem solution for various wave configurations show that neglecting the bulk elastic plastic deformations yields significant errors in the results both quantitatively and qualitatively. Copyright © 2009 John Wiley & Sons, Ltd.

SUMMARY A multiscale computational framework is presented that provides a coupled self-consistent system of equations involving molecular mechanics at small scales and quasi-continuum mechanics at large scales. The proposed method permits simultaneous resolution of quasi-continuum and atomistic length scales and the associated displacement fields in a unified manner. Interatomic interactions are incorporated into the method through a set of analytical equations that contain nanoscale-based material moduli. These material moduli are defined via internal variables that are functions of the local atomic configuration parameters. Point defects like vacancy defects in nanomaterials perturb the atomic structure locally and generate localized force fields. Formation energy of vacancy is evaluated via interatomic potentials and minimization of this energy leads to nanoscale force fields around defects. These nanoscale force fields are then employed in the multiscale method to solve for the localized displacement fields in the vicinity of vacancies and defects. The finite element method that is developed based on the hierarchical multiscale framework furnishes a two-level statement of the problem. It concurrently feeds information at the molecular scale, formulated in terms of the nanoscale material moduli, into the quasi-continuum equations. Representative numerical examples are shown to validate the model and demonstrate its range of applicability. Copyright q 2008 John Wiley & Sons, Ltd.

A method to solve steady linear groundwater flow problems using generalized Fourier Series is developed and particularized for multiple Fourier series in two-dimensional domains. It leads to a linear vector equation whose solution provides a finite number of generalized Fourier coefficients approximating the hydraulic head field. Its implementation is shown and two relevant properties are found for the system matrix. It is always symmetric and, once computed, if additional Fourier terms are needed for a better approximation of the hydraulic head field, previously computed matrix elements remain invariant, i.e. only new rows and columns are added to the system matrix. The method is demonstrated in three simple cases with different geometries and transmissivity fields, where solutions are compared with analytical and finite element method results. Thus, the method is verified as an alternative to other flow solvers. Additionally, it provides a direct way to obtain the spectral form of the flow equation solution, given a spectral representation of transmissivity, and can be easily extended to obtain continuous velocity fields and their approximated spectral expressions. Copyright © 2009 John Wiley & Sons, Ltd.

This Study is concerned with the accurate computation of the unknown forces applied on the boundary of an elastic body using its measured displacement data with noise. Vision-based minimally intrusive force-sensing using elastically deformable grasping tools is the motivation for undertaking this problem. Since this problem involves incomplete and inconsistent displacement/force of ail elastic body, it leads to an ill-posed problem known as Cauchy's problem in elasticity. Vision-based displacement measurement necessitates large displacements of the elastic body for reasonable accuracy. Therefore, We use geometrically non-linear modelling of the elastic body, which was not considered by others who attempted to solve Cauchy's elasticity problem before. We present two methods to solve the problem. The first method uses the pseudo-inverse of an over-constrained system of equations. This method is shown to be not effective when the noise in the measured displacement data is high. We attribute this to the appearance of spurious forces at regions where there should not be any forces. The second method focuses on minimizing the spurious forces by varying the measured displacements within the known accuracy of the measurement technique. Both continuum and frame elements are used in the finite element modelling of the elastic bodies considered in the numerical examples. The performance of the two methods is compared using seven numerical examples, all of which show that the second method estimates the forces with an error that is not more than the noise in the measured displacements. An experiment was also conducted to demonstrate the effectiveness of the second method in accurately estimating the applied forces.

Recent advances in the literature regularize the strain-softening response of force-based frame elements by either modifying the constitutive parameters or scaling selected integration weights. Although the former case maintains numerical accuracy for strain-hardening behavior, the regularization requires a tight coupling of the element constitutive properties and the numerical integration method. In the latter case, objectivity is maintained for strain-softening problems; however, there is a lack of convergence for strain-hardening response. To resolve the dichotomy between strain-hardening and strain-softening solutions, a numerically consistent regularization technique is developed for force-based frame elements using interpolatory quadrature with two integration points of prescribed characteristic lengths at the element ends. Owing to manipulation of the integration weights at the element ends, the solution of a Vandermonde system of equations ensures numerical accuracy in the linear-elastic range of response. Comparison of closed-form solutions and published experimental results of reinforced concrete columns demonstrates the effect of the regularization approach on simulating the response of structural members. Copyright © 2008 John Wiley & Sons, Ltd.

The recent advances in microarchitectural bone imaging disclose the possibility to assess both the apparent density and the trabecular microstructure of intact bones in a single measurement. Coupling these imaging possibilities with microstructural finite element (µFE) analysis offers a powerful tool to improve bone stiffness and strength assessment for individual fracture risk prediction. Many elements are needed to accurately represent the intricate microarchitectural structure of bone; hence, the resulting µFE models possess a very large number of degrees of freedom. In order to be solved quickly and reliably on state-of-the-art parallel computers, the µFE analyses require advanced solution techniques. In this paper, we investigate the solution of the resulting systems of linear equations by the conjugate gradient algorithm, preconditioned by aggregation-based multigrid methods. We introduce a variant of the preconditioner that does not need assembling the system matrix but uses element-by-element techniques to build the multilevel hierarchy. The preconditioner exploits the voxel approach that is common in bone structure analysis, and it has modest memory requirements, at the same time robust and scalable. Using the proposed methods, we have solved in 12min a model of trabecular bone composed of 247 734 272 elements, yielding a matrix with 1 178 736 360 rows, using 1024 CRAY XT3 processors. The ability to solve, for the first time, large biomedical problems with over 1 billion degrees of freedom on a routine basis will help us improve our understanding of the influence of densitometric, morphological, and loading factors in the etiology of osteoporotic fractures such as commonly experienced at the hip, spine, and wrist. Copyright © 2007 John Wiley & Sons, Ltd.

In his article a special form of gradient elasticity is presented that can be used to describe wave dispersion. This new format of gradient elasticity is an appropriate dynamic extension of the earlier static counterpart of the gradient elasticity theory advocated in the early 1990s by Aifantis and co-workers. In order to capture dispersion of propagating waves, both higher-order inertia and higher-order stiffness contributions are included, a fact which implies (and is denoted as) dynamic consistency. The two higher-order terms are accompanied by two associated length scales. To facilitate finite element implementations, the model is rewritten such that 0-continuity of the interpolation is sufficient. An auxiliary displacement field is introduced which allows the original fourth-order equations to be split into two coupled sets of second-order equations. Positive-definiteness of the kinetic energy requires that the inertia length scale is larger than the stiffness length scale. The governing equations, boundary conditions and the discretized system of equations are presented. Finally, dispersive wave propagation in a one-dimensional bar is considered in a numerical example. Copyright © 2007 John Wiley & Sons, Ltd.

The hydrodynamics of planing ships are studied using a finite pressure element method. In this method, a boundary value problem (BVP) is formulated based on linear planing theory; the planing ship is represented by the pressure distribution acting on the wetted bottom of the ship, and the magnitude of this pressure distribution is evaluated using a boundary element method. The pressure is discretized using overlapping pressure pyramids, known as tent functions, so that the resulting distribution is piece-wise continuous in both longitudinal and transverse directions. A set of linear algebraic equations for the determination of the pressure is then established using a collocation technique. It is found that the matrix of the linear equations is ill conditioned; this leads to oscillatory behaviour of the predicted pressure distribution if the direct solution method of LU decomposition or Gaussian elimination is used to solve the system of linear equations. In the current study, this numerical instability is analysed in detail. It is found that the problem can be addressed by adopting singular value decomposition (SVD) technique for the solution of the linear equations. Using this method, the hydrodynamic results for flat-bottomed and prismatic planing ships are calculated and a good agreement is demonstrated with Savitsky's empirical relations.

The aim of this paper is to develop non-linear recursive formulations for decoupling joint and elastic accelerations, while maintaining the non-linear inertia coupling between rigid body motion and elastic deformation in deformable mechanical systems. The inertia projection schemes used in most existing recursive formulations for the dynamic analysis of deformable mechanisms lead to dense coefficient matrices in the equations of motion. Consequently, there are strong dynamic couplings between the joint and elastic coordinates. When the number of elastic degrees of freedom increases, the size of the coefficient matrix in the equations of motion becomes large. Consequently, the use of these recursive formulations for solving the joint and elastic accelerations becomes less efficient. In this paper, the non-linear recursive formulations have been used to decouple the elastic and joint accelerations in deformable mechanical systems. The relationships between the absolute, elastic and joint variables and generalized Newton–Euler equations are used to develop systems of loosely coupled equations that have sparse matrix structure. By using the inertia matrix structure of deformable mechanical systems and the fact that joint reaction forces associated with elastic coordinates do represent independent variables, a reduced system of equations whose dimension is dependent of the number of elastic degrees of freedom is obtained. This system can be solved for the joint accelerations as well as for the joint reaction forces. The use of the approaches developed in this investigation is illustrated using deformable open-loop serial robot and closed-loop four-bar mechanical systems. Copyright © 2007 John Wiley & Sons, Ltd.

We propose the use of reduced-order models to accelerate the solution of systems of equations using iterative solvers in time stepping schemes for large-scale numerical simulation. The acceleration is achieved by determining an improved initial guess for the iterative process based on information in the solution vectors from previous time steps. The algorithm basically consists of two projection steps: (1) projecting the governing equations onto a subspace spanned by a low number of global empirical basis functions extracted from previous time step solutions, and (2) solving the governing equations in this reduced space and projecting the solution back on the original, high dimensional one. We applied the algorithm to numerical models for simulation of two-phase flow through heterogeneous porous media. In particular we considered implicit-pressure explicit-saturation (IMPES) schemes and investigated the scope to accelerate the iterative solution of the pressure equation, which is by far the most time-consuming part of any IMPES scheme. We achieved a substantial reduction in the number of iterations and an associated acceleration of the solution. Our largest test problem involved 93 500 variables, in which case we obtained a maximum reduction in computing time of 67%. The method is particularly attractive for problems with time-varying parameters or source terms. Copyright © 2006 John Wiley & Sons, Ltd.

Numerical issues arising when integrating hyperelasto-plastic constitutive equations with elastic anisotropy, stemming from anisotropic damage, and plastic anisotropy as represented by kinematic hardening are discussed. In particular, solution algorithms for the corresponding non-linear system of equations due to implicit integration algorithms are addressed. It is shown that algorithms like staggered iteration and quasi-Newton techniques are superior to a pure Newton technique when the cpu time is compared. However, the drawback of the staggered and the quasi-Newton technique is that rather small time steps must be taken to ensure convergence, which can be of importance when applying complex constitutive models in a finite element programme. Copyright © 2005 John Wiley & Sons, Ltd.

In this work, we compare the fully implicit (FI) and implicit pressure-explicit saturation (IMPES) formulations for the simulation of water injection in fractured media. The system of partial differential equations is discretized within the discrete-fracture framework using a control-volume method. A unique feature of the methodology is that there is no need for the computation of matrix–fracture transfer terms. The non-linear system of equations resulting from the FI formulation is solved with state-of-the-art Newton and tensor methods. Direct and Krylov iterative methods are employed to solve the system resulting from the Newton linearization. The performance of the FI and IMPES formulations is compared with numerical testing. Results show that the contrast between matrix and fracture properties affects the performance of both IMPES and FI formulations and that the tensor method outperforms all the Newton solvers for the near-singular Jacobian matrix. Copyright © 2006 John Wiley & Sons, Ltd.

In this paper, a new finite element method is described and applied. It is based on a theory developed to model poromechanical problems where the mechanical part is obeying a second gradient theory. The aim of such a work is to properly model the post localized behaviour of soils and rocks saturated with a pore fluid. Beside the development of this new coupled theory, a corresponding finite element method has been developed. The elements used are based on a weak form of the relation between the deformation gradient and the second gradient, using a field of Lagrange multipliers. The global problem is solved by a system of equations where the kinematic variables are fully coupled with the pore pressure. Some numerical experiments showing the effectiveness of the method ends the paper.

In general, the derivative of an eigenvector of a vibrating symmetric system is the solution of a singular problem. Further complications are encountered in dealing with asymmetric damped systems for which the left and right eigenvectors and their derivatives become distinct and complex. Several approaches have been proposed to overcome this singularity such as Nelson's method and the modal method. In the present work, a new approach is presented for calculating simultaneously the derivatives of the eigenvalues and their associated derivatives of the left and right eigenvectors for asymmetric damped systems. With the proposed method, the exact eigenderivatives can be obtained by solving a first-order linear algebraic system of equations. The method is applied on a 104 DOF ventilator–rotor system, which is used as an example of an asymmetric damped system with distinct eigenvalues. The diameter of the shaft has been chosen as the design parameter. The comparison of the computational time shows that the proposed method is more efficient than both Nelson's approach and the modal method. Copyright © 2006 John Wiley & Sons, Ltd.

We examine a bubble-stabilized finite element method for enforcing Dirichlet constraints on embedded interfaces. By ‘embedded’ we refer to problems of general interest wherein the geometry of the interface is assumed independent of some underlying bulk mesh. As such, the robust imposition of Dirichlet constraints using a Lagrange multiplier field is not trivial. To focus issues, we consider a simple one-sided problem that is representative of a wide class of evolving-interface problems. The bulk field is decomposed into coarse and fine scales, giving rise to coarse-scale and fine-scale one-sided sub-problems. The fine-scale solution is approximated with bubble functions, permitting static condensation and giving rise to a stabilized form bearing strong analogy with a classical method. Importantly, the method is simple to implement, readily extends to multiple dimensions, obviates the need to specify any free stabilization parameters, and can lead to a symmetric, positive-definite system of equations. The performance of the method is then examined through several numerical examples. The accuracy of the Lagrange multiplier is compared to results obtained using a local version of the domain integral method. The variational multiscale approach proposed herein is shown to both stabilize the Lagrange multiplier and improve the accuracy of the post-processed fluxes. Copyright © 2006 John Wiley & Sons, Ltd.

The FETI-DP and BDDC algorithms are reformulated using Block Cholesky factorizations, an approach which can provide a useful framework for the design of domain decomposition algorithms for solving symmetric positive definite linear system of equations. Instead of introducing Lagrange multipliers to enforce the coarse level, primal continuity constraints in these algorithms, a change of variables is used such that each primal constraint corresponds to an explicit degree of freedom. With the new formulation of these algorithms, a simplified proof is provided that the spectra of a pair of FETI-DP and BDDC algorithms, with the same set of primal constraints, are essentially the same. Numerical experiments for a two-dimensional Laplace's equation also confirm this result. Copyright © 2005 John Wiley & Sons, Ltd.

SUMMARY In the past two decades, considerable improvements concerning integration algorithms and solvers involved in boundary-element formulations have been obtained. First, a great deal of efficient techniques for evaluating singular and quasi-singular boundary-element integrals have been, definitely, established, and second, iterative Krylov solvers have proven to be advantageous when compared to direct ones also including non-Hermitian matrices. The former fact has implied in CPU-time reduction during the assembling of the system of equations and the latter fact in its faster solution. In this paper, a triangle-polar-co-ordinate transformation and the Telles co-ordinate transformation, applied in previous works independently for evaluating singular and quasi-singular integrals, are combined to increase the efficiency of the integration algorithms, and so, to improve the performance of the matrixassembly routines. In addition, the Jacobi-preconditioned biconjugate gradient (J-BiCG) solver is used to develop a generic substructuring boundary-element algorithm. In this way, it is not only the system solution accelerated but also the computer memory optimized. Discontinuous boundary elements are implemented to simplify the coupling algorithm for a generic number of subregions. Several numerical experiments are carried out to show the performance of the computer code with regard to matrix assembly and the system solving. In the discussion of results, expressed in terms of accuracy and CPU time, advantages and potential applications of the BE code developed are highlighted. Copyright 2006 John Wiley & Sons, Ltd.

In this work, we determine the different patterns of possible wave structures that can be observed on a thin film flowing on an inclined plane when at the free surface a shear force (surface traction) is applied. Different wave structures are obtained dependening on the selected combination of downstream and upstream boundary conditions and initial conditions. The resulting initial boundary value problems are solved numerically using the direct BEM numerical solution of the complete two-dimensional Stokes system of equations. In our numerical results, the initial discontinuous shock profiles joining uniform fluid depths are smoothed due to the two-dimensional character of the Stokes formulation, including the effect of the gravitational force as well as the interfacial surface tension force. In this way, physically feasible continuous surface profiles are determined, in which the initial uniform depths are joined by smooth moving wave structures. Numerical solutions have been attained to reproduce the different patterns of possible wave structures previously reported in the literature and extended to identify some other new structures and features defining the behaviour of the surface patterns. Copyright © 2006 John Wiley & Sons, Ltd.

This paper presents a fast formulation of the hybrid boundary node method (Hybrid BNM) for solving problems governed by Laplace's equation in 3D. The preconditioned GMRES is employed for solving the resulting system of equations. At each iteration step of the GMRES, the matrix–vector multiplication is accelerated by the fast multipole method. Green's kernel function is expanded in terms of spherical harmonic series. An oct-tree data structure is used to hierarchically subdivide the computational domain into well-separated cells and to invoke the multipole expansion approximation. Formulations for the local and multipole expansions, and also conversion of multipole to local expansion are given. And a binary tree data structure is applied to accelerate the moving least square approximation on surfaces. All the formulations are implemented in a computer code written in C++. Numerical examples demonstrate the accuracy and efficiency of the proposed approach. Copyright © 2005 John Wiley & Sons, Ltd.

The dual-primal finite element tearing and interconnecting (FETI-DP) domain decomposition method (DDM) is extended to address the iterative solution of a class of indefinite problems of the form (K-σ 2 M)u = f, and a class of complex problems of the form (K-σ 2 M+iσD)u = f, where K, M, and D are three real symmetric matrices arising from the finite element discretization of solid and shell dynamic problems, i is the imaginary complex number, and a is a real positive number. A key component of this extension is a new coarse problem based on the free-space solutions of Navier's equations of motion. These solutions are waves, and therefore the resulting DDM is reminiscent of the FETI-H method. For this reason, it is named here the FETI-DPH method. For a practically large a range, FETI-DPH is shown numerically to be scalable with respect to all of the problem size, substructure size, and number of substructures. The CPU performance of this iterative solver is illustrated on a 40-processor computing system with the parallel solution, for various a ranges, of several large-scale, indefinite, or complex-valued systems of equations associated with shifted eigenvalue and forced frequency response structural dynamics problems.

This paper presents a new method for solving any combination of linear–non-linear equations. The method is based on the separation of linear equations in terms of some selected variables from the non-linear ones. The linear group is solved by means of any method suitable for the linear system. This operation needs no iteration. The non-linear group, however, is solved by an iteration technique based on a new formula using the Taylor series expansion. The method has been described and demonstrated in several examples of analytical systems with very good results. The new method needs the initial approximations for non-linear variables only. This requires far less computation than the Newton–Raphson method. The method also has a very good convergence rate. The proposed method is most beneficial for engineering systems that very often involve a large number of linear equations with limited number of non-linear equations. Copyright © 2005 John Wiley & Sons, Ltd.

SUMMARY In this work we evaluate the performance of the Left Conjugate Direction method recently introduced by Yuan, Golub, Plemmons and Cec¶‡lio for the solution of nonsymmetric systems of linear equations arising from the implicit semi-discrete SUPG flnite element formulation of advective-difiusive and inviscid compressible ∞ows. We extend the original algorithm to accommodate restarts and typical element-by-element preconditioners. We also show how to select the flrst left conjugate vector to start LCD. Several problems are solved, accessing performance parameters such as number of iterations, memory requirements and CPU times, and results are compared with other algorithms, such as GMRES, TFQMR and Bi-CGSTAB. Copyright c ∞ 2003 John Wiley & Sons, Ltd.