Variational boundary-integral-equation approach to unilateral contact problems in elasticity

On the basis of the boundary integral equation method, three variational principles for the frictionless unilateral contact problem in elasticity are presented. Two of them are saddle-point principles for the boundary unknowns (including the contact displacements); a third one is a maximum principle for the unknown contact displacements only. A discretization by boundary elements leads to algebraic formulations in the shape either of quadratic programming problems, or of linear complementarity problems, all characterized by symmetry and sign definiteness of the coefficient matrices. The method is also applicable to contact problems between two uncompenetrable elastic solids, as well as to the crack problem of fracture mechanics.     

A numerical procedure for multiple circular holes and elastic inclusions in a finite domain with a circular boundary

This paper describes a numerical procedure for solving two-dimensional elastostatics problems with multiple circular holes and elastic inclusions in a finite domain with a circular boundary. The inclusions may have arbitrary elastic properties, different from those of the matrix, and the holes may be traction free or loaded with uniform normal pressure. The loading can be applied on all or part of the finite external boundary. Complex potentials are expressed in the form of integrals of the tractions and displacements on the boundaries. The unknown boundary tractions and displacements are approximated by truncated complex Fourier series. A linear algebraic system is obtained by using Taylor series expansion without boundary discretization. The matrix of the linear system has diagonal submatrices on its diagonal, which allows the system to be effectively solved by using a block Gauss-Seidel iterative algorithm.     

A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method

The Element free Galerkin method, which is based on the Moving Least Squares approximation, requires only nodal data and no element connectivity, and therefore is more flexible than the conventional finite element method. Direct imposition of essential boundary conditions for the element free Galerkin (EFG) method is always difficult because the shape functions from the Moving Least Squares approximation do not have the delta function property. In the prior literature, a direct collocation of the fictitious nodal values & u circ; used as undetermined coefficients in the MLS approximation, u ^h (x) [u ^h (x)=Φ·& u circ;], was used to enforce the essential boundary conditions. A modified collocation method using the actual nodal values of the trial function u ^h (x) is presented here, to enforce the essential boundary conditions. This modified collocation method is more consistent with the variational basis of the EFG method. Alternatively, a penalty formulation for easily imposing the essential boundary conditions in the EFG method with the MLS approximation is also presented. The present penalty formulation yields a symmetric positive definite system stiffness matrix. Numerical examples show that the present penalty method does not exhibit any volumetric locking and retains high rates of convergence for both displacements and strain energy. The penalty method is easy to implement as compared to the Lagrange multiplier method, which increases the number of degrees of freedom and yields a non-positive definite system matrix.     

A new boundary element approach for contact problems with friction

A new boundary element solution algorithm for two-dimensional and axisymmetric contact problems with friction, based on an independent discretization of the contacting surfaces and under static and proportional loading conditions, is presented. The solution procedure uses the element shape functions to distribute the geometry, tractions and displacements on each contact element. The contact constraints are then applied between each contacting node and the opposite contact segment. The overall boundary element matrix equations for the contacting bodies are coupled using the contact conditions at the interface without introducing any additional variables into the solution matrix. The algorithm is applied to several two-dimensional and axisymmetric frictional contact examples and the results obtained are in very good agreement with finite element and analytical solutions.     

Receding contact problem involving large displacements using the BEM

A new algorithm is presented for the boundary element analysis of the two-dimensional contact problem between elastic solids involving large displacements. The contact constraints are not applied node-on-node but node-on-element, using the element shape functions to distribute the geometry, displacements and tractions on each element in the contact zone. Thus, the discretizations performed along the two surfaces in contact need not necessarily be the same. The solution procedure is based on the updated Lagrangian approach and the resulting method is incremental. The algorithm guarantees equilibrium and compatibility at the nodes in the final deformed configuration and allows us to deal with problems undergoing large displacements without it being necessary to change the initial discretization of the boundary of the bodies. Only the frictionless static problem is dealt with, and the proposed algorithm is applied to the most representative receding contact problem: a layer pressed against an elastic foundation. The results obtained when the displacements are small are in good agreement with the analytical solution. When large displacements are considered, another nonlinearity appears and its influence will be shown in this paper.     

Inverse problem of simultaneously estimating the thermal conductivity and boundary shape

An inverse analysis is used to simultaneously estimate the thermal conductivity and the boundary shape in steady-state heat conduction problems. The numerical scheme consists of a body-fitted grid generation technique to mesh the heat conducting body and solve the heat conduction equation – a novel, efficient, and easy to implement sensitivity analysis scheme to compute the sensitivity coefficients, and the conjugate gradient method as an optimization method to minimize the mismatch between the computed temperature distribution on some part of the body boundary and the measured temperatures. Using the proposed scheme, all sensitivity coefficients can be obtained in one solution of the direct heat conduction problem, irrespective of the large number of unknown parameters for the boundary shape. The obtained results reveal the accuracy, efficiency, and robustness of the proposed algorithm.     

