WO2009116048A2 - A method for cosserat point element (cpe) modeling of nonlinear elastic materials - Google Patents

Abstract

A method for Cosserat Point Element (CPE) modeling of nonlinear elastic materials, utilizing a strain energy function characterizing the response of an eight node 3-D brick element to homogeneous and inhomogeneous deformations, which includes full coupling of bending and torsional modes of deformation. The method includes determining the constitutive coefficients which satisfy bending and torsional modes of deformation for three specific reference element geometries and combining these solutions to obtain constitutive equations for general reference element shapes which ensure that the strain energy remains a positive definite function for all inhomogeneous deformations of the CPE. The constitutive coefficients characterize a generalized CPE model that predicts the response of a 3-D brick element for finite elastic deformations which is free of hourglass instabilities and locking for nearly incompressible materials, thereby providing a robust modeling tool.

Description

A METHOD FOR COSSERAT POINT ELEMENT (CPE) MODELING OF NONLINEAR ELASTIC MATERIALS

FIELD OF THE INVENTION

The present invention relates to a generalized form for the strain energy of inhomogeneous deformations, and more particularly to a 3-D brick Cosserat Point Element (CPE) which includes full coupling of bending and torsional modes of deformation, such that the resulting constitutive coefficients ensure that the strain energy for inhomogeneous deformations remains a positive definite function of the inhomogeneous strain measures for all reference element shapes. BACKGROUND OF THE INVENTION

1. Cosserat Point Element Theory

The Cosserat Point Element (CPE) is a finite element technology for the numerical solution of large deformation problems of elastic materials. Since the CPE can be introduced into general purpose finite element codes, it can be used to solve the whole host of problems considered in industry to design elastic structures. In particular, the CPE can be used for the problem of large elastic deformations of tires, which is particularly challenging, since the solutions based on other elements typically are corrupted by unphysical hourglass instabilities.

Recently, Nadler and Rubin (2003) developed a 3-D eight-noded brick element based on the theory of a Cosserat point (Rubin, 1995, 2000). This CPE can be used to formulate the numerical solution of dynamic problems for nonlinear hyperelastic materials. Fig. 1 is a prior art sketch of a general brick CPE showing the numbering of the nodes 110 and the surfaces 120. The kinematics of the CPE are characterized by eight element director vectors 130 and the kinetics propose eight balance laws of director momentum to determine the dynamic response of the element. The locations of the nodes in the current deformed configuration are characterized by eight nodal director vectors and the element directors are related to the nodal directors by standard tri-linear shape functions. Moreover, the CPE theory considers the element as a structure and introduces a strain energy function which characterizes the response of the structure. Also, the nodal forces are related to derivatives of the strain energy function through algebraic relations in a similar manner to the relationship of the stress to derivatives of the strain energy function in the full three- dimensional theory of hyperelastic materials.

In the standard finite element procedures for hyperelastic materials the response of the element is determined by integrals over the element region which assume that the kinematic approximation is valid pointwise. This is in contrast with the CPE which needs no integration over the element region. Furthermore, it is known (Zienkiewicz and Taylor, 2005) that the standard element based on full integration is robust but that it exhibits unphysical locking for thin structures with poor element aspect ratios and for nearly incompressible materials. Special methods based on enhanced strains, incompatible modes and reduced integration with hourglass control (e.g. Belytschko et al. 1984; Simo and Rifai, 1990; Simo and Armero, 1992; Belytschko and Bindeman, 1993; Simo et al. 1993; Bonet and Bhargava, 1995; Reese and Wriggers, 1996, 2000; Reese et al. 2000; Hutter et al. 2000) have been developed to overcome these problems. However, it is also known (e.g. Reese and Wriggers, 1996, 2000; Loehnert et. al. 2005; Jabareen and M.B. Rubin, 2007a, 2007b) that these improved formulations can exhibit unphysical hourglassing in regions experiencing combined high compression with bending. Furthermore, Jabareen and Rubin (2007a) showed that some of these improved element formulations in commercial codes can exhibit inelastic response even though they attempt to model a hyperelastic material with a strain energy function. Restrictions were developed on the strain energy function which ensure that the CPE reproduces exact solutions for all homogeneous deformations (Nadler and Rubin, 2003). Consequently, the CPE automatically satisfies a nonlinear form of the patch test. Also, a functional form of the strain energy for the CPE was proposed with specific dependence on the strain energy of the three-dimensional material and the reference geometry of the CPE element. In addition, a strain energy function for inhomogeneous deformations was introduced as a quadratic function of inhomogeneous strain measures.

In the original CPE formulation (Nadler and Rubin, 2003) the coefficients in the strain energy function for inhomogeneous deformations where determined by comparing linear solutions of a rectangular parallelepiped CPE with exact solutions of the linear three-dimensional equations for pure bending, pure torsion and higher-order hourglass modes of deformation. Loehnert et al. (2005) implemented the CPE formulation into the finite element code FEAP (Taylor, 2005) and considered specific example problems which showed that the response of the original CPE was robust and locking free for thin structures. However, it was also shown that the accuracy of the original CPE degraded with increased irregularity of the reference element shape. Recently, Boerner et al. (2007) have proposed a numerical method for determining coefficients in a quadratic form of the strain energy function for inhomogeneous deformations which improve the response of the CPE for 2-D plane strain problems with irregular elements. Jabareen and Rubin (2007c) developed analytical expressions for constitutive coefficients in an improved CPE for 3-D deformations by generalizing the quadratic form of the strain energy function for inhomogeneous deformations to include additional coupling of the inhomogeneous strains active in bending modes. Functional forms for the additional coefficients were determined by considering four bending solutions for special shaped reference elements which can be represented as parallelepipeds with two right angles. Jabareen and Rubin (2007b) also developed improved coefficients for torsion by considering an exact torsion-like solution of the linear elastic equations for an isotropic elastic material.

Thus, it would be advantageous to use the new CPE technology, which overcomes the unphysical hourglass instabilities that have been observed with elements based on enhanced strains or incompatible modes, to design tires, for example, to withstand the above mentioned large elastic deformations.

References

ABAQUS Inc., Version 6.5-1, Providence RI 02909-2499. ADINA R & D Inc., Version 8.3.1, Watertown, MA 02472.

ANSYS Inc., University Advanced Version 9, Canonsburg, PA 15317.

Belytschko T, Bindeman LP (1993) Assumed strain stabilization of the eight node hexahedral element. Comp. Meth. Appl. Mech. Engrg. 105: 225-260.

Belytschko T, Ong JSJ, Liu WK, Kennedy JM (1984) Hourglass control in linear and nonlinear problems. Comp. Meth. Appl Mech. Engrg. 43: 251-276. Boerner EFL Loehnert S, Wriggers P (2007) A new finite element based on the theory of a

Cosserat point - extension to initially distorted elements for 2D plane strain. Mt. J.

Numer. Meth. Engng. 71: 454-472.

Bonet J, Bhargava P (1995) A uniform deformation gradient hexahedron element with artificial hourglass control. International Journal of Numerical Methods in

Engineering 38: 2809-2828. Cesar de Sa JMA, Areias PMA, Jorge RMN (2001) Quadrilateral elements for the solution of elasto-plastic finite strain problems. International Journal of Numerical

Methods in Engineering 51: 883-917. Crisfield MA, Moita GF, Jelenic G, Lyons LPR (1995) Enhanced lower-order element formulation for large strains. Computational Mechanics 17: 62-73. Flory P (1961) Thermodynamic Relations For High Elastic Materials. Trans. Faraday Soc.

57: 829-838.

Hutter R, Hora P, Niederer P (2000) Total hourglass control for hyperelastic materials. Comp. Meth. Appl. Meek Engrg. 189: 991-1010.

Jabareen M, Rubin MB (2007a) Hyperelasticity and physical shear buckling of a block predicted by the Cosserat point element compared with inelasticity and hourglassing predicted by other element formulations. Computational Mechanics, 40: 447-459. Jabareen M, Rubin MB (2007b) Modified torsion coefficients for a 3-D brick Cosserat point element. Computational Mechanics 41: 517-525.

Jabareen M, Rubin MB (2007c) An improved 3-D Cosserat brick element for irregular shaped elements. Computational Mechanics 40: 979-1004. Jabareen M, Rubin MB (2008) A Cosserat point element (CPE) for nearly planar problems

