U.S. patent application number 11/401832 was filed with the patent office on 2007-01-11 for solver for a restrained deformable system with released degrees of freedom.
Invention is credited to Charles-Andre De Hillerin, Sebastien Voegele.
Application Number | 20070007382 11/401832 |
Document ID | / |
Family ID | 37101477 |
Filed Date | 2007-01-11 |
United States Patent
Application |
20070007382 |
Kind Code |
A1 |
De Hillerin; Charles-Andre ;
et al. |
January 11, 2007 |
Solver for a restrained deformable system with released degrees of
freedom
Abstract
A process of using finite element analysis for determining the
behaviour of a system, submitted to predetermined load and
displacement boundary conditions with a degree of freedom relative
to a local coordinate system released at one node of the system.
For a geometrically non-linear solution scheme, a plurality of
incremental steps are carried, for which computing for a finite
element of the system a displacement vector at the nodes of the
finite element is performed by minimizing the energy of the system.
For a geometrically linear solution scheme, a single solving step
is carried, for which similar steps of computing for a finite
element and timely minimizing of the system's energy are performed.
The step of minimizing is carried out, for a node having a released
degree of freedom, in a local frame of reference, and the released
degree of freedom at the node is expressed along one axis of such
frame. Such hybrid reference frames allow expressing and taking
into account in the global formulation the constraints of the
non-released degrees of freedoms in a simple manner, without having
to factorise constraints.
Inventors: |
De Hillerin; Charles-Andre;
(Saint Cloud, FR) ; Voegele; Sebastien; (Blagnac,
FR) |
Correspondence
Address: |
CLIFFORD CHANCE US LLP
31 WEST 52ND STREET
NEW YORK
NY
10019-6131
US
|
Family ID: |
37101477 |
Appl. No.: |
11/401832 |
Filed: |
April 10, 2006 |
Current U.S.
Class: |
244/3.1 ;
244/175 |
Current CPC
Class: |
Y10S 706/919 20130101;
Y10S 706/92 20130101; G06F 30/23 20200101; G06F 2111/04 20200101;
Y10S 715/961 20130101; Y10S 715/964 20130101; G06F 2113/16
20200101; Y10S 706/902 20130101 |
Class at
Publication: |
244/003.1 ;
244/175 |
International
Class: |
F41G 7/00 20060101
F41G007/00 |
Foreign Application Data
Date |
Code |
Application Number |
Apr 8, 2005 |
EP |
05290790.4 |
Mar 21, 2006 |
EP |
06290450.3 |
Claims
1. A process for determining the behaviour of a system submitted to
predetermined load and displacement boundary conditions, using
finite element analysis, the predetermined displacement boundary
conditions including at least one degree of freedom release at one
node of the system; wherein the process comprises a plurality of
incremental steps; each incremental step comprises computing for a
finite element of the system a displacement vector at the nodes of
the finite element, by minimizing an energy of the system when
submitted to a fraction of the predetermined displacement; wherein
the step of minimizing is carried out, for a node having a released
degree of freedom, in a local frame of reference; and a released
degree of freedom at said node is expressed along one axis of the
local frame of reference.
2. The process of claim 1, wherein the system is modelled with
continuous elements for which the displacement vector has only
translational components.
3. The process of claim 1, wherein the system is modelled with
structural elements for which the displacement vector has
translational and rotational components.
4. The process of one of claims 1, wherein a node of a finite
element is provided with a nodal frame of reference; and the step
of computing further comprises computing, for a displaced node of
the finite element, an updated nodal frame of reference by rotating
the nodal frame of reference by rotational components of the
displacement vector at the node.
5. The process of claim 4, further comprising the step of
interpolating a smoothed shape for the system, based on: the
positions of the nodes for the finite elements; and, the nodal
frames of reference.
6. The process of claim 5, further comprising a step of displaying
the interpolated smooth shape.
7. The process of claim 1, wherein said minimized energy of the
system is a function of the difference between a nodal frame of
reference at a node and an elementary frame of reference of a
finite element comprising the node.
8. The process of claim 7, further comprising, during a step of
computing, computing at a node of finite element a rotation
transforming an elementary frame of reference into the nodal frame
of reference; wherein minimized energy depends on the computed
rotation.
9. The process of claim 7, wherein said minimized energy of the
system is a function of the product of a scaling factor by the
difference between a nodal frame of reference at a node and an
elementary frame of reference of a finite element comprising the
node, the process further comprising changing the scaling factor
from one incremental step to another incremental step.
10. The process of claim 9, wherein the scaling factor is changed
manually.
11. The process of claim 9, wherein the scaling factor is changed
according to a number of iterations for minimizing the energy in an
incremental step.
12. The process of claim 9, wherein the scaling factor is changed
according to the difference between an elementary frame of
reference for a finite element and the nodal frames of reference
for the finite element.
13. The process of claim 1, wherein the step of minimizing the
energy is carried out iteratively and comprises: computing a scalar
product of a residual and of a displacement variation; and,
stopping the iterations according to a value of the scalar
product.
14. A computer readable storage medium whose content enables a
computer system to determine the behaviour of a system submitted to
predetermined load and displacement boundary conditions, using
finite element analysis, the predetermined displacement boundary
conditions including at least one degree of freedom release at one
node of the system, by, employing a plurality of incremental steps,
wherein each incremental step comprises computing for a finite
element of the system a displacement vector at the nodes of the
finite element, by minimizing an energy of the system when
submitted to a fraction of the predetermined displacement; wherein
the step of minimizing is carried out, for a node having a released
degree of freedom, in a local frame of reference; and a released
degree of freedom at said node is expressed along one axis of the
local frame of reference.
15. A process for determining the behaviour of a system submitted
to predetermined load and displacement boundary conditions, using
finite element analysis, the predetermined displacement boundary
conditions including at least one degree of freedom release at one
node of the system; wherein, the process comprises a step of
computing for a finite element of the system a displacement vector
at the nodes of the finite element, by minimizing an energy of the
system when submitted to the predetermined displacement; wherein
the step of minimizing is carried out, for a node having a released
degree of freedom, in a local frame of reference, and a released
degree of freedom at said node is expressed along one axis of the
local frame of reference.
16. The process of claim 15, wherein the system is modelled with
continuous elements for which the displacement vector has only
translational components.
17. The process of claim 15, wherein the system is modelled with
structural elements for which the displacement vector has
translational and rotational components.
18. The process of one of claims 15, wherein a node of a finite
element is provided with a nodal frame of reference ; and the step
of computing further comprises computing, for a displaced node of
the finite element, an updated nodal frame of reference by rotating
the nodal frame of reference by rotational components of the
displacement vector at the node.
19. The process of claim 18, further comprising the step of
interpolating a smoothed shape for the system, based on: the
positions of the nodes for the finite elements; and, the nodal
frames of reference.
20. The process of claim 19, further comprising a step of
displaying the interpolated smooth shape.
21. The process of claim 15, wherein said minimized energy of the
system is a function of the difference between a nodal frame of
reference at a node and an elementary frame of reference of a
finite element comprising the node.
22. The process of claim 21, further comprising, during a step of
computing, computing at a node of finite element a rotation
transforming an elementary frame of reference into the nodal frame
of reference; wherein minimized energy depends on the computed
rotation.
23. The process of claim 20, wherein said minimized energy of the
system is a function of the product of a scaling factor by the
difference between a nodal frame of reference at a node and an
elementary frame of reference of a finite element comprising the
node, the process further comprising changing the scaling factor
from one incremental step to another incremental step.
24. The process of claim 23, wherein the scaling factor is changed
manually.
25. The process of claim 23, wherein the scaling factor is changed
according to a number of iterations for minimizing the energy in an
incremental step.
26. The process of claim 23, wherein the scaling factor is changed
according to the difference between an elementary frame of
reference for a finite element and the nodal frames of reference
for the finite element.
27. The process of claim 15, wherein the step of minimizing the
energy is carried out iteratively and comprises: computing a scalar
product of a residual and of a displacement variation; and,
stopping the iterations according to a value of the scalar
product.
28. A computer readable storage medium whose content enables a
computer system to determine the behaviour of a system submitted to
predetermined load and displacement boundary conditions, using
finite element analysis, the predetermined displacement boundary
conditions including at least one degree of freedom release at one
node of the system, by, employing a step of computing for a finite
element, by minimizing an energy of the system when submitted to
the predetermined displacement; wherein the step of minimizing is
carried out, for a node having a released degree of freedom, in a
local frame of reference; and a released degree of freedom at said
node is expressed along one axis of the local frame of reference.
Description
CROSS REFERENCE TO RELATED APPLICATIONS
[0001] This application claims priority to European Patent
Application No. 06290450.3, filed on Mar. 21, 2006, which claims
priority to European Patent Application No. 05290790.4, filed on
Apr. 8, 2005.
FIELD OF THE INVENTION
[0002] The invention relates to the field of computer programs and
systems, more specifically to part design programs and such
constraints-based systems.
BACKGROUND OF THE INVENTION
[0003] A number of systems and programs are offered on the market
for the design of parts or assemblies of parts, such as the one
provided by the applicant under the trademark CATIA. These
so-called computer-aided design (CAD) systems allow a user to
construct and manipulate complex three-dimensional (3D) models of
parts or assemblies of parts. These systems and programs use
various constraints for defining models. The set of constraints is
solved by the system when the model is edited. The program or
system used for solving the set of constraints is generally called
a "solver." Such solvers are used in CAD/CAM/CAE systems or more
generally in any system using constraints for defining objects of
any kind. A solver is adapted to the objects being designed and to
the type of constraints applied to these objects. A solver such as
the one used in CATIA is adapted for the design of solid objects,
with constraints comprising dimensional constraints, stresses,
contacts between objects and the like.
[0004] There is a need for a solver making it possible to simulate
the deformed shapes of deformable solid physical systems subjected
to predetermined load and displacement boundary conditions.
[0005] Some flexible systems and objects usually designed in
CAD/CAM/CAE systems are characterized by their high deformability
("geometrical non-linearity") in the sense that distinct points can
independently undergo arbitrarily large rotations under the action
of loads. In other words, for such flexible systems, the deformed
shape of the system may substantially differ from the non-deformed
shape of the system. On the other hand, geometrically linear
solvers in existing CAD/CAM/CAE systems work under the assumption
that the deformed shape of the object is close to the non-deformed
shape.
[0006] More generally, systems undergoing large rotations, such
that the rotation angles in radians cannot be approximated by the
tangent--typically 5 degrees, but this figure may vary depending on
the amount of modelling error tolerated--are usually considered as
"geometrically non-linear systems." Alternatively, one could say
that a system is non-linear when a "large strain" definition such
as the Green-Lagrange measure: ij GL = 1 2 .function. [
.differential. u i .differential. X j + .differential. u j
.differential. X i + .differential. u k .differential. X i .times.
.differential. u k .differential. X j ] ##EQU1## cannot be replaced
by the "small strain" definition of the strain tensor components:
ij LIN = 1 2 .times. ( .differential. u i .differential. X j +
.differential. u j .differential. X i ) ##EQU2## that is, when the
product term in the Green-Lagrange measure cannot be neglected.
[0007] In general, there is a need in solvers to allow the partial
release of degrees of freedom at some restrained points of the
system. For instance, one may provide a sliding point, at which the
system displacement is restrained in two directions while being
allowed to slide along one direction, and possibly to rotate around
the sliding direction; in this example, the system has one degree
of freedom in translation and one degree of freedom in rotation at
that point. As discussed in detail below, the prior art (e.g.: Use
of the FEM for the Design of Flexible Parts, C. A. de Hillerin,
Proceedings of NAFEMS Word Congress 1999, pp. 345-356) teaches
factorising the constraints to account for such situations.
[0008] This need, which is present for slender body flexible
systems, which are non-linear systems, also arises for linear
systems, in which the deformed shape of the system is assumed to be
close to the non-deformed shape.
[0009] Another feature of slender body flexible systems is that
elements of the system in a finite element analysis are subject to
displacements having translational and rotational components; the
elements are so-called "structural" elements; other examples of
structural elements include beams and shells. The need to allow the
partial release of degrees of freedom at constrained points is also
present in systems modelled using so-called "continuous" elements,
where the displacement vector for an element of the system only has
translational components. Examples of continuous systems include
bars, membranes and solids; the distinction between structural and
continuous system is discussed in Finite Elements Procedures in
Engineering Analysis, K-J. Bathe, Prentice-Hall 1982. Complex
systems may, of course, contain both types of elements in the same
model.
[0010] The prior art solution of factorising the constraints may
involve long and complex computation, especially when the released
degrees of freedom do not correspond to the axes of the global
reference frame. Thus, there is a need for a solver allowing
partial release of degrees of freedom at some restrained points of
the system, while making the computation simpler.
SUMMARY OF THE INVENTION
[0011] According to one embodiment of the present invention, there
is provided a solution to this problem. This solution can be
applied to both structural and continuum element formulations and
to both non-linear and linear solution schemes.
[0012] For continuum systems, a first embodiment of the invention
relates to the way in which, within the framework of a static
finite element solution, one can treat a system of linear
constraint equations, where the number of constraint relations
involving the global translational degrees of freedom of a partly
restrained node is 1 or 2. These constraint relations represent the
partial restraining of the specified nodal translational
displacement along one or two oblique directions non-aligned with
the global coordinate axes.
[0013] For such situations, this first embodiment suggests
replacing the standard <<constraints factorization>>
procedure by a new procedure.
[0014] The standard "constraints factorization" procedure consists
in generating a set of coupled constraint equations, written in
terms of global displacement components, then transforming it into
an equivalent form, from which all dependent degrees of freedom are
eliminated, then substituting this form in the global system of
equilibrium equations (assumed to be constructed from elementary
stiffness matrices and load vectors rotated in the global reference
frame prior to assembly), and finally projecting this system on the
nullspace of the set of coupled constraint equations, thereby
eliminating all dependent degrees of freedom and unknown Lagrange
multipliers and transforming the global system of equilibrium
equations into a more compact and invertible system of
equations.
[0015] In the new procedure according to the first embodiment, the
global components of the translational displacement of each partly
restrained node are first converted to local components relative to
an orthonormal reference frame constructed from the local oblique
restraint directions. In local components, the constraint relations
involving the translational displacement of such a partly
restrained node are then expressed in a simpler, uncoupled form.
The construction of the corresponding system of equilibrium
equations requires the computation of "hybrid" elementary stiffness
matrices and load vectors prior to assembly.
[0016] The "hybrid" elementary stiffness matrices and load vectors
are obtained by "hybrid" matrix and vector rotation, i.e. by
selectively rotating diagonal blocks of the stiffness matrix
corresponding to different nodes in different reference frames and
by appropriately propagating the effect on non-diagonal blocks,
with an analogous treatment for the load vector.
[0017] Once the appropriate system of equilibrium equations has
been assembled, the uncoupled constraint equations are taken into
account by simple row and column elimination operations, upon
direct inspection of the assembled stiffness matrix and load
vector. The resulting system of equations is termed "hybrid" since
its unknowns are not displacement components relative to a unique
reference frame. Nevertheless, all dependent degrees of freedom and
unknown Lagrange multipliers have been eliminated and the
transformed system of equilibrium equations is again compact and
invertible.
[0018] The process applies not only to non-linear solutions but
also to linear solutions, for which the absence of incremental
steps implies that restraints have to be accounted for only
once.
[0019] For structural elements (including beam, plate and shell
elements, in addition to the "cable" elements described in the
examples below), the main difference with continuum elements
consists in the fact that structural elements do exhibit rotational
nodal degrees of freedom in addition to translational nodal degrees
of freedom. For structural elements, a second embodiment of the
invention relates to the way in which, within the framework of a
static finite element solution, one can treat a system of linear
constraint equations, where the number of constraint relations
involving the global translational degrees of freedom of a partly
restrained node is 1 or 2, and/or the number of constraint
relations involving the global rotational degrees of freedom of a
partly restrained node is 1 or 2. These constraint relations
represent the partial restraining of the specified nodal
translational and/or rotational displacement along oblique
directions non-aligned with the global coordinate axes.
[0020] For such situations, the second embodiment of the invention
suggests replacing the standard <<constraints
factorization>> procedure by a new procedure, similar to the
one described above. Aside from this specific difference (which is
related to the fact that structural systems with linearized
rotation kinematics require updating in the case of non-linear
incremental solutions), everything else remains the same regarding
the hybridation process for the treatment of oblique constraints
with partial relaxation.
[0021] Specifically, the invention provides in a first embodiment a
process for determining the behaviour of a system submitted to
predetermined load and displacement boundary conditions, using
finite element analysis, the predetermined displacement boundary
conditions including at least one degree of freedom release at one
node of the system wherein the process comprises a plurality of
incremental steps, wherein each incremental step comprises
computing for a finite element of the system a displacement vector,
at the nodes of the finite element, by minimizing an energy of the
system when submitted to a fraction of the predetermined
displacement; wherein the step of minimizing is carried out, for a
node having a released degree of freedom, in a local frame of
reference, and wherein a released degree of freedom at said node is
expressed along one axis of the local frame of reference.
[0022] In a second embodiment, the invention provides a process for
determining the behaviour of a system submitted to predetermined
load and displacement boundary conditions, using finite element
analysis, the predetermined displacement boundary conditions
including at least one degree of freedom release at one node of the
system, wherein the process comprises a step of computing for a
finite element of the system a displacement vector, at the nodes of
the finite element, by minimizing an energy of the system when
submitted to the predetermined displacement; wherein the step of
minimizing is carried out, for a node having a released degree of
freedom, in a local frame of reference, and wherein a released
degree of freedom at said node is expressed along one axis of the
local frame of reference.
[0023] In the first and second embodiments, one may provide one or
more of the following additional features: [0024] the system is
modelled with continuous elements for which the displacement vector
has only translational components; [0025] the system is modelled
with structural elements for which the displacement vector has
translational and rotational components; [0026] a node of a finite
element is provided with a nodal frame of reference and the step of
computing further comprises computing, for a displaced node of the
finite element, an updated nodal frame of reference by rotating the
nodal frame of reference by rotational components of the
displacement vector at the node; [0027] the process further
comprises the step of interpolating a smoothed shape for the
system, based on the positions of the nodes for the finite elements
and on the nodal frames of reference; [0028] the process further
comprises a step of displaying the interpolated smooth shape;
[0029] the minimized energy of the system is a function of the
difference between a nodal frame of reference at a node and an
elementary frame of reference of a finite element comprising the
node; [0030] the process further comprises, during a step of
computing, computing at a node of finite element a rotation
transforming an elementary frame of reference into the nodal frame
of reference, minimized energy depending on the computed rotation;
[0031] the minimized energy of the system is a function of the
product of a scaling factor by the difference between a nodal frame
of reference at a node and an elementary frame of reference of a
finite element comprising the node, the process further comprising
changing the scaling factor from one incremental step to another
incremental step; the scaling factor may be changed manually or
according to a number of iterations for minimizing the energy in an
incremental step; else, the scaling factor may be changed according
to the difference between an elementary frame of reference for a
finite element and the nodal frames of reference for the finite
element; [0032] the step of minimizing the energy is carried out
iteratively and comprises computing a scalar product of a residual
and of a displacement variation and stopping the iterations
according to a value of the scalar product.
[0033] Last, the invention provides a computer program product
residing on a computer-readable medium, comprising computer program
code means adapted to run on a computer all the steps of such a
process.
[0034] A system embodying the invention will now be described, by
way of non-limiting example, and in reference to the accompanying
drawings.
BRIEF DESCRIPTION OF THE FIGURES
[0035] FIG. 1 is a schematic view of the flexible system, in the
starting configuration and in the final configuration, showing a
number of intermediate configurations;
[0036] FIG. 2 is a partial view representation of the finite
element model for the flexible system, in the starting
configuration and in the first interim configuration;
[0037] FIG. 3 is an enlarged view of a finite element of FIG.
2;
[0038] FIGS. 4 and 5 are view of a finite element with various
reference frames;
[0039] FIG. 6 is a view of two finite elements with the frames used
in case where a node of the system has a degree of freedom;
[0040] FIGS. 7 to 22 show various examples of the process.
DETAILED DESCRIPTION OF THE INVENTION
[0041] The description of slender-body flexible systems in
reference to FIGS. 1-5 is useful for understanding the general
operation of a solver, even though it does not address specifically
the problem of limited degrees of freedom, but rather co-rotational
updating in the prior art incremental total Lagrangian formulation.
For the sake of clarity, the process disclosed in the NAFEMS World
Congress paper (thereafter the NAFEMS paper) is discussed again
here. In the first part of the present description, one considers
the case where the predetermined reference displacement of the
flexible system at both ends is fully determined.
[0042] The purpose of the process is to determine the behaviour of
the flexible system, when submitted to predetermined displacement
boundary conditions. Reference 2 in FIG. 1 designates the starting
configuration of the flexible system. At both ends of the system,
frames of reference 4 and 6 are represented. These frames of
reference are local frames of reference and are representative of
the orientation of the flexible system at both ends. One
understands that the word "displacement" thus includes not only the
position of the ends of the flexible system, but also the
orientation of the flexible system. Displacement of a point thus
includes six degrees of freedom.
[0043] Reference 12 in FIG. 1 designates the final configuration of
the flexible system, while references 14 and 16 designate the
reference frames at both ends of the flexible system. At the first
end--reference frames 6 and 16--the displacement includes a
displacement in position, the axis of reference frames 6 and 16
being parallel. At the second end--reference frames 4 and 14--the
displacement includes a displacement in position and orientation,
the axis of reference frame 14 being turned compared to the axis of
reference frame 4.
[0044] FIG. 1 further shows interim configurations 20, 22 and 24 of
the flexible system. The number of interim configuration in FIG. 1
is solely selected for the purposes of explanation. For passing
from one configuration to the next one--which represents one
incremental step--a fraction of the predetermined displacement is
applied to the flexible system. In the example of FIG. 1, with
three interim configurations and therefore four incremental steps,
one passes from one configuration to the next one by applying to
the flexible system one fourth of the predetermined displacement.
For instance, in the first incremental step, one passes from the
initial configuration 2 to the first interim configuration 20, by
applying to the ends of the flexible system one fourth of the
predetermined displacement--in translation and in rotation.
[0045] Each interim configuration and the final configuration of
the flexible system is computed using geometrically non-linear
finite element analysis; specifically, the flexible system is
represented as a series of finite elements. For the sake of
explanation, the first incremental step is considered: in this
first incremental step, one passes from the starting configuration
2 to the first interim configuration 20. FIG. 2 shows a partial
view of the starting configuration 2 with the reference frame 4 and
some of the elements. It shows a number of elements 30, 35, 40, 45;
each element extends between two nodes 31, 36, 41, 46 and 51. Node
31 is actually the end of the flexible system. Each element 30, 35,
40, 45 is associated with an elementary frame of reference,
referenced 32, 37, 42 and 47 on FIG. 2. FIG. 2 further shows a
partial view of the first interim configuration 20, each of the
nodes; elements and frames of reference having a reference number
increased by 30.
[0046] In each incremental step, one minimizes the total potential
energy of the flexible system, when submitted to the fraction of
the predetermined displacement. The total potential energy of the
flexible system is computed as the sum of the total potential
energy of each of the finite elements forming the flexible system.
This is carried out thanks to an iterative process discloses in the
NAFEMS article. This iterative process may also be carried out as
discussed below.
[0047] In each incremental step, one uses a total Lagrangian
formalism with co-rotational update. As discussed in section
4.3.2.1 of the NAFEMS article, displacements within a given
incremental step--or "palier" in the wording used in the NAFEMS
article--are measured relative to the configuration corresponding
to the solution obtained in the previous incremental step. In other
words, in one incremental step, minimisation of the total potential
energy of the flexible system is carried out, for each element, in
the elementary reference frame associated to the element in the
configuration of the flexible system at the beginning of the
incremental step. In the example of FIG. 2, the configuration at
the beginning of the incremental step is the starting configuration
2. One will therefore iteratively minimize the total potential
energy of the flexible system for obtaining the first interim
configuration, in the frames of reference corresponding to the
starting configuration--that is in elementary frames of reference
32, 37, 42 and 47 of FIG. 2.
[0048] However, for the next incremental step, one will consider
the updated frames of reference, corresponding to the first interim
configuration. In other words, the following incremental step
passes from first interim configuration 20 to second interim
configuration 22. Computation of the second incremental step is
carried out in the elementary frames of references 62, 67, 72 and
77 associated with the first interim configuration. This makes it
necessary to compute new frames of references at the end of each
incremental step, before the next incremental step. The NAFEMS
article does not provide any detail on the computation of the
updated reference frame for each element. In other words, even
through co-rotational update is suggested, no explanation is
provided as regards the carrying out of this co-rotational
update.
[0049] The invention suggests proceeding with the co-rotational
update for one given finite element, based on the components of the
displacement vectors computed at the nodes of the finite element,
without considering the bending rotational components of the
displacement vectors.
[0050] The solution of the invention has the advantage of ensuring
that in the updated elementary frame of reference, one of the axes
is parallel to the line joining the nodes of the finite element.
This property is maintained in each incremental step and makes all
computations similar, hence simpler and easier.
[0051] FIG. 3 is an enlarged view of finite element 40 and 70 of
FIG. 2. It shows the nodes 41 and 46 in the starting configuration,
with finite element 40 extending between nodes 41 and 46. The
elementary frame of reference 42 of finite element 40 has axes x, y
and z. The x-axis is parallel to the line joining nodes 41 and
46.
[0052] Minimizing the total potential energy of the flexible system
provides a displacement vector for each node of the flexible
system. The displacement vector comprises translational components
(u.sub.x, u.sub.y, u.sub.z) and rotational components
(.theta..sub.x, .theta..sub.y, .theta..sub.z). The translational
components (u.sub.x, u.sub.y, u.sub.z) are representative of the
displacement in position of each node; the translational component
is easily understood and in the example of FIG. 3 corresponds to
the displacement in position of node 41 to node 71, or of node 46
to 76.
[0053] The rotational components (.theta..sub.x, .theta..sub.y,
.theta..sub.z) are representative of the displacement in
orientation of each node; this rotational component is
representative of the fact that the system is a flexible system, so
that the orientation of the system at a given node would change in
one incremental step.
[0054] The rotational components include bending rotational
components, which are representative of the bending applied to the
flexible system at a given node in the configuration obtained at
the end of the incremental step.
[0055] The rotational components further include a torsional
rotational component, which is representative of the torsion
applied to the flexible system at the given node.
[0056] If the rotational components (.theta..sub.x, .theta..sub.y,
.theta..sub.z) at node 41 are expressed in the frame of reference
42 of element 40--in which the x-axis is parallel to the line
extending from node 41 to node 46, the torsional component is the
component .theta..sub.x measured along axis x. The bending
rotational components are components .theta..sub.y and
.theta..sub.z. In other words, bending rotational components are
representative of the rotation of a node along axes transverse to
the x-axis of the elementary frame of reference, which is parallel
to the line joining the nodes of the finite element.
[0057] The torsional rotational component is .theta..sub.x and is
representative of the torsion applied at the node the along the
x-axis of the elementary frame of reference.
[0058] In a first example, the updated elementary frame of
reference is computed solely based on the translational components
of the displacement vectors at the nodes of the finite element.
This example is adapted to flexible systems which deform in two
dimensions only; in such cases, the flexible system is not
submitted to any torsional rotational constraint, but only bends in
the 2D plane. An example of such a plane deformation is provided by
a horizontal beam clamped at both ends, acted upon by gravity and
subjected to an imposed horizontal compressive displacement at one
end
[0059] In this case, the updated elementary frame of reference may
be computed as follows. One considers, in the configuration at the
beginning of the incremental step, a finite element P.sub.1P.sub.2
joining two nodes P.sub.1 and P.sub.2. The elementary frame for
this finite element has three axes x, y and z. The x-axis is
parallel to the line (P.sub.1P.sub.2) joining the two nodes. The
y-axis is perpendicular to the x-axis in the 2D plane. The z-axis
is perpendicular to the 2D-plane.
[0060] Let (.theta..sub.x1, .theta..sub.y1, .theta..sub.z1) be the
translational components of the displacement vector computed at the
first node P.sub.1 of the finite element and measured in the
elementary reference frame of the finite element. Let
(.theta..sub.x1, .theta..sub.y1, .theta..sub.z1) be the rotational
components of the displacement vector computed at the first node
P.sub.1 of the finite element. Similarly, let (u.sub.x2, u.sub.y2,
u.sub.z2) and (.theta..sub.x2, .theta..sub.y2.theta..sub.z2) be the
rotational and translational components of the displacement vector
for the second node P.sub.2, measured in the elementary reference
frame of the finite element. P.sub.1' and P.sub.2' are the
displaced nodes, computed in the incremental step.
[0061] Since the flexible system in the 2D case does not have any
torsional constraint, .theta..sub.x1=.theta..sub.x2=0. The
torsional rotational constraints at the first and second nodes are
zero. The flexible system moves within the 2D plane, so that
u.sub.z1=u.sub.z2=0.
[0062] The bending rotational component of the displacement vectors
at the first and second nodes are thus expressed by .theta..sub.z1
and .theta..sub.z2.
[0063] In this example, the updated elementary frame of reference
is computed as follows. The x-axis of the elementary frame of
reference is set parallel to the line (P.sub.1'P.sub.2') joining
displaced nodes P.sub.1' and P.sub.2'. The x-axis is oriented in
the direction from P.sub.1' to P.sub.2'. The y-axis is
perpendicular to the x-axis, in the 2D plane where the flexible
system is moving. The z-axis is again perpendicular to the 2D
plane. For computing the updated elementary frame of reference, the
rotational components of the displacement vectors are not taken
into account.
[0064] In a second example, the updated elementary frame of
reference is computed based on the translational components and on
the torsional rotational component of the displacement vectors at
the nodes of the finite element. This example is adapted to
flexible systems which deform in three dimensions. This would for
instance apply to cables, for which a torsional constraint applied
to the cable results in a 3D deformation of the cable, and more
generally to any situation where a straight neutral line deforms
into a non-planar curve.
[0065] The displacement vectors for the nodes P.sub.1 and P.sub.2
are referenced as in the first example; in this example, there is
no reason to assume that .theta..sub.x1=.theta..sub.x2=0 or
u.sub.z1=u.sub.z2=0. The updated elementary frame of reference is
computed as follows. One first considers the torsional component
.theta..sub.x1 and .theta..sub.x2 of nodes P.sub.1 and P.sub.2. The
updated elementary frame of reference is computed by [0066]
applying to the elementary frame of reference a rotation along the
x-axis--that is along the line (P.sub.1 P.sub.2) joining nodes
P.sub.1 and P.sub.2--having an angle equal to the mean torsional
rotational component (.theta..sub.x1+.theta..sub.x2)/2 of the nodes
P.sub.1 and P.sub.2, [0067] applying to the rotated frame a
rotation equal to the minimal rotation that turns a unitary vector
directed from node P.sub.1 to node P.sub.2 into a unitary vector
directed from node P.sub.1' to node P.sub.2'; this rotaiton is
computed based on the translational components (u.sub.x1, u.sub.y1,
u.sub.z1 ) and (u.sub.x2, u.sub.y2, u.sub.z2) of the displacement
vectors, and does not take into account the rotational components
(.theta..sub.x1, .theta..sub.y1, .theta..sub.z1) and
(.theta..sub.x2, .theta..sub.y2, .theta..sub.z2) of the
displacement vectors. In other words, the elementary frame of
reference used during the incremental step is first rotated by an
amount corresponding to the mean torsional rotational component
applied to the node, and the rotated frame is then further rotated
so that the x-axis of the updated elementary frame of reference is
parallel to the line (P.sub.1'P.sub.2') joining the transformed
nodes P.sub.1' and P.sub.2'.
[0068] As in the first example, the x-axis of the updated frame of
reference is parallel to the line joining the displaced nodes
P.sub.1' and P.sub.2'; however, contrary to the first example, the
y- and z-axes are rotated to take into account the torsional
displacement applied to the nodes, so that the xy-plane has no
reason of remaining constant.
[0069] In both examples, the updated elementary frame of reference
is computed without considering the bending rotational components
of the displacement vectors at the node of the finite element.
Also, in both examples, the x-axis of the elementary frame of
reference remains parallel to the line joining the nodes of the
finite element.
[0070] The process disclosed above provides, for each finite
element, an updated frame of reference. This updated frame of
reference may be used for the iterative computation in the next
incremental step. The update of the elementary frames of reference
ensures that the approximation used in the finite element analysis
remains valid. Specifically, the cable kinematics u(x, y,
z)=u.sup.0(x)-y..theta..sub.z(x)+z..theta..sub.y(x) v(x, y,
z)=v.sup.0(x)-z..theta..sub.y(x) w(x, y,
z)=w.sup.0(x)+y..theta..sub.z(x)
[0071] disclosed in section 4.1.1.1 of the NAFEMS article remains
valid while rotations .theta..sub.y and .theta..sub.z are small
enough, so that sin .theta..apprxeq..theta.; this assumption is
valid when the elementary frames of reference are updated as
discussed above.
[0072] The co-rotational update of the elementary frames of
reference does not take into account the bending components of the
displacement vectors at each of the nodes. The bending components
of the displacement vectors are taken into account thanks to the
process of updating nodal reference frames.
[0073] In addition, the incremental step provides a finite element
representation of the flexible system, which is a
slope-discontinuous representation of the system. For obtaining a
smooth--and slope-continuous--representation of the flexible
system, the invention suggests using nodal frames of reference. A
nodal frame of reference is associated with each node of the
flexible system. The nodal frames are updated separately from the
elementary frames of reference, by taking into account the bending
rotational components of the displacement vectors.
[0074] FIG. 3 further shows the nodal frames of reference. In the
starting configuration 2, nodal frames of reference are defined as
follows: [0075] the x-axis of each nodal frame of reference is
tangent to the flexible system; [0076] the y-axis of each nodal
frame of reference is representative of the torsion of the cable
system along the neutral line; in other words, one selects a given
arbitrary direction of the y-axis at one end of the flexible
system, e.g. at the first node 31. For the next node, the y-axis is
determined as being perpendicular to the y-axis, taking into
account the torsional constraint applied to the flexible system.
Thus, if the flexible system is without any torsion in a plane, the
y-axes for the various nodal frames of reference either lie all in
the plane or are all orthogonal to the plane.; [0077] the z-axis of
each nodal frame of reference is such that all axes of a frame are
perpendicular.
[0078] At each incremental step, nodal frames of reference are
updated as follows. The nodal frame of reference at node P.sub.1 is
updated by being rotated by the rotational components
(.theta..sub.x1, .theta..sub.y1, .theta..sub.z1) of the
displacement vector at the node. The origin of the updated nodal
frame of reference is at the displaced node P.sub.1'. This is
represented in FIG. 3, with nodal frames of reference 43 and 48 for
nodes 41 and 46 respectively, as well as nodal frames of reference
73 and 78 for nodes 71 and 76 respectively. FIG. 3 further shows,
in interrupted lines with two dots, the smoothed configuration of
the flexible system. The update of the nodal frames of reference
takes into account the bending and torsional rotational components
of the displacement vector.
[0079] Providing the nodal frames of reference makes it possible,
at each incremental step, to compute the shape of the flexible
system. The shape of the flexible system is computed, based on the
positions of the nodes and on the nodal reference frames.
Specifically, the shape of the flexible system is interpolated
using the position of the nodes and the nodal frames of reference
at the various nodes. The following rules of interpolation may be
used: [0080] the flexible system should pass through each node;
[0081] at each node, the flexible system should be tangent to the
x-axis of the nodal frame of reference; [0082] at each node, the
y-axis and z-axis are representative of the torsional constraints
applied to the cable.
[0083] For interpolating the shape of the flexible system, one may
for instance use state of the art interpolation tools, such as
cubic spline interpolation.
[0084] The step of interpolating is preferably carried out before
displaying to a user the deformed shape of the flexible system.
Thus, the user may be provided with a smoothed view of the deformed
flexible system, and not with a finite element view of the flexible
system.
[0085] In each incremental step, the initial stresses applied to
the flexible system at the beginning of the step may be taken into
consideration, thanks to the nodal frames of reference. As
explained above, the elementary frames of reference are updated
without taking into account the bending rotational components of
the displacement vectors applied to the nodes of finite elements.
On the other hand, the nodal frames of reference are updated based
on all rotational components of the displacement vectors applied to
the nodes. After updating the elementary frames and reference and
the nodal frames of reference, the angular difference at one node
between the elementary frame of reference--on one side of the other
of the node--and the nodal frame of reference is representative of
the amount of the bending rotational components. This difference is
also representative of the amount of initial stresses applied to
the flexible system.
[0086] The invention also suggests using the nodal frames of
reference as a means of measuring the initial stresses applied to
the flexible system at the beginning of an incremental step. This
is carried out by comparing, at each node of each finite element of
the flexible system, the nodal frame of reference to the elementary
frame of reference of the finite element. The comparison--or the
difference between nodal frames of reference and elementary frames
of reference--is representative of the initial stresses applied to
the flexible system. One possible comparison consists in computing
the rotation transforming an elementary frame of reference into a
nodal frame of reference. The angle of this rotation is
representative of the nodal initial stresses and may be used for
computing the corresponding contribution to the energy to be
minimized during the incremental step. One may consider either the
differences at one given node between the nodal frame of reference
and the elementary frame of reference for the two finite elements
comprising the node. Alternatively, one may consider the
differences between the elementary frame of reference for a given
element and the nodal frames of reference of the nodes forming the
finite element. Both examples amount to the same result or to
proportional results, when the difference is summed over the whole
flexible system.
[0087] For instance, for plane motion in the x-y plane, one may use
the following formula for the elementary strain energy: W int = (
EI zz / 2 .times. L 3 ) .function. [ 12 .times. ( ( .DELTA. .times.
.times. v ) 2 - L .times. .times. .DELTA. .times. .times. v
.function. ( ( .theta. z .times. .times. 1 + .theta. z .times.
.times. 1 0 ) + ( .theta. z .times. .times. 2 + .theta. z .times.
.times. 2 0 ) ) + 4 .times. L 2 ( ( .theta. z .times. .times. 1 +
.theta. z .times. .times. 1 0 ) 2 + ( .theta. z .times. .times. 1 +
.theta. z .times. .times. 1 0 ) .times. ( .theta. z .times. .times.
2 + .theta. z .times. .times. 2 0 ) + ( .theta. z .times. .times. 2
+ .theta. z .times. .times. 2 0 ) 2 ] ##EQU3##
[0088] In this formula, the superscripted terms are representative
of the angular difference, at the element nodes, between the nodal
frames of reference and the corresponding elementary frame of
reference at the beginning of the incremental step. The
non-superscripted terms are representative of the current iteration
nodal displacements.
[0089] The use, in the energy to be minimized, of terms
representative of the difference between nodal frames of reference
and elementary frames of reference, may cause instabilities in the
computation of the behaviour of the flexible system; such
instabilities appear as wide changes in the shape of the flexible
system from one incremental step to the next. The use of such terms
may also make it more difficult to find a converged solution during
an incremental step.
[0090] For avoiding or limiting these problems, the corresponding
term or terms may be multiplied by a scaling factor .mu..sub.IS;
the scaling factor is comprised between 0 and 1. The exemplified
formula for the elementary strain energy would then become W int =
( EI zz / 2 .times. L 3 ) [ 12 .times. ( ( .DELTA. .times. .times.
v ) 2 - L .times. .times. .DELTA. .times. .times. v .function. ( (
.theta. z .times. .times. 1 + .mu. IS .times. .theta. z .times.
.times. 1 0 ) + ( .theta. z .times. .times. 2 + .mu. IS .times.
.theta. z .times. .times. 2 0 ) ) + 4 .times. L 2 ( ( .theta. z
.times. .times. 1 + .mu. IS .times. .theta. z .times. .times. 1 0 )
2 + ( .theta. z .times. .times. 1 + .mu. IS .times. .theta. z
.times. .times. 1 0 ) .times. ( .theta. z .times. .times. 2 + .mu.
IS .times. .theta. z .times. .times. 2 0 ) + ( .theta. z .times.
.times. 2 + .mu. IS .times. .theta. z .times. .times. 2 0 ) 2
##EQU4##
[0091] When the scaling factor .mu..sub.IS is 1, initial stresses
are taken into account; when the scaling factor .mu..sub.IS is 0,
the initial stresses are not taken into account. The scaling factor
.mu..sub.IS is preferably changed over the incremental steps,
according to the results. One would typically decrease the scaling
factor when no converged solution may be found during an
incremental step, or when the flexible system proves to be
unstable. This may be carried out by the user; alternatively, the
scaling factor could be adapted automatically, e.g. according to
the current number of iterations in the incremental step. In
another example, the scaling factor depends on the value of the
difference between the mean nodal frame of reference of an element
(average of the nodal frames of reference at its end nodes) and its
elementary frame of reference; this example uses the fact that a
small difference is also indicative of the convergence of the
results of the iterative process in the incremental step; in other
words, the smaller the difference, the better the finite element
simulates the behaviour of the flexible system.
[0092] FIGS. 4 and 5 show examples of a finite element at the end
of an incremental step. Both figures are in 2D for the sake of
explanation and show the finite element 70, nodal frames 73 and 78,
elementary frame 72, as well as a mean nodal frame 80. As explained
above, the mean nodal frame is computed as the average of the nodal
frames; the origin is the middle of the nodes--thus being the same
origin as the elementary frame. Each vector of the mean nodal frame
is the normalized sum of the corresponding vectors of the nodal
frames. In FIG. 4, the angular difference between the elementary
frame 72 and the mean nodal frame 80 is small--some degrees in the
example. This is representative of a good convergence and of a
correct selection of the .mu..sub.IS scaling factor. In FIG. 5, the
angular difference between the elementary frame 72 and the mean
nodal frame 80 is important--around 45.degree.--which is
representative of a too high value of the .mu..sub.IS scaling
factor. In a three-dimensional example, the difference between the
mean nodal frame and the elementary frame of reference may be
computed as the angle of the rotation transforming the elementary
frame of reference into the mean nodal frame. A value of this angle
less than 10.degree. is representative of a correct value of the
scaling factor and of appropriate convergence.
[0093] As discussed above, in each incremental step, one proceeds
with an iterative search for equilibrium solutions. Each iteration
is a two-phase process; the first phase of the process consists in
finding a descent direction over the energy surface; the second
phase of the process consists in line-searching the descent
direction for finding an optimal step. The finding of an optimal
step requires a stopping criterion; section 4.3.2.4 of the NAFEMS
paper suggests evaluating at each iteration the residual vector and
the total energy. Separate convergence checks are then performed on
the norm of the residual and on the monotony of the total energy
variation. This stopping criterion, based on both the norm of the
residual and on the monotony of the energy variation, is effective;
however, it may be useful to provide less stringent stopping
criterion; this is especially useful where the flexible system
passes through an unstable state, before finding an
equilibrium.
[0094] The invention suggests using, as a stopping criterion, a
residual energy norm defined as the scalar product between the
residual and the displacement variation. Using such a stopping
criterion allows the system to pass more easily through unstable
states. This is particularly visible when the system undergoes a
bifurcation and passes brutally from one given configuration to a
very different configuration. Using the proposed stopping criterion
makes it easier to find a solution during the iterative minimizing
of energy carried out in each incremental step; accordingly,
computations are shorter. The use of this stopping criterion is
independent of the corotational update disclosed above.
[0095] In the description of the process up to now, the
predetermined reference displacement of the flexible system at both
ends therefore is fully determined--in other words, position and
orientation of the flexible system are fully constrained. In the
invention, one considers the case where the flexible system has
some degrees of freedom, even at its ends. One example, for
flexible systems, is where one point of the flexible system--often
one end thereof--is fixed in position, but free in rotation; this
means that the point is fixed in translation, but has three degrees
of freedom in rotation. Another example is a sliding point for the
flexible system: the system is free to slide along one direction,
and possibly to rotate around the sliding direction; in this
example, the flexible system has one degree of freedom in
translation and one degree of freedom in rotation.
[0096] This problem exists for slender body flexible systems, which
are non-linear and structural systems. This problem also exists for
linear structural systems, such as relatively rigid or non flexible
beams or shells.
[0097] This problem is also present for continuous system, be they
linear or non-linear. Linear continuous systems include complex
solid shapes, such as the ones modelled in usual CAD, e.g. massive
airplane fixtures, engine caps, cylinder blocks, crankshafts, or
other complex shape massive mechanical parts. Non-linear continuous
systems include massive deformable objects such as rubber silent
blocks used in the automotive industry or joints used for
windsurfing masts. The displacement vector for each finite element
only has translational components--although the modelled deformed
shape of the macroscopic object will possibly include rotation as a
whole.
[0098] Solutions exist in the prior art, for linear systems. Such
solutions are disclosed in Analyse des structures par elements
finis, J. F. Imbert, CEPADUES Editions, 1979-1984, page 276 or in
MSC NASTRAN Version 64, Handbook for linear Analysis, 1985, The
MacNeal-Schwendel Corp., page 4.4.1, Chapter 4.4, Multipoint
Constraint Operations. These documents suggest factorising the
constraints, that is expressing some of the constraints or degrees
of freedom as a function of the other degrees of freedom.
Specifically, let X be the column vector of degrees of freedom in a
given linear system. The generic problem to be solved may be
written as K.X+A.sup.T..lamda.=F.sub.app and A.X=b (1) where K is
the stiffness matrix, A is the Jacobian matrix of the displacement
constraints, .lamda. is the column vector of Lagrange multipliers,
b is the column vector of constraint right-hand sides, and
F.sub.app is the column vector of forces applied to the system
(A.sup.T denotes the transpose of matrix A). The prior art solution
consists in factorising the constraints, that is finding an
independent subset X.sub.a of degrees of freedom, along with
matrices N and M such that X=N.X.sub.a.+M (2)
[0099] This makes it possible to replace the value of X in eq. (1).
Assuming A.N=0, eq. (1) when pre-multiplied by N.sup.T provides
N.sup.T.K.N.X.sub.a=N.sup.T.(F.sub.app-K.M) which may be solved in
the independent degrees X.sub.a of freedom using prior art
solution. Once X.sub.a. (which includes the set of relaxed degrees
of freedom), is found, the full displacement may be recovered
thanks to eq. (2).
[0100] Besides the fact that the automatic partition of the set of
degrees of freedom into independent and dependent subsets is an
expensive process even in the case of linear problems, this
solution causes some additional problems in the case of non-linear
problems. First, in a non-linear system, the frames of reference
need to be updated, as explained above. This means that the
constraints will have to be factorised again each time the frames
of reference are updated; in other words, the value of matrix N in
eq. (2) needs to be computed again after each incremental step.
Second, in a computer system, where numbers are stored on a limited
number of digits, inverting a matrix, as required for factorising
the constraints, may lower the robustness of the computations.
Last, the stiffness matrix K would vary in each incremental step,
which requires repeatedly computing N.sup.T.K.N; this proves to be
costly.
[0101] In a linear system, co-rotational update is not needed; the
invention still avoids the need for factorising the constraints,
and only requires a simple change of frame of reference.
[0102] Thus, there is a requirement for allowing constraints to be
relaxed in a non-linear system, involving update of reference
frames. The invention suggests using hybrid reference frames for
minimizing energy, with local reference frames for those nodes of
the system with constraints consisting of partly relaxed and partly
imposed degrees of freedom.
[0103] Local reference frames are selected so that the released
degrees of freedom are expressed along the axes of the local
reference frame. Finding the subset of unconstrained degrees of
freedom X.sub.a may then be simply carried out by selecting in
equation (1) the lines that correspond to relaxed (non-imposed)
degrees of freedom (in other words, by eliminating from the
stiffness matrix the lines and columns corresponding to
imposed,--i.e.: non-relaxed-degrees of freedom).
[0104] This provides a solution to the problems discussed above.
The solution of the invention only requires the carrying out of
specifically selected changes of reference frames at each
incremental step; however, a change of reference frame is anyhow
necessary at each step for minimizing the energy; as discussed
below, the change to a set of hybrid reference frames is not more
complex from the computation standpoint.
[0105] The solution is exemplified in reference to FIG. 6, which
shows two finite elements 35 and 40, with their respective nodes
36, 41 and 41, 46. FIG. 6 further shows elementary reference frames
37 and 42 for the first and second finite elements 35 and 40;
frames 37 and 42 are respectively noted R.sub.1 and R.sub.2, for
the sake of simplicity. Let X.sub.1 be the column vector of degrees
of freedom at point 36, while X.sub.2 respectively X.sub.3 are the
column vectors of degrees of freedom at points 41 and 46.
[0106] In the absence of relaxed degrees of freedom, the stiffness
matrix would first be computed for each finite element 35 and 40,
in the elementary reference frame or in any frame relative to the
finite element. Once the stiffness matrix is computed for each
element in the corresponding elementary frame, the various matrices
are expressed in a global frame of reference, and may then be
summed, as explained in eq. (3) to (6) below.
[0107] Specifically, the equilibrium equation K.X=F is written, for
the first element 35 in the local frame R.sub.1 [ K 11 K 12 K 21 K
22 ] .function. [ X 1 X 2 ] 1 = [ F 1 F 2 ] ( 3 ) ##EQU5##
Similarly, for the second finite element, in the corresponding
local frame R.sub.2 for the second element, one obtains [ K 22 K 23
K 32 K 33 ] .function. [ X 2 X 3 ] 2 = [ F 2 F 3 ] ( 4 ) ##EQU6##
where each K.sub.ij is actually a 6.times.6 matrix if each X.sub.k
has six degrees of freedom. In eq. (3) and (4), the subscripts 1
and 2 to the matrices are representative of the fact that the
matrices are expressed in the local frames of reference R.sub.1 and
R.sub.2.
[0108] If no constraints are relaxed, one would simply apply a
change of reference frame to eq. (3) and eq. (4), so that they all
read in a global reference frame R.sub.G, referenced 82 in FIG. 6.
Assuming the matrix for passing from local reference frame R.sub.1
into global reference frame R.sub.G is noted P.sub.1, one has [ X 1
X 2 ] 1 = [ P 1 P 1 ] .function. [ X 1 X 2 ] G ##EQU7## where the
subscript G indicates that the degrees of freedom are expressed in
the global reference frame R.sub.G. Equation (3) thus changes into
[ P 1 P 1 ] T .function. [ K 11 K 12 K 21 K 22 ] .function. [ P 1 P
1 ] .function. [ X 1 X 2 ] G = [ P 1 P 1 ] T .function. [ F 1 F 2 ]
1 ##EQU8## where P.sub.1 is the change of reference frame matrix
for passing from local reference frame R.sub.1 into global
reference frame R.sub.G. Eq. (3) thus reads in frame R.sub.G [ P 1
T .times. K 11 .times. P 1 P 1 T .times. K 12 .times. P 1 P 1 T
.times. K 21 .times. P 1 P 1 T .times. K 22 .times. P 1 ]
.function. [ X 1 X 2 ] G = [ P 1 T .times. F 1 P 1 T .times. F 2 ]
( 5 ) ##EQU9## while eq. (4) reads in frame R.sub.G [ P 2 T .times.
K 22 .times. P 2 P 2 T .times. K 32 .times. P 2 P 2 T .times. K 32
.times. P 2 P 2 T .times. K 33 .times. P 2 ] .function. [ X 2 X 3 ]
G = [ P 2 T .times. F 2 P 2 T .times. F 3 ] ( 6 ) ##EQU10## P.sub.2
being the change of reference frame matrix for passing from local
reference frame R.sub.2 into global reference frame R.sub.G. One
obtains the global equilibrium equation by summing equations (5)
and (6) for all finite elements. This provides, in the example of
FIG. 6 with only two finite elements [ P 1 T .times. K 11 .times. P
1 P 1 T .times. K 12 .times. P 1 P 1 T .times. K 21 .times. P 1 P 1
T .times. K 22 .times. P 1 + P 2 T .times. K 22 .times. P 2 P 2 T
.times. K 32 .times. P 2 P 2 T .times. K 32 .times. P 2 P 2 T
.times. K 33 .times. P 2 ] .function. [ X 1 X 2 X 3 ] G = [ .times.
P 1 T .times. F 1 P 1 T .times. F 2 + P 2 T .times. F 2 P 2 T
.times. F 3 ] ( 6 ) ##EQU11## There is some ambiguity in these
formula--for instance, F.sub.2 denotes the force applied to the
second node 41 expressed in R.sub.1 and at the same time the force
applied to the second node 41 expressed in R.sub.2. The ambiguity
is raised thanks to the use of matrices P.sub.1 and P.sub.2; thus,
for the sake of simplicity, one does not use additional subscripts
for avoiding the ambiguity. Equation (6) shows the form of the
stiffness matrix , as computed in a global reference frame. The
example of FIG. 6 generalises to more than two finite
elements--with a summation of the stiffness matrices computed in
each local frame of reference. Due to the necessary co-rotational
update of the frames of reference, the computation of the changes
of frames P.sub.i for each element of the flexible system is
carried out at each incremental step.
[0109] The invention, in this embodiment, suggests using a hybrid
frame of reference for minimizing energy; for nodes with some
imposed and some released degrees of freedom, one uses a local
frame instead of using a global frame. The local frame is selected
so that the degrees of freedom are expressed along the axes of the
local frame.
[0110] Assume for instance that at least a degree of freedom is
released in node 41 which is otherwise fixed. Let R.sub.L be a
frame of reference where the released degree of freedom is
expressed along one axis; this frame is referenced 83 in FIG. 6.
For instance, if the relaxed constraint or degree of freedom is a
tangential one, whereby the system is free to slide in one
direction, one could use the tangential nodal frame as local frame
of reference R.sub.L. The same nodal frame could be used if the
system were free to slide along one direction and rotation around
the same direction. In such a case, for matrix X.sub.1, one would
use the change of frame P.sub.1 to the same global reference frame
R.sub.G; similarly, for matrix X.sub.3, one would use the change of
frame P.sub.2 to the same global reference frame R.sub.G. However,
for matrix X.sub.2, the frame is changed to the local frame
R.sub.L. eq. (3) above changes into [ P 1 P 1 .times. L ] T
.function. [ K 11 K 12 K 21 K 22 ] .function. [ P 1 P 1 .times. L ]
.function. [ X 1 .times. G X 2 .times. L ] = [ P 1 P 1 .times. L ]
T .function. [ F 1 F 2 ] 1 ##EQU12## where the last subscript,
where applicable, denotes the reference frame used for expressing
X.sub.1 and X.sub.2 and where [0111] P.sub.1L denotes the matrix
for changing frames of reference, from local frame of reference
R.sub.1 to local frame of reference R.sub.L and [0112] P.sub.2L
denotes the matrix for changing frames of reference, from local
frame of reference R.sub.2 to local frame of reference R.sub.L. One
understands that P.sub.1L could simply be the identity matrix, if
the local frame R.sub.L was the local frame R.sub.1. This may be
the case if frame R.sub.1 is adapted for expressing the degree(s)
of freedom along one axis.
[0113] Thus, in the hybrid frame of reference, eq. (3) and (4) read
respectively [ P 1 T .times. K 11 .times. P 1 P 1 T .times. K 12
.times. P 1 .times. L P 1 .times. L T .times. K 21 .times. P 1 P 1
.times. L T .times. K 22 .times. P 1 .times. L ] .function. [ X 1
.times. G X 2 .times. L ] = [ P 1 T .times. F 1 P 1 .times. L T
.times. F 2 ] .times. .times. and ( 7 ) [ P 2 .times. L T .times. K
22 .times. P 2 P 2 .times. L T .times. K 32 .times. P 2 P 2 T
.times. K 32 .times. P 2 .times. L P 2 T .times. K 33 .times. P 2 ]
.function. [ X 2 .times. l X 3 .times. G ] = [ P 2 .times. L T
.times. F 2 P 2 T .times. F 3 ] ( 8 ) ##EQU13## This leads to the
following stiffness matrix for the two finite elements [ P 1 T
.times. K 11 .times. P 1 P 1 T .times. K 12 .times. P 1 .times. L P
1 .times. L T .times. K 21 .times. P 1 P 1 T .times. K 22 .times. P
1 .times. L + P 2 .times. L T .times. K 22 .times. P 2 .times. L P
2 .times. L T .times. K 32 .times. P 2 P 2 T .times. K 32 .times. P
2 .times. L P 2 T .times. K 33 .times. P 2 ] .function. [ X 1
.times. G X 2 .times. L X 3 .times. G ] = .function. [ P 1 T
.times. F 1 P 1 .times. L T .times. F 2 + P 2 .times. L T .times. F
2 P 2 T .times. F 3 ] ( 8 ) ##EQU14##
[0114] eq. (8) demonstrates that the stiffness matrix and the
forces applied to the system may be expressed in a hybrid frame of
reference, with a computation which is substantially the same as
the one used in the case where constraints are not relaxed. The
only difference is that one has to find a local reference frame
R.sub.L and to compute two additional matrices P.sub.1L and
P.sub.2L for passing from the local reference frames R.sub.1 and
R.sub.2 into the new local reference frame at partly constrained
node(s) where a constraint is relaxed.
[0115] However, once eq. (8) is computed, minimizing energy does
not require anything like the computation of matrices N and M of
eq. (2) given above. Indeed, the local reference frame R.sub.L is
selected so that the relaxed degree of freedom is expressed along
one axis of the reference frame. This means that the selection of
the subset X.sub.a may be carried out simply by selecting all
constraints that are relaxed. In other words, one simply deletes in
eq. (8) the lines and columns of the stiffness matrix and the lines
of the load vector corresponding to the axis or axes of the local
frame of reference on which the constraint is not relaxed.
[0116] Specifically, in the example of FIG. 6, one considers for
each node six degrees of freedom, three in translation (u, v, w)
and three in rotation (.theta..sub.x, .theta..sub.y,
.theta..sub.z). Each of the X.sub.i is written as X i = [ u i v i w
i .theta. xi .theta. yi .theta. zi ] ##EQU15## The stiffness matrix
in eq. (8), expressed in the hybrid frame of reference, is thus
24.times.24 matrix. Assume that there is one degree of freedom at
node 41, the node being free to rotate around axis u--that is,
.theta..sub.x2 is relaxed. This means that u.sub.2, v.sub.2,
w.sub.2, .theta..sub.y2 and .theta..sub.z2 are fixed, while
.theta..sub.x2 and all variables in X.sub.1 and X.sub.3 need to be
computed. To this effect one proceeds as disclosed above in
reference to FIGS. 2 and 3 for the computation of displacement in
the nodes apart from the end nodes. More specifically: equation 8
is first detailed as: [ k 11 k 12 k 13 k 14 k 15 k 16 k 17 k 18 k
19 k 110 k 111 k 112 k 21 k 22 k 23 k 24 k 25 k 26 k 27 k 28 k 29 k
210 k 211 k 212 k 31 k 32 k 33 k 34 k 35 k 36 k 37 k 38 k 39 k 310
k 311 k 312 k 41 k 42 k 43 k 44 k 45 k 46 k 47 k 48 k 49 k 410 k
411 k 412 k 51 k 52 k 53 k 54 k 55 k 56 k 57 k 58 k 59 k 510 k 511
k 512 k 61 k 62 k 63 k 64 k 65 k 66 k 67 k 68 k 69 k 610 k 611 k
612 k 71 k 72 k 73 k 74 k 75 k 76 k 77 k 78 k 79 k 710 k 711 k 712
k 713 k 714 k 715 k 716 k 717 k 718 k 81 k 82 k 83 k 84 k 85 k 86 k
87 k 88 k 89 k 810 k 811 k 812 k 813 k 814 k 815 k 816 k 817 k 818
k 91 k 92 k 93 k 94 k 95 k 96 k 97 k 98 k 99 k 910 k 911 k 912 k
913 k 914 k 915 k 916 k 917 k 918 k 101 k 102 k 103 k 104 k 105 k
106 k 107 k 108 k 109 k 1010 k 1011 k 1012 k 1013 k 1014 k 1015 k
1016 k 1017 k 1018 k 111 k 112 k 113 k 114 k 115 k 116 k 117 k 118
k 119 k 1110 k 1111 k 1112 k 1113 k 1114 k 1115 k 1116 k 1117 k
1118 k 121 k 122 k 123 k 124 k 125 k 126 k 127 k 128 k 129 k 1210 k
1211 k 1212 k 1213 k 1214 k 1215 k 1216 k 1217 k 1218 k 137 k 138 k
139 k 1310 k 1311 k 1312 k 1313 k 1314 k 1315 k 1316 k 1317 k 1318
k 147 k 148 k 149 k 1410 k 1411 k 1412 k 1413 k 1414 k 1415 k 1416
k 1417 k 1418 k 157 k 158 k 159 k 1510 k 1511 k 1512 k 1513 k 1514
k 1515 k 1516 k 1517 k 1518 k 167 k 168 k 169 k 1610 k 1611 k 1612
k 1613 k 1614 k 1615 k 1616 k 1617 k 1618 k 177 k 178 k 179 k 1710
k 1711 k 1712 k 1713 k 1714 k 1715 k 1716 k 1717 k 1718 k 187 k 188
k 189 k 1810 k 1811 k 1812 k 1813 k 1814 k 1815 k 1816 k 1817 k
1818 ] .function. [ u 1 .times. G v 1 .times. G w 1 .times. G
.theta. x .times. .times. 1 .times. G .theta. y .times. .times. 1
.times. G .theta. z .times. .times. 1 .times. G u 2 .times. L v 2
.times. L w 2 .times. L .theta. x .times. .times. 2 .times. L
.theta. y .times. .times. 2 .times. L .theta. z .times. .times. 2
.times. L u 3 .times. G v 3 .times. G w 3 .times. G .theta. x
.times. .times. 3 .times. G .theta. y .times. .times. 3 .times. G
.theta. z .times. .times. 3 .times. G ] .function. [ F x .times.
.times. 1 .times. G app F y .times. .times. 1 .times. G app F z
.times. .times. 1 .times. G app M x .times. .times. 1 .times. G app
M y .times. .times. 1 .times. G app M z .times. .times. 1 .times. G
app F x .times. .times. 2 .times. L app F y .times. .times. 2
.times. L app F z .times. .times. 2 .times. L app M x .times.
.times. 2 .times. L app M y .times. .times. 2 .times. L app M z
.times. .times. 2 .times. L app F x .times. .times. 3 .times. G app
F y .times. .times. 3 .times. G app F z .times. .times. 3 .times. G
app M x .times. .times. 3 .times. G app M y .times. .times. 3
.times. G app M z .times. .times. 3 .times. G app ] ##EQU16##
[0117] In this equation, k.sub.ij, with 1.ltoreq.i.ltoreq.18 and
1.ltoreq.j.ltoreq.18 are the various coefficients of the stiffness
matrix, computed as discussed above. The first subscript to u, v,
w, .theta..sub.x, .theta..sub.y and .theta..sub.z is representative
of the node; the second subscript to u, v, w, .theta..sub.x,
.theta..sub.y and .theta..sub.z is representative of the frame in
which the displacement is expressed, G standing for the global
frame of reference R.sub.G and L standing for the local frame of
reference R.sub.L.
[0118] As regards the applied load, F.sub.x.sup.app,
F.sub.y.sup.app and F.sub.z.sup.app are representing the components
of the force applied to the node, while M.sub.x.sup.app,
M.sub.y.sup.app and M.sub.z.sup.app are representing the components
of the torque applied to the node; The first subscript (Arabic
subscript) is representative of the node to which the load is
applied; the second subscript is representative of the frame in
which the displacement is expressed, G standing for the global
frame of reference R.sub.G and L standing for the local frame of
reference R.sub.L.
[0119] When considered in conjunction with the five constraint
conditions: { u 2 .times. L = u _ 2 .times. L v 2 .times. L = v _ 2
.times. L w 2 .times. l = w _ 2 .times. L .theta. y .times. .times.
2 .times. L = .theta. _ y2L .theta. z .times. .times. 2 .times. L =
.theta. _ z .times. .times. 2 .times. L ##EQU17## expressed in the
local coordinate system R.sub.L, the equation becomes: [ k 11 k 12
k 13 k 14 k 15 k 16 k 110 k 21 k 22 k 23 k 24 k 25 k 26 k 210 k 31
k 32 k 33 k 34 k 35 k 36 k 310 k 41 k 42 k 43 k 44 k 45 k 46 k 410
k 51 k 52 k 53 k 54 k 55 k 56 k 510 k 61 k 62 k 63 k 64 k 65 k 66 k
610 k 101 k 102 k 103 k 104 k 105 k 106 k 1010 k 1013 k 1014 k 1015
k 1016 k 1017 k 1018 k 1310 k 1313 k 1314 k 1315 k 1316 k 1317 k
1318 k 1410 k 1413 k 1414 k 1415 k 1416 k 1417 k 1418 k 1510 k 1513
k 1514 k 1515 k 1516 k 1517 k 1518 k 1610 k 1613 k 1614 k 1615 k
1616 k 1617 k 1618 k 1710 k 1713 k 1714 k 1715 k 1716 k 1717 k 1718
k 1810 k 1813 k 1814 k 1815 k 1816 k 1817 k 1818 ] .function. [ u 1
.times. G v 1 .times. G w 1 .times. G .theta. x .times. .times. 1
.times. G .theta. y .times. .times. 1 .times. G .theta. z .times.
.times. 1 .times. G .theta. x .times. .times. 2 .times. L u 3
.times. G v 3 .times. G w 3 .times. G .theta. x .times. .times. 3
.times. G .theta. y .times. .times. 3 .times. G .theta. z .times.
.times. 3 .times. G ] .function. [ F x .times. .times. 1 .times. G
TOT F y .times. .times. 1 .times. G TOT F z .times. .times. 1
.times. G TOT M x .times. .times. 1 .times. G TOT M y .times.
.times. 1 .times. G TOT M z .times. .times. 1 .times. G TOT M x
.times. .times. 2 .times. L TOT F x .times. .times. 3 .times. G TOT
F y .times. .times. 3 .times. G TOT F z .times. .times. 3 .times. G
TOT M x .times. .times. 3 .times. G TOT M y .times. .times. 3
.times. G TOT M z .times. .times. 3 .times. G TOT ] ##EQU18##
[0120] Once energy is minimized (thus obtaining the X.sub.a set of
degrees of freedom, which includes the relaxed rotation about the
local u axis at the middle node, that is, the local degree of
freedom .theta..sub.x2, along with all the global degrees of
freedom of the other two nodes), the full vector of degrees of
freedom X may be recovered as follows: [ u 1 .times. G v 1 .times.
G w 1 .times. G .theta. x .times. .times. 1 .times. G .theta. y
.times. .times. 1 .times. G .theta. z .times. .times. 1 .times. G u
2 .times. L v 2 .times. L w 2 .times. L .theta. x .times. .times. 2
.times. L .theta. y .times. .times. 2 .times. L .theta. z .times.
.times. 2 .times. L u 3 .times. G v 3 .times. G w 3 .times. G
.theta. x .times. .times. 3 .times. G .theta. y .times. .times. 3
.times. G .theta. z .times. .times. 3 .times. G ] .function. [ u 1
.times. G v 1 .times. G w 1 .times. G .theta. x .times. .times. 1
.times. G .theta. y .times. .times. 1 .times. G .theta. z .times.
.times. 1 .times. G 0 0 0 .theta. x .times. .times. 2 .times. L 0 0
u 3 .times. G v 3 .times. G w 3 .times. G .theta. x .times. .times.
3 .times. G .theta. y .times. .times. 3 .times. G .theta. z .times.
.times. 3 .times. G ] .function. [ 0 0 0 0 0 0 u _ 2 .times. L v _
2 .times. L w _ 2 .times. L 0 .theta. _ y .times. .times. 2 .times.
L .theta. _ z .times. .times. 2 .times. L 0 0 0 0 0 0 ]
##EQU19##
[0121] Notice that in this restitution formula the displacement
vector components are expressed partly in the local reference
frame, that is, at the middle node, and partly in the global
reference frame, that is, at the other two nodes; it is of course
now possible to entirely express the displacement vector in the
global reference frame.
[0122] Once all components of the X.sub.a vector are computed,
co-rotational update as well as update of the nodal frames of
reference may be carried out as discussed above, if needed. The
only difference is that the update uses, for nodes having a degree
of freedom, a combination of the computed components and of the
constrained components. In the example of FIG. 6, u.sub.2, v.sub.2,
w.sub.2, .theta..sub.y2 and .theta..sub.z2 are fixed, while
.theta..sub.x2 is computed at the same time as the components of
X.sub.1 and X.sub.2. For the update of the reference frame and of
the nodal frame, one may proceed as disclosed above, using
displacement components (u.sub.2, v.sub.2, w.sub.2, .theta..sub.x2,
.theta..sub.y2, .theta..sub.z2), without taking into account the
fact that these components result from [0123] a fraction of the
predetermined displacement boundary conditions; [0124] the
minimizing of energy during the incremental step. If update is not
needed, e.g. in the case of linear systems or in the case of
continuum systems, one may proceed directly to the next step.
[0125] The invention makes it possible to relax some degrees of
freedom at constrained nodes of the finite element representation
of the system. It uses the fact that a careful selection of the
frames in which the energy is minimized makes it possible to avoid
any factorisation of constraints.
[0126] Examples of the process are given below, in reference to
FIG. 7 to 22. These four examples use a non-linear solver core
(FORTRAN code) integrated in CATIA V5, of Dassault Systemes,
Suresnes France. For each example the simulation consists of 75
incremental steps (45 steps from start to initial configuration, 30
steps from initial to final configuration). Only a few of these
steps are displayed in the figures. All runs were performed on an
IBM Thinkpad T42 (1.8 GHz Pentium IV processor, 1.0 GB of RAM).
[0127] The first example is represented in FIGS. 7 to 10. The
system is a single segment; both end points are fully imposed. The
system is analysed as a 12 elements model. Full run in carried out
in 0.28 sec CPU time; in each incremental step of the
initial-to-final phase, convergence is attained in 6 iterations on
the average.
[0128] Node positions and nodal reference frames are represented in
all figures. FIG. 7 shows the start configuration, in which the
shape of the system is fully specified by the user.
[0129] FIG. 8 shows the initial configuration; end point positions
and orientations are specified by the user, and represented in the
figure as light reference frames. Deformed shape of the system is
computed by the program.
[0130] FIG. 9 shows the final configuration; end point positions
and orientations are specified by the user and are represented as
thick reference frames. Deformed shape is computed by the program.
The difference between the end points positions and orientations in
FIGS. 8 and 9 represents the predetermined displacement boundary
conditions to which the system is submitted.
[0131] FIG. 10 shows three intermediate configurations; end point
positions and orientations are interpolated, and represent a
fraction of the predetermined displacement boundary conditions.
Deformed shape is computed by the program.
[0132] The second example is represented in FIGS. 11 to 14. It is
similar to the first example, except for the fact that the rotation
degree of freedom (around a local axis orthogonal to the
deformation plane) is released at the left end point, both in the
initial and in the final configurations.
[0133] Full run is carried out in 0.27 sec CPU time; as in the
first example, convergence is attained in 6 iterations per
increment on the average.
[0134] FIGS. 11-14 respectively correspond to FIGS. 7-10, and are
not described again. As for the first example, node positions and
nodal reference frames are represented in all figures.
[0135] The third example is represented in FIGS. 15 to 18. In this
example, the system comprises three segments of different
characteristics; one end point and one intermediate point fully are
fully imposed. This appears on the figure, where one clearly sees
that the slenderer section of the system is not constrained and
bends downwards under the action of gravity.
[0136] Finite element analysis is carried out with a 3.times.12
elements model.
[0137] Full run is carried out in 6.45 sec CPU time; convergence is
attained in 40 iterations per increment on the average.
[0138] FIGS. 15-18 respectively correspond to FIGS. 7-10, and as
not described again. Node positions and nodal reference frames are
not represented in all figures.
[0139] The fourth example is represented in FIGS. 19 to 22; this
example is identical to the third example, except for the fact that
all three rotation degrees of freedom are released at the
intermediate point.
[0140] Full run is carried out in 5.2 sec CPU time; in each
incremental step, convergence is attained in 25 iterations per
increment on the average.
[0141] FIGS. 19-22 respectively correspond to FIGS. 15-18, and are
not described again. Note the difference in the final
configuration, due to the different constraints applied to the
system.
[0142] The process is not limited to the examples provided above.
For instance, in the example of FIG. 1, the predetermined
displacement is considered at both ends of the flexible system. One
may also consider predetermined displacement--fully constrained or
not--at more than two points in the system.
[0143] In the process disclosed in reference to FIGS. 1 to 3, the
co-rotational update is carried out after each incremental step--or
before the next incremental step. It is also possible not to carry
out the co-rotational update at each step; specifically, one may
only carry out the co-rotational update when necessary, that is
when the change of angle between each elementary frame of reference
and the corresponding updated elementary frame of reference is
higher than a predetermined value. Practically speaking, this
predetermined value may be between 0 and 30 degrees with low
values, e.g.: 1 sec, corresponding to quasi-systematic
updating.
[0144] A scaling factor is disclosed above for controlling
convergence during an incremental step; such a scaling factor could
also be used for the purposes of controlling the behaviour of a
flexible system. For instance, applying a value of zero to the
scaling factor at one given incremental step results in simulating
a visco-elastic flexible system: all initial stresses are
relaxed.
[0145] The various embodiments disclosed above may be combined, for
a more efficient simulation of the flexible system. One may also
carry out the various embodiments independently. For instance, the
stopping criterion may be used independently of nodal frames. The
stopping criterion discussed above may also be used independently
of the co-rotational update discussed in reference to FIGS. 3 and
4; the stopping criterion could also be used for non-structural
(i.e.: continuous) elements; such non-structural or continuous
elements are actually "solid" elements, which do not possess nodal
rotation degrees of freedom.
[0146] As discussed above, the use of hybrid reference frames
disclosed in reference to FIG. 6 may be used independently of the
co-rotational update discussed previously. One could also use such
hybrid reference frames in a linear solver.
* * * * *