Variation of mixed modes stress intensity factors of an inclined semi-elliptical surface crack

In this paper, a singular integral equation method is applied to calculate the stress intensity factor along crack front of a 3D inclined semi-elliptical surface crack in a semi-infinite body under tension. The stress field induced by displacement discontinuities in a semi-infinite body is used as the fundamental solution. Then, the problem is formulated as a system of integral equations with singularities of the form r ^–3. In the numerical calculation, the unknown body force doublets are approximated by the product of fundamental density functions and polynomials. The results show that the present method yields smooth variations of mixed modes stress intensity factors along the crack front accurately for various geometrical conditions. The effects of inclination angle, elliptical shape, and Poisson's ratio are considered in the analysis. Crack mouth opening displacements are shown in figures to predict the crack depth and inclination angle. When the inclination angle is 60 degree, the mode I stress intensity factor F _I has negative value in the limited region near free surface. Therefore, the actual crack surface seems to contact each other near the surface.     

Use of an extremal functional in solving for an unknown bottom surface given a free surface profile

In this paper the free surface fluid flow is induced by an obstacle on the bottom of the channel whose exact shape and location are unknown a priori. The inviscid fluid flow is assumed to be steady, incompressible and irrotational under the influence of gravity. A boundary integral technique, which is based on the combination of the boundary integral method (BIM), the variational principle technique (VPT) and the application of a minimization technique (MT) is developed. Two minimization techniques, namely the extremal pressure method (EPM) and the extremal energy method (EEM), are extensively used in identifying the unknown bottom surfaces. To illustrate these techniques the free surface profile to be applied in the inverse analysis has been generated following a direct formulation when the solid bottom boundary is monotonically decreasing, monotonically increasing, when it possesses a single hump and when it possesses a single depression. It is found that in all four cases both extremum techniques can be very accurately used to identify the shape and position of the solid obstacle.     

The computation of the J integral with the traction boundary integral equation

The solutions of the displacement boundary integral equation (BIE) are not uniquely determined in certain types of boundary conditions. Traction boundary integral equations that have unique solutions in these traction and mixed boundary cases are established. For two‐dimensional linear elasticity problems, the divergence‐free property of the traction boundary integral equation is established. By applying Stokes' theorem, unknown tractions or displacements can be reduced to computation of traction integral potential functions at the boundary points. The same is true of the J integral: it is divergence‐free and the evaluation of the J integral can be inverted into the computation of the J integral potential functions at the boundary points of the cracked body. The J integral can be expressed as the linear combination of the tractions and displacements from the traction BIE on the boundary of the cracked body. Numerical integrals are not needed at all. Selected examples are presented to demonstrate the validity of the traction boundary integral and J integral. Copyright © 2004 John Wiley & Sons, Ltd.     

Image reconstruction for a partially immersed imperfectly conducting cylinder by genetic algorithm

This article presents a computational approach to the imaging of a partially immersed imperfectly conducting cylinder. An imperfectly conducting cylinder of unknown shape and conductivity scatters the incident transverse magnetic (TM) wave in free space while the scattered field is recorded outside. Based on the boundary condition and the measured scattered field, a set of nonlinear integral equations, and the inverse scattering problem are reformulated into an optimization problem. We use genetic algorithm (GA) to reconstruct the shape and the conductivity of a partially immersed imperfectly conducting cylinder. The genetic algorithm is then used to find out the global extreme solution of the cost function. Numerical results demonstrated that, even when the initial guess is far away from the exact one, good reconstruction can be obtained. In such a case, the gradient‐based methods often get trapped in a local extreme. In addition, the effect of random noise on the reconstruction is investigated. © 2009 Wiley Periodicals, Inc. Int J Imaging Syst Technol, 19, 299–305, 2009     

Green's functions for Kirchhoff plates of irregular shape

A method is proposed for the construction of Green's functions for the Sophie Germain equation in regions of irregular shape with mixed boundary conditions imposed. The method is based on the boundary integral equation approach where a kernel vector function B satisfies the biharmonic equation inside the region. This leads to a regular boundary integral equation where the compensating loads and moments are applied to the boundary. Green's function is consequently expressed in terms of the kernel vector function B, the fundamental solution function of the biharmonic equation, and kernel functions of the inverse regular integral operators. To compute moments and forces, the kernel functions are differentiated under the integral sign. The proposed method appears highly effective in computing both displacements and stress components.     