(including thickness changes) in nonlinear elasticity. Submitted to International Journal for Engineering Science 46: 986-1010.

Loehnert S, Boerner EFL Rubin MB, Wriggers P (2005) Response of a nonlinear elastic general

Cosserat brick element in simulations typically exhibiting locking and hourglassing.

Computational Mechanics 36: 255-265.

Nadler B, and Rubin MB, (2003) A new 3-D finite element for nonlinear elasticity using the theory of a Cosserat point, Mt. J. Solids and Structures 40: 4585-4614. Reese S, Wriggers P (1996) Finite element calculation of the stability behaviour of hyperelastic solids with the enhanced strain methods. Zeitschrift fur angewandte Mathematik and Mechanik 76: 415-416.

Reese S, Wriggers P (2000) A stabilization technique to avoid hourglassing in finite elasticity. Int. J. Numer. Meth. Engng. 48:79-109.

Reese S, Wriggers P, Reddy BD (2000) A new locking free brick element technique for large deformation problems in elasticity. Computers and Structures 75: 291-304. Rubin MB (1995) Numerical solution of two- and three-dimensional thermomechanical problems using the theory of a Cosserat point, J. of Math, and Physics (ZAMP) 46, Special Issue, S308-S334. In Theoretical, Experimental, And Numerical Contributions To The

Mechanics Of Fluids And Solids, Edited by J. Casey and M. J. Crochet, Brikhauser Verlag,

Basel (1995). Rubin MB (2000) Cosserat Theories: Shells, Rods and Points. Solid Mechanics and its

Applications, Vol. 79, Kluwer, The Netherlands. Simo JC, Rifai MS (1990) A class of mixed assumed strain methods and the method of incompatible modes. Mt. J. Numer. Meth. Engng. 29: 1595-1638. Simo JC, Armero F (1992) Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. Int. J. Numer. Meth. Engng. 33: 1413-1449. Simo JC, Armero F, Taylor RL (1993) Improved versions of assumed enhanced strain tri- linear elements for 3D finite deformation problems. Comp. Meth. Appl. Mech. Engrg.

110: 359-386. Sokolnikoff IS (1956) Mathematical Theory of Elasticity. McGraw-Hill, New York. Taylor

RL (2005) FEAP - A Finite Element Analysis Program, Version 7.5, University of California,

Berkeley. Zienkiewicz OC, Taylor RL (2005) The Finite Element Method for Solid and Structural

Mechanics. Sixth edition, Elsevier.

SUMMARY OF THE INVENTION

Accordingly, it is a principal object of the present invention to provide a 3-D brick Cosserat Point Element (CPE), which is accurate, robust and user friendly. It is another principal object of the present invention to provide a CPE, which can be used with confidence for thin structures, such as shells and beams.

It is one other principal object of the present invention to provide a CPE that does not exhibit unphysical hourglass instabilities, which have been observed with elements based on enhanced strains or incompatible modes.

A method is disclosed for Cosserat Point Element (CPE) modeling of nonlinear elastic materials, utilizing a strain energy function characterizing the response of an eight node 3-D brick element to homogeneous and inhomogeneous deformations, which includes full coupling of bending and torsional modes of deformation.. The method includes determining the constitutive coefficients which satisfy bending and torsional modes of deformation for three specific reference element geometries and combining these solutions to obtain constitutive equations for general reference element shapes which ensure that the strain energy remains a positive definite function for all inhomogeneous deformations of the CPE. The constitutive coefficients characterize a generalized CPE model the response of a 3-D brick element for finite elastic deformations which is free of hourglass instabilities and locking for nearly incompressible materials, thereby providing a robust modeling tool.

The present invention provides a generalized form for the strain energy of inhomogeneous deformations, developed for a 3-D brick CPE, which includes full coupling of bending and torsional modes of deformation. The constitutive coefficients, which depend on the reference geometry of the element, are determined by solving eighteen bending problems and six torsion problems on special elements that are parallelepipeds with two right angles. The resulting constitutive coefficients ensure that the strain energy for inhomogeneous deformations remains a positive definite function of the inhomogeneous strain measures for all reference element shapes. A number of exemplary problems are considered which show that the generalized CPE produces results as accurate as enhanced strain and incompatible elements for thin structures and is free of hourglass instabilities typically predicted by these enhanced elements in regions experiencing combined high compression with bending.

Exact linear solutions of bending and simple torsion of these special elements are reexamined and it is shown that by including more general coupling of bending and torsional modes it is possible to match solutions of the two additional bending modes not handled in (Jabareen and Rubin, 2007c) as well as the simple torsion solutions. Matching these additional two bending modes significantly improves the response of the generalized CPE for problems of bending of a rhombic plate, which was shown to be inaccurate by Erhlich (2007).

Section 2 of the summary presents the basic equations of the CPE. Section 3 of the summary describes the procedure to use eighteen bending solutions and six torsions solutions to determine the functional forms of the constitutive coefficients in the strain energy function for inhomogeneous deformations. Section 4 of the detailed description, considers a number of illustrated examples that demonstrate the response of the generalized CPE. Section 5 presents conclusions.

2. Basic equations of a Cosserat Point Element (CPE)

Here attention is focused on an eight noded 3-D brick CPE as developed by Nadler and Rubin (2003) and modified by Jabareen and Rubin (2007b,c). Moreover, for the present discussion it is sufficient to confine attention to the equilibrium equations. Within the context of this theory the kinematics of the CPE are characterized by

where D₁ are constant reference directors associated with the reference configuration and d_j are present directors associated with the deformed present configuration. These element director vectors are related to the nodal director vectors [which locate the nodes of a brick element (see Fig. 1) relative to a fixed origin] by a constant matrix and details can be found in (Nadler and Rubin, 2003).

Without loss in generality, the three-dimensional position vector X associated with material points in the reference configuration of the element can be expressed in the form

where θ* (i=l,2,3) are convected material coordinates having the units of lengths, H₁ are constant lengths characterizing the element defined by

and N¹ are shape functions

Moreover, the region P₀ associated with the reference configuration is mapped to the region

P associated with the present configuration, which is bounded by SP which is characterized by the six surfaces SP_j (J=l,2,..,), such that

The element directors and nodal directors are related using a tri-linear form like (2.2) for the position vector x of material points in the present configuration. This causes the displacement field to be continuous across element boundaries. However, it is important to emphasize that the constitutive equations in the CPE approach are determined by a strain energy function for the structure and do not depend on the point-wise validity of this approximate expression for x*.

Now, in the absence of body force, the equilibrium forms of the balances of director momentum become

where m* are director couples due to surface tractions on the boundaries of the element, and t¹ are intrinsic director couples which require constitutive equations. In this regard, it is mentioned that director couples are kinetic quantities conjugate to the element directors and they do not necessarily have the units of mechanical moment. For the CPE it is convenient to introduce the deformation measures

