Finite Element Modelling Group
This page uses CSS to present the content in the best possible manner.
If you can see this message, then CSS is not enabled in your browser,
and the page will not appear as the designer intended.
Finite Element Modelling Group
Introduction
The areas of application and the potential of the finite
element method are enormous. Since it's inception in the 1950's, it has
developed to be the most powerful tool used in the solution of partial
differential equations, which occur in almost all of the applied sciences
and engineering fields. It's development has not been paralleled by any
other numerical analysis procedure and it has made many other numerical
analysis techniques and experimental testing methods redundant. However,
despite the time and money invested in an effort to mature the method,
there are many avenues to be explored before the method's true elegance
and utility is exploited.
The time taken to achieve equilibrium for any given analysis
is a function of the number of degrees of freedom (DOF) of the model being
analysed. Therefore, if the DOF used to model a particular problem can
be reduced while maintaining the same accuracy of the results, an important
step in the development of the method has been achieved. Any steps that
involve converting the complex real object geometry & loading to a
more efficient finite element model (and hence reducing the number of DOF
to be solved for) are known as idealisation.
NAFEMS [1] suggests the following categories of idealisation:
MATHEMATICAL MODEL SELECTION
Behavioural assumptions introduced in order to yield
a solvable set of expressions but having the side effect of introducing
approximations. An example of this is assuming linear elastic behaviour
in the finite element method.
PHYSICAL SCALES
Behavioural considerations at different size scales but
critical to the prediction of behaviour at all scales. For example, composite
materials and soils, where the micro-mechanical behaviour of the individual
constituents determines the overall behaviour of the material.
SYMMETRY
By taking advantage of the reflective, rotational or
translational symmetry within a component, only part of the component need
be modelled, so long as appropriate constraints are applied along the symmetry
lines.
DOMAIN SPECIFICATION OR SPATIAL EXTENT
The domain over which the mathematical model is applied.
Domain idealisation focuses on a particular area of the model to be analysed
and establishes the boundary condition locations e.g. sub-modelling. Figure
1 below shows an example of an analysis which employs sub-modelling.
As in Figure 3, a large percentage of finite element analyses
can make use of reduced dimensional element types that are defined in terms
of a reduced geometric representation with associated element properties
that account for the dimensions not included. These element types tend
to produce more computationally efficient models, thus reducing analysis
time and cost.
Because of the nature of the geometry in Figure 3, the complete model could be dimensionally reduced to 1D beam elements that contain all the geometric properties of the 3D model. However this case does not occur very often when modelling engineering components. In most practical models, there are indeed areas of the model which are ideal candidates for dimensional reduction, but they are usually bounded by areas that contain complex geometry. Therefore, to dimensionally reduce these regions of constant geometry, a method to couple or join the newly created reduced sections to the original complex sections of the model is needed.
The process of joining entities of different geometric dimension is known as multi (or mixed) dimensional coupling. Figure 4 illustrates the concept, and it is seen that there is a marked reduction in the number of elements required to carry out the analysis (each line in (ii) consists of 1 element).
Mc Cune [], [], concentrated on 3D-2D and 2D-1D coupling. The logical next step in the process is to reduce regions of constant cross-section in 3D solids to their 1D equivalent beam element. In order to analyse stress concentrations around discontinuities in geometry and loading, the ends of this 1D beam must then be coupled to the new surfaces created by the dimensional reduction process, Figure 5. By doing so, a marked reduction in the number of DOF required to analyse the model can be achieved without any loss of accuracy.
In this paper the coupling of a 1D beam element to a 3D
continuum is described. Alhough the procedure has only been implemented
for linear elements, it illustrates the concept and may be further developed
into a general and robust process. The reader is introduced to the area
by firstly dealing with the simple case for axial force in section 3. This
will be followed by the bending moment and torsion cases in sections 4
& 5 respectively. In section 6 the shear force case will be discussed,
and section 7 provides some model examples and the results from performing
the coupling of the mixed dimensions. Finally sections 8 and 9 conclude
the work presented and indicate the plans for further development.
| Beam Cross-Sectional Area | A |
| Beam Displacements | u, v, w |
| Beam Forces | Fx, Fy, Fz |
| Beam Rotations | x, y, z |
| Bending moments | Mx, My, Mz |
| Centroid of element face | xc, yc |
| Continuum Displacements | U, V, W |
| Moment of Inertia | I |
| Poisson's Ratio | |
| Shape functions | N |
| Shear Modulus | G |
| Stress (direct) | x, y, z |
| Stress (shear) | xy, xz, yz |
| Stress function | |
| Temperature | t |
| Young's Modulus | E |
Axial Force
The following derivation is based on an outcome of Reissner's transverse plate bending theory [] by which simple relationships are shown to be obtainable between the translational and rotational degrees of freedom in a plate and the displacements in the 3D continuum. To couple a beam and a solid, the first step is to equate the work done on either side of the interface between dimensions.
Equating the work done by the axial force acting on the 1D beam with the work done by the surface stresses of the 3D body at the interface, the following equation results:
... 3-1If the 3D region is long and slender, then the axial stress is uniform over the cross-section and is given by:
... 3-2In the 3D model, the axial displacement at any point, in terms of the nodal displacements{W} and shape functions[N], is:
This implies:
... 3-4Since this must be true for any axial force, the beam displacement w and the displacements of the 3D continuum nodes on the interface {W} must be related by:
... 3-5Displacement compatibility between the 1D beam element and the adjacent 3D continuum elements can therefore be enforced as a multipoint constraint equation of the form
... 3-6This can be applied in the ABAQUS commercial finite element package as a *EQUATION command [].
Bending Moment
Again Reissner's relationships can be used to equate the work done on either side of the coupling interface. Considering the moment about the x-axis and equating the work done at the interface gives the following equation:
... 4-1From simple beam theory, assuming for simplicity that the co-ordinate system is aligned with the principal axes,
... 4-2Substituting also for the displacements in the z-direction in the 3D finite element model,
On substitution:
... 4-4This is the general solution for any type of element.
Consideration of the bending moment acting about the y-axis gives a similar equation:
... 4-5
Torsion
The distribution of shear stress on the cross-section of a beam subjected to a torsional moment is commonly solved by introducing a single stress function. If a function (x,y), the Prandtl stress function, is assumed to exist such that:
... 5-1; 
... 5-2then the stress function must satisfy the differential equation:
... 5-3where is the twist per unit length of the beam and G is the shear modulus. The stress function therefore must satisfy Poisson's equation []. A conductive heat transfer case also may be represented as a form of Poisson's equation.
Therefore, the variation of the Prandtl stress function over any cross-section can be found using the facilities available in standard finite element packages for conductive heat transfer. The shear stress on the cross-section can be inferred from the resulting temperature gradients. The cost of a 2 dimensional heat transfer analysis of the cross-section is much less than that of a 3D analysis of the stress.
The total torque generated by a given twist can be found from the torsion
analysis [9] as:
The coupling equation is again formed by equating the work done at the interface.
On the 3D side of the interface, the work done is:
... 5-5On the 1d side, the work done is:
... 5-6Therefore,
... 5-7Replacing the integral over the whole cross-section with the sum of the integrals over the element faces lying on the interface and the continuum displacements with the 3D finite element displacements:
... 5-8The total torque Mz in response to any arbitrary twist can be found from Equation 5-4, as can the shear stresses xz and yz on each element face. Evaluating the [C] and [D] matrices reduces to summing integrals of the form:
... 5-9The shear stress xz can be written in terms of shape functions [N] and nodal values of the stress function {} from the 2D analysis of the cross-section as:
... 5-10so that the integral over an element face becomes:
... 5-11where {} are the nodal temperatures found in the heat transfer analysis. The complete multipoint constraint equation can then be written in the form:
This equation can also be applied as a linear constraint equation in
the finite element model (e.g. applied in Abaqus as a *EQUATION command).
The effect of this equation is to couple the displacements of the 3D continuum
nodes on the interface to the twisting rotation of the beam node such that
the distribution of shear stress on the interface is the same as that given
by the St. Venant torsion analysis of the beam cross-section.
Shear Force
The approach in sections 3-5 may be extended to ensure the proper distribution of shear stress in response to a given shear force using an analysis of the cross-section similar to that for torsion []. The stresses on the cross-section at any point (x,y) due to a shear force Fx on the cross-section are given in terms of a stress function as
... 6-1; 
... 6-2where must obey the governing equation
... 6-3with the boundary condition
... 6-4on the boundary of the section.
Figure 8 illustrates the shear stress xz derived from the results of a 2D thermal analysis of the cross-section. Given the shear stresses on the cross-section, the appropriate multipoint constraint equation can be generated as before by equating the work done on both sides of the interface. This yields
... 6-5This procedure has been demonstrated successfully for 2D - 1D coupling of beams and plates [6].
Results
Sample results for the cases considered are given in their respective sections below. They demonstrate the concept graphically and give an indication of the possible benefits of carrying out the reduction & coupling.
AXIAL FORCE
The stress plot of the coupling for the Axial force case is shown in Figure 9 on the right. The stress is constant across the section as expected, and is the same as the analytical result.
BENDING MOMENT
The stress plots of the coupling for the bending moment cases are shown in Figure 10 and Figure 11 below. Again the results compare favourably with analytical calculations. Note that there is no disturbance to the stress contours at the interface between the 1D beam and the 3D solid elements.
TORSION
As can be seen from the contours of maximum shear stress in Figure 12 and Figure 13 the coupling equations produce stress contours which are axisymmetric and as close to the analytic solution for a thin walled tube as can be obtained with constant stress elements and nodal averaging.
Figure 14 and Figure 15 show the von Mises stress obtained for coupling a square bar to a beam element of equivalent properties. Again there is no disturbance to the stress contours at the dimensional interface and the stresses on the interface are effectively identical to those obtained by the St. Venant torsion analysis or at sections in the interior of the 3D model. There is some disturbance at the rigidly fixed end due to the method by which the fixing was defined.
Discussion
Using the procedure described here, the analysis of a slender region within a 3D solid can be reduced to a 1D analysis plus 6 simple analyses of the cross-section. These determine the stress distribution on the cross-section in response to each component of force and moment.
The cross-sectional analyses can also be used to generate 6 multipoint constraint equations linking the displacements and rotations of the beam element to the nodal displacements of the 3D continuum elements representing the material adjacent to the slender area of the model represented by the beam element. This greatly facilitates mixed-dimensional modelling, where 3D details such as joints, changes in section or loading are properly represented, but dimensionally-reduced beam elements are used to model slender parts (economically and accurately), where 3D elements are expensive and potentially ill-conditioned.
The coupling does not introduce spurious stresses at the interface, so that measures based on the implied stress jumps between the 3D and the reduced dimensional models can be used as an a posteriori estimate of the idealisation error introduced by dimensional reduction [5], [6]. The cross-sectional analyses can also be used to determine section properties for beam elements and to transfer beam force and moment results to stresses on the cross-section for visualisation.
The ideas presented here have already been applied to beam-plate-solid
coupling [5] and should be extensible to the coupling of arbitrary combinations
of beam, plate, shell and solid elements of any formulation.
Conclusions