Layerwise mixed element with sublaminates approximation and 3D zig‐zag field,for analysis of local effects in laminated and sandwich composites

The analysis of damaged sandwich composites and thick laminates remains a still challenging task, owing to their intricate distributions of displacements and stresses across the thickness. The arduous task is capturing these distributions with an affordable computational effort. The use of sublaminate zig‐zag models appears promising since they allow to group several physical layers into a sublaminate without violating the continuity of interlaminar stresses. To contribute to research in this field, a mixed, sublaminate element is developed in which displacements and stresses are approximated by two independent zig‐zag models. Like for any zig‐zag model, the interfacial contact conditions are fulfilled a priori suitably choosing the continuity functions and constants. In addition to the contact conditions on interlaminar shears, here those on the transverse normal stress and stress gradient are also fulfilled, as required by the elasticity theory. This uncustomary feature implies the presence of unwise terms which need C² continuous shape functions. To overcome this drawback, C¹ and C² terms are eliminated from the displacement field, as allowed by the mixed approach, then the interlaminar stresses obtained by the C⁰‐made displacement model are cast into a separate zig‐zag representation which restores their continuity at the interfaces. As nodal d.o.f. the three displacements and the three interlaminar stresses at the upper and lower faces are assumed, in order to fulfil the contact conditions assembling the sublaminates. Standard, C⁰, serendipity interpolation polynomials are used to represent both internal displacements and stress fields. As a consequence of this choice, the intra‐element stress equilibria are fulfilled in an approximate integral form. This makes easier the development of the element, but does not compromise its accuracy and convergence rate, as shown by former applications. A postprocessing procedure is developed which improves the estimation of interlaminar stresses and helps to reduce the number of sublaminates. The result is that the analysis can be carried out using a single computational layer across the thickness and a reasonably fine in‐plane discretization, even with distinctly different material properties of layers and in presence of damage. The overall computational time is reduced to an half of that required by a conventional sublaminate model, like present one deprived of zig‐zag terms and postprocessing procedure. The numerical applications concern thermally loaded and piezoelectrically actuated thick laminates and damaged sandwich panels. Copyright © 2006 John Wiley & Sons, Ltd.     

A SINGULAR INTEGRAL EQUATION SOLUTION FOR THE LINEAR ELASTIC CRACK OPENING DISPLACEMENT OF AN ARBITRARILY SHAPED PLANE CRACK: PART I FINITE-PART INTEGRAL SOLUTIONS

Abstract— The aim of the paper is to compute the local crack face displacements of a linear elastic body containing an arbitrarily shaped plane crack. From the crack face displacements the local stress intensity factors can be derived.
The boundary value problem for a plane crack of arbitrary shape, embedded in a linear elastic medium, has been treated by several authors by the singular integral equation (SIE) approach. Their computations lead to a set of hyper-singular integral equations for the Cartesian components of the unknown crack face displacements. To solve these equations the authors present a discretization procedure based on six-node triangular finite elements. A total set of 24 finite-part integrals defined over a triangular area can be developed. These 2D-finite-part integrals can be split into both a 1D-regular and a 1D-finite-part-integral by means of the polar coordinates so that they can be solved in closed form. Finally, the investigation of the SIEs is reduced to a discrete set of linear algebraic equations for the unknown nodal point values. The necessary steps will be demonstrated in detail. The derived closed-form solutions will be offered in the text and in the appendices.     

Obstacle identification using the TRAC algorithm with a second-order ABC

We consider obstacle identification using wave propagation. In such problems, one wants to find the location, shape, and size of an unknown obstacle from given measurements. We propose an algorithm for the identification task based on a time-reversed absorbing condition (TRAC) technique. Here, we apply the TRAC method to time-dependent linear acoustics, although our methodology can be applied to other wave-related problems as well, such as elastodynamics. There are two main contributions of our identification algorithm. The first contribution is the development of a robust and effective method for obstacle identification. While the original paper presented criteria for accepting or rejecting regions that enclose the obstacle, we use these criteria to develop an algorithm that automatically identifies the location of the obstacle. The second contribution is the utilization of an improved absorbing boundary condition (ABC) for the identification. We use the second-order Engquist-Majda ABC, and we implement it with a finite element scheme. To our knowledge, this is the first time that the second-order Engquist-Majda ABC is employed with the finite element method, as this boundary condition does not naturally fit in finite element schemes in its original form. Numerical experiments for the algorithms are presented.     

An advanced numerical method for computing elastodynamic fracture parameters in functionally graded materials