where {D¹, d¹} (i= 1,2,3) are the reciprocal vectors of {D_i? dfi, respectively, and {V, V¹} characterize the reference configuration of the element and are defined in (Nadler and Rubin, 2003, Appendix B). Specifically, when the reference CPE is a rectangular parallelepiped then V represents the volume of the CPE and V¹ vanish. Nadler and Rubin (2003) developed restrictions on the strain energy of an elastic CPE which ensure that the element satisfies a nonlinear form of the patch test for all uniform homogeneous nonlinear elastic anisotropic materials. In particular, the specific (per unit mass) strain energy Σ of the CPE can be expressed in terms of the specific strain energy Σ* of the three-dimensional material and the specific strain energy Ψ associated with inhomogeneous deformations, such that

where the kinematic quantities are defined by

Here, it is convenient to introduce the alternative variables b_j (i=l,2,...,12) by

Then, following the work in Jabareen and Rubin (2007b,c) the strain energy Ψ for inhomogeneous deformations is expressed in the form

where m is the mass of the element, {μ, v} are the shear modulus and Poisson's ratio associated with the small deformation response, B_ij is a constant symmetric matrix with the response to higher order hourglassing being uncoupled from that to bending and torsion so that B_ij has the following zero components

and the values of

are given by the expressions associated with higher order hourglassing in (Nadler and Rubin, 2003). It then can be shown (Nadler and Rubin, 2003) that the constitutive equations for a hyperelastic CPE become

and the remaining quantities tⁱ are determined by

These constitutive equations automatically satisfy a nonlinear form of the patch test (Nadler and Rubin, 2003). Next, for elastically isotropic materials it is convenient to use the work of Flory (1961) to

introduce the dilation

, the pure measures of distortion and the scalar measures of

distortion

Then, the strain energy function Σ* can be written in the form

•

so that

Also, for the example problems considered later is specified in terms of the small deformation bulk modulus K and shear modulus μ, such that

where P_Q is the constant reference density of the material, the mass m of the CPE is given by

and K and Young's modulus E associated with the small deformation response satisfy the equations

Furthermore, for later reference it can be shown (Nadler and Rubin, 2003) that the director couples m* are related to the traction vector t^* applied to the boundary d? by the integrals

where da ^f is the current element of area. Also, it can be shown (Loehnert et al. 2005) that d^1/2T is equal to the volume integral of the three-dimensional Cauchy stress T^*

where dv^* is the current element of volume.

3. Determination of the constitutive coefficients

In this paper the constitutive coefficients By in (2.11) are determined by matching solutions of the linearized equations for the CPE with exact solutions of the linear theory of elasticity for special element shapes. Specifically, the classical pure bending solution (e.g. Sokolnikoff, 1956) of the three dimensional equations of elasticity for a rectangular parallelepiped can be written in the form

In this solution, X^* locates a material point in the reference configuration, u^* is the displacement vector, T^* is the stress tensor, the constant γ controls the magnitude of the bending field and the components of the tensors are referred to the right-handed orthonormal base vectors

. Also, the simple torsion-like solution in (Jabareen and Rubin,

2007b) can be expressed in the form

where the constant ω is the twist per unit length in the ej direction and the constant φ controls the warping of the cross-section with unit normal e-[.

In order to determine values of the coefficients EL- it is convenient to consider the following three elements shapes, which are defined in terms of another fixed rectangular Cartesian triad β_j by

Element

Element

Element

where the metric D- is defined by

Each of these element shapes is a parallelepiped with two right angles. Now, from (Nadler and Rubin, 2003) it follows that for these element shapes the reference geometry is characterized by

and the position vector X^H in (2.2) and the Cartesian coordinates X_j' can be expressed in the forms

Consequently, the exact displacements u^* and stresses T^* in (3.1) and (3.2) can be rewritten as functions of the convected coordinates θ¹ (i=l,2,3).

Within the context of the linear theory of a CPE (Nadler and Rubin, 2003) the director displacements δ_{ are defined such that

and for the special elements defined by (3.3) the linearized forms of the inhomogeneous strains β_j become

As explained in (Nadler and Rubin, 2003), the values δ^* of the element director displacements δ_j which correspond to the exact displacement field u* need to be properly defined. Specifically, for these element shapes the values

are determined by the equations in (Nadler and Rubin, 2003) which connect δ \ to integrals over the reference element region of derivatives of u* with respect to the convected coordinates. In particular, for the exact solutions (3.1) and (3.2) and the element shapes (3.3) it can be shown that these expressions yield

so that when δ_j are replaced by the exact values δ \ the linearized values of κ\ vanish

and the linearized forms of the constitutive equations (2.13) and (2.14) reduce to

where d¹ have been replaced by the reference values D*. Also, the values of m* in (2.21) associated with the exact solutions (3.1) and (3.2) are given by

where 9P₀ is the reference boundary of the CPE, N^* is the unit outward normal to and

dA^ is the reference element of area. It then follows that within the context of the linearized theory, the equations of equilibrium associated with the bending (3.1) and torsion (3.2) solutions reduce to three vector equations

Analytical expressions for B₁- can be developed by matching the solutions (3.1) and (3.2) for each of the element shapes (3.3). Then, the resulting coefficients are combined in a manner that ensures that B-- is a positive definite tensor. Specifically, with reference to the element shape El in (3.3a) consider six bending solutions associated with specifications of the orientation of e- relative to D₁

Also, consider two torsion solutions associated with specifications

For each bending and torsion solution the exact displacements u^* are used to determine the exact values δ-. Then, the linearized values of κ| are determined using (2.9) and (3.8) with δ_j replaced by δ_j and the resulting constitutive equations for t¹ are determined by (3.11).

Similarly, the values of m* are determined by using the exact stress T^* in the equations (3.12). Also, the values of the warping constant φ which corresponds to nearly pure torsion (Jabareen and Rubin, 2007b) are determined by

In these expressions it can be seen that {m⁶, D₁) are associated with the cross-sectional coordinates and {m⁵, D₂) are associated with the cross-sectional coordinates

There has thus been outlined, rather broadly, the more important features of the invention in order that the detailed description thereof that follows hereinafter may be better understood. Additional details and advantages of the invention will be set forth in the detailed description, and in part will be appreciated from the description, or may be learned by practice of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

For a better understanding of the invention with regard to the embodiments thereof, reference is now made to the accompanying drawings, in which like numerals designate corresponding elements or sections throughout, and in which: Fig. 1 is a prior art sketch of a general brick CPE showing the numbering of the nodes and the surfaces;

Fig. 2 is a schematic illustration of the cross-section of a parallelepiped element El, constructed according to the principles of the present invention; Fig. 3 is a schematic illustration of a shear load on a thin cantilever beam, where the irregular element mesh is based on the distorted center cross-section, constructed according to the principles of the present invention;

Figs. 4a-4d show the response of a shear load on a thin cantilever beam for small deformations, constructed in accordance with the principles pf the present invention; Fig. 5 shows a sketch of a thin slanted cantilever beam with dimensions (4.4) and with the slanting angle θ, constructed in accordance with the principles of the present invention; Fig. 6a shows the displacement component U^₃ of point A in the e₃ direction as a function of θ for the most refined mesh (n=5), and is plotted in accordance with the principles of the present invention; Figs. 6b and 6c show that (G) and (F) converge to the same values and that (G) is slightly more accurate than (F) for n=l and large values of θ;

Fig. 7 shows a shear load on a thin twisted cantilever beam where the element mesh is based on the distorted center cross-section, and is constructed in accordance with the principles of the present invention; Figs. 8a - 8e are graphs for shear loading of a thin twisted cantilever beam with a pre- twist θ and the mesh {20nx5nχn} (small deformations), plotted in accordance with the principles of the present invention;

Figs. 9a - 9d are graphs for the results of large deformation lateral torsional buckling of a thin cantilever beam with a small pre-twist of θ = 0.1 deg. using the mesh {20n^χn^χn}, plotted according to the principles of the present invention;

Figs. 10a - 1Od are graphs for large deformation shear loading of a thin cantilever beam with a pre-twist Θ = 30 deg. using the mesh {20nxnxn}, plotted according to the principles of the present invention;

Fig. 11 is a schematic illustration of a point load on the corner of a thin partially clamped rhombic plate, constructed according to the principles of the present invention; Figs. 12a - 12c are graphs of the response to a point load on the corner of a thin partially clamped rhombic plate for small deformations, plotted according to the principles of the present invention;

Fig. 13 is a sketch of a point load on the center of a thin fully clamped square plate with an irregular element mesh, constructed according to the principles of the present invention;

Figs 14a and b are graphs of the response to a point load on the center of a thin fully clamped square plate with small deformations, plotted according to the principles of the present invention;

Figs. 15a and 15b show the deformed shapes of a thin partially clamped rhombic plate subjected to a point load on its corner for two different angles θ = 0 and 60° and the same load value P, constructed according to the principles of the present invention;

Figs. 16a - 16d are graphs of the response to a point load on a partially clamped rhombic plate for large deformations, plotted according to the principles of the present invention;

Fig. 17 shows a sketch of one eighth of a thin circular cylindrical shell that is subjected to a pair of opposing point loads P, constructed according to the principles of the present invention;

Figs. 18a and 18b are graphs of the force P versus radial displacement u_r of the point A under the load for different mesh refinements, plotted according to the principles of the present inventions showing that {G, F} tend to converge to the same solutions; Fig. 19 shows the deformed shape of one eighth of the circular shell predicted by (G) for n=5 and P = 1.73 kN, with no enhancement of the displacements, constructed according to the principles of the present invention;

Fig. 20 is a schematic illustration of plane strain indentation of a rigid plate into a nearly incompressible block showing the boundary conditions and definition of element irregularity for this problem, constructed according to the principles of the present invention;

Figs. 21a and 21b are graphs of the response of plane strain indentation of a rigid plate into a nearly incompressible block, plotted according the principles of the present invention;

Figs. 22a and 22b are graphs representing plane strain indentation of a rigid plate into a nearly incompressible block, and showing the error E in the displacement u_C2 of the point C for two cases of element irregularity and for the mesh {8nχ4nxl} with n=5 and ^UA₂ = - 0.1 m, plotted according to the principles of the present invention;

Figs. 23a - 23e represent plane strain indentation of a rigid plate into a nearly incompressible block showing graphs of nonlinear load curves using the regular mesh {8nx4nxl} for different values of n, plotted according to the principles of the present invention;

Figs. 24a - 24f show deformed shapes for plane strain indentation of a rigid plate into a nearly incompressible block for the regular mesh {8nx4nxl} with n = 5, constructed according to the principles of the present invention; Fig. 25 is a schematic representation of indentation of a rigid plate into a block showing the boundary conditions and definition of element irregularity, where symmetry conditions are used so that only one fourth of the block is meshed, and constructed according to the principles of the present invention;

Fig. 26 is a graphical representation of the response to indentation of a rigid plate into a nearly incompressible block, where the displacement u_E3 of the point E as a function of element irregularity for the mesh {4nx4nx4n} with n = 3 and u^ = - 0.1 m is plotted according to the principles of the present invention; and

Figs. 27a - 27d are graphic representations of the response to indentation of a rigid plate into a nearly incompressible block, plotted according to the principles of the present invention.

DETAILED DESCRIPTION OF AN EXEMPLARY EMBODIMENT

The principles and operation of a method and an apparatus according to the present invention may be better understood with reference to the drawings and the accompanying description, it being understood that these drawings are given for illustrative purposes only and are not meant to be limiting.

With reference now to Fig. 2, a schematic illustration of the cross-section of a parallelepiped element El 200, constructed according to the principles of the present invention, it can be seen that the torsion expressions (3.16a) and (3.16b) are similar to those in (Nadler and Rubin, 2003), since the cross-section normal to D₁ 210 has lengths and the cross-section normal to D₂ 220 has lengths .

For each bending solution the value of γ can be eliminated in the resulting equations of equilibrium (3.13) and the value of ω can be eliminated from each of the equations of equilibrium associated with the torsion solutions. Also, the values (3.16) are used in the resulting torsion equations. It therefore follows that each of the solutions (B1)-(B6), (Tl) and (T2) yield nine scalar equations of equilibrium which total 72 scalar equations to determine the values of By as functions of H_j and D₁₂. Some of these scalar equations are trivially satisfied and others are redundant. In particular, using a symbolic program like Maple it can be shown that the equations associated with the bending solutions (3.14) can be solved for By (i, j=l, 2,..5) in terms of

Then, these expressions can be substituted into the equations for torsion to determine all of the values of B- (Lj= 1,2,...,9) in terms of one of the values (B₇₇, B₈₈, B₉₉}, say B₇₇.

Next, it is noted that that the strain energy function Ψ in (2.11) will be positive definite provided that the coefficient matrix By is positive definite. The results of the calculations just described yield a form for By which separates into two matrices [By (i,j=l,2,3,4) and By (i,j=5,6,..9)]. In particular, it can be shown that

which is independent of the lengths H₁ and is positive for the full range of Poisson's ratio and D₁₂

Moreover, the det(By) for (i,j=5,6,...,9) is a linear function OfB₇₇ which has the value

It was found that if the value of B₇₇ is specified so that

then the expressions for B- are quite simple and

which is positive for the range of values (3.19).

Similar procedures can be used to define bending and torsion solutions for the element shapes E2 and E3 and the resulting equations can be solved for B₁- to determine the dependence on the metrics D₁₃ and D₂₃. Next, introducing the auxiliary variables {λ₁₂, λ₁₃, λ₂₃} defined by

it is possible to denote the values of B- associated with the solutions of the three elements

E1-E3 in (3.3) by for El, by for E2, and by for E3. Moreover, the matrix BJj is

defined so that it yields a strain energy function Ψ equivalent to that obtained in (Nadler and

Rubin, 2003) for a rectangular parallelepiped, when the value of the torsion function b^*(l) is taken to be 1/2 as suggested in (Jabareen and Rubin, 2007b). Then, the general expression which combines these solutions is given by

Now, using the definitions (3.23) it follows that each of the coefficients is non-negative and that at least one of them is positive.

Also, each of the matrices , given in Appendix A, is positive definite so

that the combined matrix is also positive definite for all reference element

shapes.

4.0 Examples For planar problems (in the D₁-D₂ plane with D₁₃=D₂₃=O) it was shown in (Jabareen and Rubin, 2008) that b₅ = b₆ = b₇ = b₈ = b₉ = b₁₀ = b_n = 0 . (4.1)

Moreover, since the values of By (i, j = 1, 2, 3, 4) for these planar problems reduce to those of the improved CPE in (Jabareen and Rubin, 2007c), the results of the improved and generalized CPE will be identical for all planar problems. Consequently, the results of the example problems of a Cook's membrane, a Kirsch problem and buckling of a block considered in (Jabareen and Rubin, 2007a, c) are identical to those that would be predicted by the generalized CPE of Sections 2 and 3. In particular, it follows that the generalized CPE is free of the unphysical hourglassing that is predicted by the elements in ABAQUS, ADINA, ANSYS and FEAP (Taylor, 2005) based on enhanced strain and incompatible mode methods.

In this section a number of example problems are considered to examine the predictions of the generalized CPE which show that it is more accurate than the improved CPE in (Jabareen and Rubin, 2007c). For all of the example problems the full nonlinear equations are solved using the constitutive equations (2.13) and (2.14) with strain energy specified by (2.8), (2.13), (2.14) and (2.18), even when the loads are small and the deformations remain small. Also, for irregular shaped elements the values of V¹ in (2.7) and (2.13) can be nonzero so the response is examined for more general conditions than those used in the last section to develop expressions for the constitutive coefficients By. Unless otherwise stated the material is taken to be compressible with K = I GPa, μ = 0.6 GPa, v = 0.25 . (4.2)

For the special examples which consider a nearly incompressible material the material constants are specified by

K = 1000 GPa , μ = 0.6 GPa , v « 0.4997 . (4.3)

Furthermore, it was shown in (Jabareen and Rubin, 2007a, c) that the enhanced strain element in FEAP produces results similar to the enhanced strain or incompatible mode elements in ABAQUS, ADINA and ANSYS so that comparison with these types of elements will be limited to the element in FEAP.

For the calculations with the improved CPE presented here use is made of the modified torsion coefficients discussed in (Jabareen and Rubin, 2007b). In addition, for the improved CPE the matrices By associated with the individual metrics {D₁₂, D₁₃, D₂₃) are combined using an expression of the form (3.24) instead of using the method described in (Jabareen and Rubin, 2007c) for ensuring that B^ remains positive definite.

Furthermore, in the following figures the symbols {G, I, F, Ql, QlPO, HO9} denote predictions of: the generalized CPE developed here; the improved CPE developed in (Jabareen and Rubin, 2007c); the enhanced strain, full integration, mixed element for nearly incompressible materials and the mixed higher order 9 node quadrilateral in FEAP, respectively. Also, it is noted that when the elements are rectangular parallelepipeds with (D₁₂=D₁₃=D₂₃=O) the predictions of (G) and (I) are identical.

4.1 Shear load on a thin cantilever beam for small deformations

Fig. 3 is a schematic illustration of the shear load on a thin cantilever beam, where the irregular element mesh is based on the distorted center cross-section, constructed according to the principles of the present invention. Fig. 3 shows a sketch of a thin cantilever beam with dimensions L 311 = 200 mm , H 312 = W 313 = 10 mm , (4.4) in the C₁ 310, e₂ 320, e₃ 330 directions, respectively, which is fully clamped at its end X₁ = 0, and is subjected to a shear force P 300 (modeled by a uniform shear stress) applied in the e₂ direction to its end X₁ = L. The lateral surfaces are traction free. The mesh {20nxnxn} is defined by distorting the middle cross-section in its reference configuration (using the parameters a_1} a₂, a₃, a₄ shown in Fig. 3), with 1On elements on each side of this cross-section and n elements in each of the e₂ 320 and e₃ 330 directions. Two cases of element distortion are considered

and the parameter a/H defines the element irregularity. Both of these cases cause the middle surface to remain planar with the normal to that surface being in the C₁-C₂ P^{lane for Case I} and in the C₁-C₃ plane for Case II. The value

of the e₂ component of the displacement of point A (see Fig. 3) predicted by (G) with the most refined mesh (n = 5) and a regular mesh (a/H=0) is considered to be exact and the error E associated with the predictions u^ of other calculations for the same value of P is defined by

Figs. 4a - 4d show the results for a shear load on a thin cantilever beam for small deformations, plotted in accordance with the principles pf the present invention. Figs. 4a and 4c show errors in the displacement of the point A in the e₂ direction versus variation in the distortion parameter a/H 410 and Figs. 4b and 4d show the errors versus n for the mesh {20nxnxn} defined for two cases of element distortion.

Fig. 4a shows the error as a function of the irregularity parameter a/H for Case I with n = 1. Ideally the response should be nearly insensitive to the value of a/H. This figure shows that the predictions of (G) and (I) are identical and are slightly more accurate than those of (F) for large values of irregularity. Fig. 4b shows that the three elements converge to the same value for Case I with the refined mesh (n=5) and large irregularity a/H = 2 420. The results for Case II shown in Figs. 4c and 4d indicate that (G) is again slightly more accurate than (F), and that they both converge to the same solution. In contrast, Figs. 4c and 4d for fixed a/H = 2 show that (I) predicts significant errors and converges slowly for large irregularities. It will be shown later that this deficiency for Case II causes (I) to be inaccurate for out-of -plane bending of a rhombic plate.

4.2 Shear load on a thin slanted cantilever beam for small deformations

Fig. 5 shows a sketch of a thin slanted cantilever beam 500 with dimensions (4.4) and with the slanting angle θ 510, constructed in accordance with the principles of the present invention. The boundary conditions are the same as those for the previous example except that the shear load P 520 is applied in the e₃ 530 direction to emphasize differences between the predictions of (G) and (I). Again, the mesh is taken to be {20nxnxn} with 2On elements in the axial direction of the beam. All of the elements have parallelogram cross-sections in the C₁-C₂ plane, with sides parallel to the ends of the beam.

Fig. 6a shows the displacement component u_A3 610 of point A (see Fig. 5) in the e₃ direction as a function of θ 620 for the most refined mesh (n = 5), and is plotted in accordance with the principles of the present invention. The error E 630 in u^ is defined in a similar manner to (4.7) with the exact value u^ taken to be that predicted by (G) for each value of θ 620 with n = 5 and with the load P given by (4.6). Figs 6b and 6c show that (G) and (F) converge to the same values and that (G) is slightly more accurate than (F) for n = 1 and large values of θ 620. Also, it can be seen that (I) predicts significant errors with slow convergence.

4.3 Shear load on a thin twisted cantilever beam for small deformations

Fig. 7 shows a shear load on a thin twisted cantilever beam where the element mesh is based on the distorted center cross-section, and is plotted in accordance with the principles of the present invention. The problem of shear loading of a thin twisted beam provides a severe test of the accuracy of an element formulation since the elements have irregular shapes and the response couples torsion and bending modes of deformation. In its unstressed reference configuration the twisted beam has length L 711 and a rectangular cross-section with height H 712 and width W 713 given by

L=200mm, H=10mm, W=2mm (4.8)

Also, each of the cross-sections is twisted by the angle θ such that the position of a material point X^* in the reference configuration is given by

where Θ controls the magnitude of the twist. Furthermore, the influence of irregular element meshes is explored by applying twist θ 710 to the cross-sections of the beam shown in Fig. 3 (and alternatively in Fig. 7) using the element irregularity specified by Case I in (4.5). The surface X^'=0 is fully clamped, a shear force P 720 (modeled by a uniform shear stress) is applied to the end X^=L in the constant e₂ direction (defined by its direction in the reference configuration) and the remaining lateral surfaces are traction free. Also, the element mesh for this problem is specified by {20nxnxn} with 2On elements along the length of the beam.

Figs. 8a - 8e are graphs of the results for shear loading of a thin twisted cantilever beam with a pre-twist Θ and the mesh {20nx5nxn} (small deformations), plotted in accordance with the principles of the present invention. The influence of element irregularity is shown in

Fig. 8a 801 for n = 5, Figs. 8b 802 and 8c 803 for n = 1 and in Figs. 8d 804 and 8e 805 for n = 2. Figure 8a shows the component Ug₂ of the displacement of the point B (see Fig. 7) in the e₂ direction

versus the pre-twist Θ for the most refined mesh with n = 5, no element irregularity (a/H=0) and the load given by P = 0.01 N . . (4.11)

This value of load is used for all calculations in this subsection. It can be seen that {G, I, F} converge to the same results. Thus, the value u_B*₂ of u_B2 predicted by (G) for n = 5 is considered to be exact and the error E of other calculations is defined in a similar manner to (4.7). Figs. 8b, c show the influence of element irregularity for the coarse mesh with n = 1 and Figs. 8d, e consider a finer mesh with n = 2. From these figures it can be seen that (G) is as accurate as (F) and is much more accurate than (I), especially for irregular shaped elements.

4.4 Lateral torsional buckling of a thin cantilever beam for large deformations For this example the beam and element mesh are characterized by the same parameters as used in the previous sub-section (4.3) with the shear force P being applied in the e₂ direction.

To investigate rotation of the beam's end it is convenient to consider the difference in the displacements of the points A and B shown in Fig. 7. Specifically, the quantity Δu is defined by

Figs. 9a - 9d 901-904, are graphs for results of large deformation lateral torsional buckling of a thin cantilever beam with small pre-twist of Θ = 0.1° using the mesh {20nxnxn}, plotted according to the principles of the present invention. The influence of element irregularity is shown in both Figs. 9a and 9b for n = 2 and a/H = 0 and 2 and in Figs. 9c, d for n = 3. Again it is emphasized that the direction e^ of the load is held constant during loading. The buckling process is triggered by the small pre-twist

Θ = 0.1 deg , (4.13)

which smoothes out the bifurcation that would occur for a perfect beam with Θ = 0.

Figs. 10a - 1Od 1001 - 1004 are graphs of the predictions for large deformation shear loading of a thin cantilever beam with a pre-twist for Θ = 30 deg. using the mesh {20nxnxn}, plotted according to the principles of the present invention. Fig, 10 shows the results for shear loading with a large pre-twist Θ = 30° . (4.14)

Moreover, the curves in Figs. 9 and 10 denoted by (E) are predicted by (G) with n = 5 and a/H = 0 and are considered to be exact.

The results in Fig. 9 show that for n = 2 (Figs. 9a and 9b) the predictions are not yet converged and are sensitive to element irregularity with large errors being predicted by (I). Figs 9c and 9d show that for n = 3 the predictions are reasonably converged and that the sensitivity of element irregularity is reduced for {G, F} but that (I) still predicts large errors. The results in Fig. 10 for shearing of a thin cantilever beam with a large pre-twist again show that the predictions of {G, F} are relatively accurate but that the predictions of (I) are inaccurate for the irregular shaped elements even for the mesh with n = 3 (Fig. 1Od). This result is consistent with that in Jabareen and Rubin, 2007c, Fig. 8b, which showed that the error of out-of-plane bending of a beam with the element irregularity of Case I [(4.5) here] does not have a zero slope as a/H approaches zero. Furthermore, this error causes (I) to predict inaccurate results for out-of-plane bending of a rhombic plate, as will be shown in the next sub-section. 4.5 Point load on corner of a thin partially clamped rhombic plate for small deformations

Fig. 11 is a schematic illustration a point load on the corner of a thin partially clamped rhombic plate, constructed according to the principles of the present invention. Dimensions are L = 500 mm 1101 , H = 10 mm 1102 , (4.15)

The value L corresponds to the actual length of each edge and the value of P is given by

P 1103 = 1 N . (4.16)

The mesh used for the plate is defined by {IOnxiOnxn} with n elements through the thickness.

Figs. 12a - 12c are graphs of the response to a point load on the corner of a thin partially clamped rhombic plate for small deformations, plotted according to the principles of the present invention. Figs. 12a - 12c are graphs showing the component u^ of the displacement of the point A in the e₃ direction as a function of θ for the most refined mesh n = 5. The error E in this displacement is defined in a similar manner to (4.7) with the exact value U^₃ taken to be that predicted by (G) for each value of θ with n = 5 and the load P given by (4.16). Figs. 12b and 12c show that {G, F} predict nearly the same values, that (I) predicts significant errors for n = 1, especially for the angle θ = 45°, and that (I) exhibits slow convergence.

4.6 Point load on the center of a thin fully clamped square plate with an irregular element mesh for small deformations

Fig. 13 is a sketch of one quarter of a thin fully clamped square plate that is loaded by a point force P 1310 at its center, constructed according to the principles of the present invention. Dimensions are given in (4.8). Only one quarter of the plate is modeled and the value of force P 1310 given by (4.16) corresponds to one quarter of the load that would be applied to the center of the entire plate. Irregular elements are specified by moving the center point of the quarter section to the position characterized by the lengths {a_1? a₂}, as shown in Fig. 13, and defined by the two cases

The quarter section of the plate is meshed by {IOnxiOnxn} with each subsection being meshed by {5nx5nxn} and with n elements through the thickness. The error E in the displacement component u_A3 of point A 1320 in the e₃ 1330 direction is defined in a similar manner to (4.7) with the exact value U^₃ taken to be that predicted by (G) for regular elements (a/L=0) with n = 5

Figs. 14a and b are graphs of the point load on the center of a thin fully clamped square plate with small deformations, plotted according to the principles of the present invention. Errors in the displacement of the point A in the e₃ direction versus the distortion parameters

(a) 4a/L and; (b) the angle θ for two cases of element irregularity with the mesh {10x10x1}. Figs 14a and b show the error for n=l as a function of the irregularity parameters 4a/L for Case 1 1410 in Fig. 14a and as a function of θ/(2π) for Case

II 1420 in Fig. 14b, respectively. From these Figs. 14a and 14b it can be seen that {G, F} are relatively insensitive to the magnitude and type of element irregularity but that (I) predicts significant errors for irregular elements.

4.7 Point load on the corner of a thin partially clamped rhombic plate for large deformations Figs. 15a and 15b show the deformed shapes of a thin partially clamped rhombic plate subjected to a point load on its corner for two different angles θ = 0 1510 and 60° 1520 and the same load value P 1520, constructed according to the principles of the present invention.

The plate is fully clamped on two edges and the other edges and major surfaces are traction free. The dimensions are given by (4.8) as shown in Fig. 11, with L now being the length of the plate's edge, and the point force P is

P = 1 kN . (4.19)

The mesh is specified by {IOnxiOnxl} and the exact value U₃ ^* of the displacement of the corner in the e₃ direction is determined by the most refined solution (G) with n=5 ujj = 0.21084 m for θ = 0 deg , U₃ = 0.39306 m for θ = 60 deg . (4.20)

Figs. 16a - 16d are graphs of the response to a point load on a partially clamped rhombic plate for large deformations, plotted according to the principles of the present invention. Figs. 16a and 16c show the load P 1610 versus displacement curves for n = 2 and the convergence curves for two values of the angle θ. Comparison of Figs. 16a and 16c shows that the rhombic plate with angle θ = 60 deg is more flexible than that for θ = 0 deg, and that {G, F} predict nearly the same values, whereas (I) predicts significant errors for the angle θ = 60 deg. Also, Fig. 16d shows that the convergence properties of (G) are slightly better than those of (F) for the case when θ = 60 deg.

4.8 A pair of opposing point loads applied to a complete circular cylindrical shell for large deformations

Fig. 17 shows a sketch of one eighth of a thin circular cylindrical shell that is subjected to a pair of opposing point loads P 1710, constructed according to the principles of the present invention. The entire shell has length 2L, middle surface radius R 1720, and thickness H

1730

L = 300 mm , R = 300 mm , H = 3 mm . (4.21)

These are exemplary measurements, as are all others for all other drawings, herein. All nodes (except for one) on the circular edges of the shell are allowed to move freely in the axial direction but their radial and circumferential positions are fixed. Also, the eighth region of the shell is modeled by the mesh {IOnxiOnxl} in the axial, circumferential and radial directions, respectively. Figs. 18a and 18b are graphs of the force P versus radial displacement u_r of the point A under the load for different mesh refinements, plotted according to the principles of the present invention. Curves are presented for (G) for the most refined mesh of n = 10 and for {G, F} for n = 3 1810 in Fig. 18a and n = 5 1820 in Fig. 18b. It can be seen from Fig. 18a that for the coarser mesh (n = 3) the load-deflection curve exhibits ratcheting due to localized limit points whereas for the more refined mesh (n = 5) in Fig. 18b the load-deflection curve is smooth. It can also be seen that the predictions of {G, F} tend to converge to the same solutions.

Fig. 19 shows the deformed shape of one eighth of the circular shell 1910 predicted by

(G) for n = 5 and P = 1.73 kN 1920, with no enhancement of the displacements, constructed according to the principles of the present invention. In particular, it is noted that the inability of the coarse mesh to capture the high curvature of the middle of the shell far away from the load is most likely the cause of the ratcheting shown in Fig. 18a.

4.9 Plane strain indentation of a rigid plate into a nearly incompressible block for large deformations

Crisfield, et al, (1995) and Cesar de Sa, et al, (2001) considered the example of plane strain indentation of a rigid plate into a block and showed limitations of enhanced strain elements for elastic and elastic-plastic response.

Fig. 20 is a schematic illustration of plane strain indentation of a rigid plate into a nearly incompressible block showing the boundary conditions and definition of element irregularity for this problem, constructed according to the principles of the present invention. The block has length 2L, height L and depth W. Material points on its sides and bottom remain in contact with a rigid container and are allowed to slide freely. The top surface of the block is loaded by a rigid plate (AB) of length L which makes perfect contact with the block so that material points in contact with the rigid plate move only vertically. The remaining half of the block's top surface is traction free. The dimensions of the block are given by L = W = I m (4.22) Irregular meshes are defined by dividing the block into four subsections, with the central node moving to the position characterized by the lengths {a_χ, a₂}, as shown in Fig. 20, and as defined by two cases

The entire block is meshed by {8nx4nxl} with 4n elements in the C₁ direction 2010 and 2n elements in the e₂ direction 2020 in each of the subsections. The point C 2030 is located on the free top surface at a distance 0.25 L from the corner B 2040 of the rigid plate. Also, the material is considered to be nearly incompressible.

Figs. 21a and 21b are graphs of the response of plane strain indentation of a rigid plate into a nearly incompressible block, plotted according the principles of the present invention. Convergence of the error E in the displacement u_C2 of the point C using the regular (a = 0) mesh {8nx4nxl} for u^ = — 0.1 m versus: n in Fig. 21a and versus the number of degrees of freedom (DOF) in Fig. 21b. The converged value u£₂ of the displacement of the point C in the e₂ direction predicted by (G) for a regular mesh with n = 20 is considered to be exact and is given by uj₂ = 0.071895 m for U^ = - 0.1 m with n = 20 . (4.24)

The error of E in the values u_C2 predicted by calculations of other elements and meshes is defined by an expression similar to (4.7). Fig. 21 shows the convergence of this error predicted by {G, QlPO, HO9}. This error is plotted relative to n for the mesh {8n^χ4n^χl} in Fig. 21a and is plotted relative to the degrees of freedom (DOF, calculated for plane strain response) in Fig. 21b. From Fig. 21a it is not clear if (QlPO) exhibits a locking behavior by converging to a value different from (G) or whether the convergence rate is very slow. To validate the converged value of (G) for n=20, calculations were also performed using the mixed higher order element (HO9) with the mesh {8nx4nxl} up to n = 10. In particular, it can be seen in Fig. 21b that (HO9) tends to converge to the value predicted by (G). Figs. 22a and 22b are graphs representing plane strain indentation of a rigid plate into a nearly incompressible block, and showing the error E in the displacement u_C2 of the point C for two cases of element irregularity and for the mesh {8nx4nxl} with n = 5 and u^ = - 0.1 m, plotted according to the principles of the present invention. Since there is a strain concentration near the edge of the plate it is expected that a non-fully converged solution will be sensitive to element irregularity. In particular, it can be seen from Fig. 22a that

(QlPO) is more sensitive to element irregularity than (G) for positive values of Case I 2210.

This causes the elements near the plate's edge B to be more irregular. The results in Fig. 22b for case II 2220 show that the error reduces slightly for increasing positive values of a, and cause the elements near the plate's edge B to be more refined.

Figs. 23a — 23e represent plane strain indentation of a rigid plate into a nearly incompressible block showing graphs of nonlinear load curves using the regular mesh {8nχ4nχl} for different values of n, plotted according to the principles of the present invention. Fig. 23a shows nonlinear load curves of (G) using the regular mesh {8n^χ4n^χl} for different values of n. Fig. 23b shows nonlinear load curves of (QlPO) using the regular mesh {8nx4nxl} for different values of n. Figs. 23c, 23d, and 23e show nonlinear load curves of (G) and (QlPO) using the regular mesh {8nx4nχl} for n = 3, 5 and 10, respectively. These figures again show that (G) predicts more flexible response than (QlPO) for the coarser meshes.

Figs. 24a - 24f show deformed shapes for plane strain indentation of a rigid plate into a nearly incompressible block for the regular mesh {8n^χ4nxl} with n = 5, constructed according to the principles of the present invention. The left column shows the results for (G) and the right column shows the results (QlPO), each for 3 different displacements: u^ = -0.15 (2410 and 2420), -0.20 (2430 and 2440) and -0.25 (2450 and 2460), respectively. In particular, it can be seen that the flexibility of (G) allows the elements near the plate's corner to roll around the corner more easily than allowed by (QlPO). Since the flexibility of (G) has been validated relative to the mixed higher order element (H09) it is concluded that the stiffness shown by (QlPO) is unphysical. 4.10 Indentation of a rigid plate into a nearly incompressible block for large deformations

Fig. 25 is a schematic representation of indentation of a rigid plate into a block showing the boundary conditions and definition of element irregularity, where symmetry conditions are used so that only one fourth of the block is meshed, and constructed according to the principles of the present invention. Fig. 25 shows a sketch of one fourth of a nearly incompressible block that has total length 2L, height L and depth 2L with

L = I m . (4.25)

The bottom (X₃=O) and exterior lateral surfaces

of the block remain in contact with and slide freely on smooth rigid planes. The block's top surface (X_j=L) is loaded by a rigid plate (ABCD) 2510, which makes perfect contact with the material points so that these points can only move vertically in the e₃ direction 2520. The remaining portion of the block's top surface is traction free. Irregular elements are generated by moving the nodes of the center plane by the displacements {a_l5 a^, a₃, a₄} as shown in Fig. 25 with

The mesh for the one fourth region is specified by {4nx4n^χ4n} with 2n elements below and above the distorted center surface and with {nxn} elements under the rigid plate. Moreover, the vertical reference locations of material points on these distorted surfaces are described by a bi-linear form of the coordinates (X^, X^). Furthermore, the point E 2530 is located at a distance L/4 from the edge of the rigid plate.

Based on the results of the previous example it is expected that the deformation will be concentrated near the edges of the plate so that a refined mesh will be required to obtain an accurate solution. Mesh refinement of this 3-D problem is beyond the capacity of the hardware being used to obtain the solution. Therefore, attention will be focused on the robustness of the solutions with relative coarse meshes.

Fig. 26 is a graphical representation of the response to indentation of a rigid plate into a nearly incompressible block, where the displacement u_E3 of the point E as a function of element irregularity for the mesh {4nx4nx4n} with n = 3 and u^ = - 0.1 m (for example) is plotted according to the principles of the present invention. Fig. 26 explores the sensitivity of the elements {G, Q IPO} to irregularity of the reference element shape. Since the converged value of the displacement u_E3 is not known it is convenient to define the difference of u_E3 relative to the value Ug₃ predicted for a regular mesh (a = 0). Specifically, the value of Ug₃ for each element for the regular mesh {4nx4nx4n} with n = 3 and U^₃ = - 0.1 m is given by

^UE₃ = ~ 0.0078595 m for (G) , Ug₃ = - 0.0071327 m for (QlPO) . (4.27)

The difference Δ is then defined by a formula of the type (4.7)

for each of the elements using its value of u_E3 so that Δ vanishes for each element when a =

0. The results in Fig. 26 demonstrate that (G) 2610 can be used for three-dimensional irregularity of the reference element shapes and that the predictions of (G) 2610 are similar to those of (QlPO) 2620. Since the solution is not fully converged for this mesh it is expected that the solution will be sensitive to the value of the irregularity parameter a. In particular, the results in Fig. 26 are consistent with the observation that the mesh is more refined under the rigid plate for positive values of a.

Figs. 27a - 27d 2710, 2720, 2730 and 2740, respectively, are graphic representations of the response to indentation of a rigid plate into a nearly incompressible block, plotted according to the principles of the present invention. Fig. 27 shows nonlinear load curves using the regular mesh {4nx4nχ4n} for different values of n. From these figures it can be seen that the results predicted by {G, Q1P0} are similar. In particular, the unphysical stiffness exhibited by (QlPO) in the previous plane strain problem is not detected in this three-dimensional problem.

5. Conclusions

A generalized form (2.11) for the strain energy of inhomogeneous deformations of a Cosserat point element (CPE) has been developed which includes full coupling of bending and torisonal modes. The dependence of the constitutive coefficients By on the metric D_j-

(3.4) have been determined by matching exact small deformation solutions for pure bending (eighteen solutions) and simple torsion (6 solutions). These coefficients are then used with nonlinear strain measures to characterize the response of general irregular element shapes to large deformations.

Ideally, for a fully converged solution the response of a structure to a specified load should be insensitive to irregularities in the element shapes used to mesh the structure. The results here indicate that the main features of this desired response for general shaped elements can be obtained by properly modeling pure bending and simple torsion of parallelepipeds with two right angles. Also, it is recalled (Jabareen and Rubin, 2007c) that insensitivity to element irregularity can best be exhibited by focusing attention on a thin structure that is loaded so that inhomogeneous deformations (like bending and torsion) dominate homogeneous deformations. In particular, plots like Fig. 4c for (I) clearly emphasize undesirable sensitivity to element irregularity. In contrast with standard finite elements, the nodal forces in the generalized CPE are determined by algebraic expressions in terms of derivatives of a strain energy function and no integration is needed over the element region. A number of example problems (also see Jabareen and Rubin, 2007a,c) have been considered which show that the generalized CPE is as accurate as elements based on enhanced strain and incompatible modes and is as robust as elements based on full integration. The plane strain example of indentation of a rigid plate into a nearly incompressible block showed that the flexibility exhibited by the CPE is physical and that the mixed QlPO element predicts unphysical stiffness. Furthermore, the generalized CPE can be used to model 3-D bodies, thin shells and rods and nearly incompressible materials. In addition, the generalized CPE is free of hourglass instabilities that are observed in other element formulations in regions experiencing combined high compression with bending. Consequently, the generalized CPE is truly a robust user friendly element that can be used with confidence to model problems in nonlinear elasticity.

Having described the present invention with regard to certain specific embodiments thereof, it is to be understood that the description is not meant as a limitation, since further modifications will now suggest themselves to those skilled in the art, and it is intended to cover such modifications as fall within the scope of the appended claims. Appendix A: Values of the constitutive coefficients B_ij

The equations of the bending and torsion problems discussed in Section 3 for the element shapes can be solved for the values of B_ij in the strain energy function (2.11) for inhomogeneous deformations and the results were reported in (Jabareen and Rubin, 2008a). Specifically, the nonzero components of the upper diagonal of the symmetric matrices

in (3.24) for a general shaped element are specified by

with {λ₁₂, λ₁₃, λ₂₃} defined by (3.23).

Appendix B: Tire Deformation

Insurance companies are particularly interested in tire safety. In order to improve tire safety, tire deformation methodology is used to design safer tires.

Tire deformation methodology is known in the art.

US Pat. No. 7,472,587, by Loehndorf, et al., Tire Deformation Detection, discloses embodiments related to systems, methods, and apparatuses for detecting tire deformation. In one embodiment, a tire deformation detection system comprises a deformation detection structure, a transmitter, and receiver. US Pat. No 4,475,384, by Christie, Tire Sidewall Deformation Detection Techniques, describes an improved tire manufacturing system, including means for testing the extent of sidewall deformation of a tire. The tire is rotated and measured for lateral run-out of both sidewalls. The resulting data is used to increase the accuracy by which unacceptable sidewall bulges and valleys are detected. US Pat. No 7,412,879, by Serra, et al., Method for Monitoring Tyre Deformations and

Monitoring System Therefor, teaches a method for monitoring deformations in a tire of a vehicle wheel, the wheel including the tire and a rim, the method including: providing the wheel with at least two magnetic-field emitters mutually disposed so as to yield a composite magnetic field and monitoring variations in the magnetic field and correlating the variation with tire deformations.

US Pat. No. 7,497,112 by Kobayakawa, A Method for Estimating a Tire Running Condition and an Apparatus for Effecting the Method and a Tire with Sensors Disposed Therein, describes a method for estimating accurately and stably the change of forces exerted on the tire or change of contact condition with the ground of the tire and an apparatus for putting it in operation.

The various geometric models described by the drawings in the detailed description of the present invention can be correlated with the following definitions in tire deformation, in order to apply the present invention to tire deformation modeling and design. These definitions follow: Slip Angle is the angle between a rolling wheel's actual direction of travel and the direction towards which it is pointing (i.e., the angle of the vector sum of wheel translational velocity v and sideslip velocity u). This slip angle results in a force perpendicular to the wheel's direction of travel ~ the cornering force.

Tire load sensitivity describes the behavior of tires under load. The maximum horizontal force developed should be proportional to the vertical load on the tire. Cornering force is the sideways force produced by a vehicle tire during cornering.

Load Transfer is the imaginary "shifting" of weight around a motor vehicle during acceleration and deceleration.

Weight Transfer is the redistribution of weight supported by each tire during acceleration or deceleration. Weight transfer is a crucial concept in understanding vehicle dynamics. Slip ratio When a vehicle is being driven along a road in a straight line its wheels rotate at virtually identical speeds. The vehicle's body also travels along the road at this same speed.

Scrub Radius is the lateral distance measured in front or rear view between the center of the tire contact patch and the intersection of the steering axis with the ground. Scrub radius and kingpin inclination determine the moment arm about the steering axis for longitudinal (braking and acceleration) forces acting at the tire contact patch.

Tire uniformity parameters

Axes of measurement Tire forces are divided into three axes: radial, lateral, and tangential.

The radial axis runs from the tire center toward the tread, and is the vertical axis running from the roadway through the tire center toward the vehicle. This axis supports the vehicle's weight. The lateral axis runs sideways across the tread.

Radial force variation acts upward to support the vehicle, radial force variation describes the change in this force as the tire rotates under load. As the tire rotates and spring elements with different spring constants enter and exit the contact area, the force will change. RFV, as well as all other force variation measurements, can be shown as a complex waveform. This waveform can be expressed according to its harmonics by applying Fourier Transform (FT).

Lateral force variation acts side-to-side along the tire axle, lateral force variation describes the change in this force as the tire rotates under load.

Tangential force variation acts in the direction of travel, and describes the change in this force as the tire rotates under load. Conicity is based on lateral force behavior. It is the characteristic that describes the tire's tendency to roll like a cone. PIysteer is a parameter based on lateral force behavior. It is the characteristic that is usually described as the tire's tendency to "crab walk", or move sideways while maintaining a straight-line orientation.

Radial runout describes the deviation of the tire's roundness from a perfect circle. RRO can be expressed as the peak-to-peak value as well as harmonic values.

Lateral runout describes the deviation of the tire's sidewall from a perfect plane. LRO can be expressed as the peak-to-peak value as well as harmonic values.

Sidewall Shear Displacement is the displacement of sidewalls from their at rest center line.

Tread Shear Displacement is the amount the tread particles are displaced from their at rest center. By examining the picture under "slip," one can notice that the tread shear displacement is represented by the contact pitch, which is offset.

Total Shear Displacement is the combined amount the tread, sidewall/carcass and belt are displaced from their at rest center line.

Tread Shear Angle the angle of tread particle displacement within the contact patch. Sidewall Shear Angle is the angle of sidewall, carcass and belt (radial tire) displacement from their at rest center line.

Total Shear Angle the combined tread and sidewall/carcass angle of shear.

Trailing Edge Slip Percentage is the percentage of tread particles exceeding the coefficient of friction at the trailing edge of the contact patch in relation to the contact patch's length. Tread Shear Angle the angle in which, during tire deformation, the tread particles are moved in the different of the lateral forces.

Tread Particles do not act independently of each other but rather push and pull one against the next much like the fibers of a cleaning brash. In the case of tire tread depth, a shorter fiber of equal strength will product a higher force for an equal amount of shear.

Claims

We claim:

1. A method for Cosserat Point Element (CPE) modeling of nonlinear elastic materials, utilizing a strain energy function characterizing the response of an eight node 3-D brick element to homogeneous and inhomogeneous deformations, which includes full coupling of bending and torsional modes of deformation, said method comprising: determining the constitutive coefficients which satisfy bending and torsional modes of deformation for three specific reference element geometries; combining said constitutive coefficients to obtain constitutive equations for general reference element shapes which ensure that the strain energy remains a positive definite function for all inhomogeneous deformations of the CPE, said constitutive coefficients characterizing a generalized CPE model to predict the response of a 3-D brick element for finite elastic deformations which is free of hourglass instabilities and locking for nearly incompressible materials, thereby providing a robust modeling tool.

2. The method of claim 1, wherein said modeling is applied to general three dimensional structures.

3. The method of claim 1, wherein said modeling is applied to thin structures.

4. The method of claim 3, wherein thin structures are shells.

5. The method of claim 3, wherein thin structures are beams.

6. The method of claim 1, wherein said modeling is applied to the problem of large elastic deformations of tires.

7. The method of claim 1, further comprising considering a number of exemplary problems which show that the generalized CPE produces results as accurate as enhanced strain and incompatible elements for thin structures.

8. The method of claim 8, wherein said considering produces results free of hourglass instabilities typically predicted by these enhanced elements in regions experiencing combined high compression with bending.

9. The method of claim 1, further comprising reexamining exact linear models of bending and simple torsion of these special elements.

10. The method of claim 9, further comprising including more general coupling of bending and torsional modes, thereby enabling matching models for the two additional bending modes not previously handled.

11. The method of claim 10, further comprising matching these additional two bending modes, thereby significantly improving the response of the generalized CPE for problems of bending of a rhombic plate.

12. The method of claim 1, performed such that no integration is needed over the element region.

13. The method of claim 1, wherein said reference geometry of the CPE is used to generate a mesh of the body being modeled.

14. The method of claim 1, wherein said reference geometry is an eight node 3-D brick element.

15. The method of claim 1, wherein said reference geometry of the body is a cylinder.

16. The method of claim 1, wherein said reference geometry of the body is a circular shell.

17. The method of claim 1, wherein said reference geometry of the body is a rhombic plate.

18. The method of claim 1, wherein said reference geometry of the body is a thin cantilever beam.

19. The method of claim 18, wherein said thin cantilever beam is a thin twisted cantilever beam.

20. A software program (code generator) to generate the code for implementing the method for Cosserat Point Element (CPE) modeling of claim 1, wherein: a general reference geometry of the structure of a 3-D brick element is selected; and the constitutive coefficients of the strain energy function associated with said 3-D brick element are determined, such that the resulting constitutive coefficients ensure that the strain energy function for inhomogeneous deformations remains a positive definite function for said 3-D brick element, said constitutive coefficients characterizing a generalized CPE free of hourglass instabilities, thereby providing a robust modeling tool.

21. A computer program product (code generator) comprising a computer usable medium having computer readable code embodied therein for execution on a general purpose computer, said computer usable medium storing instructions that, when executed by the computer, cause the computer to perform a method for Cosserat Point Element (CPE) modeling of nonlinear elastic materials, utilizing a strain energy function characterizing the response of an eight node 3-D brick element to homogeneous and inhomogeneous deformations, which includes full coupling of bending and torsional modes of deformation, said method comprising: determining the constitutive coefficients which satisfy bending and torsional modes of deformation for three specific reference element geometries; combining these solutions to obtain constitutive equations for general reference element shapes which ensure that the strain energy remains a positive definite function for all inhomogeneous deformations of the CPE, said constitutive coefficients characterizing a generalized CPE model the response of a

3-D brick element for finite elastic deformations which is free of hourglass instabilities and locking for nearly incompressible materials, thereby providing a robust modeling tool.