A new meshless method for computing the dynamic stress intensity factors (SIFs) in continuously non-homogeneous solids under a transient dynamic load is presented. The method is based on the local boundary integral equation (LBIE) formulation and the moving least squares (MLS) approximation. The analyzed domain is divided into small subdomains, in which a weak solution is assumed to exist. Nodal points are randomly spread in the analyzed domain and each one is surrounded by a circle centered at the collocation point. The boundary-domain integral formulation with elastostatic fundamental solutions for homogeneous solids in Laplace-transformed domain is used to obtain the weak solution for subdomains. On the boundary of the subdomains, both the displacement and the traction vectors are unknown generally. If modified elastostatic fundamental solutions vanishing on the boundary of the subdomain are employed, the traction vector is eliminated from the local boundary integral equations for all interior nodal points. The spatial variation of the displacements is approximated by the MLS scheme.     

A level set method for shape and topology optimization of large‐displacement compliant mechanisms

A parameterization level set method is presented for structural shape and topology optimization of compliant mechanisms involving large displacements. A level set model is established mathematically as the Hamilton–Jacobi equation to capture the motion of the free boundary of a continuum structure. The structural design boundary is thus described implicitly as the zero level set of a level set scalar function of higher dimension. The radial basis function with compact support is then applied to interpolate the level set function, leading to a relaxation and separation of the temporal and spatial discretizations related to the original partial differential equation. In doing so, the more difficult shape and topology optimization problem is now fully parameterized into a relatively easier size optimization of generalized expansion coefficients. As a result, the optimization is changed into a numerical process of implementing a series of motions of the implicit level set function via an existing efficient convex programming method. With the concept of the shape derivative, the geometrical non‐linearity is included in the rigorous design sensitivity analysis to appropriately capture the large displacements of compliant mechanisms. Several numerical benchmark examples illustrate the effectiveness of the present level set method, in particular, its capability of generating new holes inside the material domain. The proposed method not only retains the favorable features of the implicit free boundary representation but also overcomes several unfavorable numerical considerations relevant to the explicit scheme, the reinitialization procedure, and the velocity extension algorithm in the conventional level set method. Copyright © 2008 John Wiley & Sons, Ltd.     

基于PVM的网络并行子结构共轭梯度法

网络并行环境是近年来国际上并行环境的一个重要方向,ＰＶＭ是当前最流行的支持异构或同构型网络并行计算的软件平台之一。本文采用子结构共轭梯度法研究了基于ＰＶＭ的网络并行有限元,该方法将有限元网格划分为ｎ个子结构,再将ｎ个子结构的数据分送给网上ｎ台可用微机,ｎ台微机并行形成和组集ｎ个子结构的劲度矩阵和荷载列阵,然后采用预条件共轭梯度法并行求解结点位移,最后ｎ台微机并行对ｎ个子结构进行应变和应力分析。该方法不需形成结构的总体劲度矩阵和荷载列阵,可同时迭代求出所有结点位移,且比一般的迭代法收敛要快。算例表明此种并行子结构共轭梯度法在网络上能获得较高的并行加速比。     

SHAPE OPTIMIZATION IN SHELL THEORY: Design Sensitivity of the Continuous Problem

We consider a non-shallow shell made of an isotropic homogeneous material, working in linear elastic conditions, subjected to a given load. Our aim is to change the shape of the shell so that it resists better towards a given criterion. By shape, we mean essentially the midsurface of the shell. The thickness could be added without any difficulty. The important aspect that we study here is the midsurface.

This problem is worked by gradient type methods. We prove that if the criterion depends on the displacement field through a differentiable function, then it depends on the shape in a differentiable manner, because the displacement field is a differentiable function of the shape. Then we present an analytical formula giving the exact gradient of the criteria before any discretization. After that, we explain how to compute numerically an approximation to this exact gradient. Then we give numerical results.     

Shape optimization using the boundary element method with substructuring

Applications of the boundary element method for two- and three-dimensional structural shape optimization are presented. The displacements and stresses are computed using the boundary element method. Sub-structuring is used to isolate the portion of the structure undergoing geometric change. The corresponding non-linear programming problem for the optimization is solved by the generalized reduced gradient method. B-spline curves and surfaces are introduced to describe the shape of the design. The control points on these curves or surfaces are selected as design variables. The design objective may be either to minimize the weight or a peak stress of the component by determining the optimum shape subject to geometrical and stress constraints. The use of substructuring allows for problem solution without requiring traditional simplifications such as linearization of the constraints. The method has been successfully applied to the structural shape optimization of plane stress, plane strain and three-dimensional elasticity problems.