The pin jointed bar is the simplest element, we will see how it is used in a stress analysis. The ends of the element are assumed to be pin jointed so it only transmits tensile or compressive axial forces. Referring to the diagram below it has two nodes i and j ; it lies in the x axis and therefore only experiences axial displacements. Assuming that the displacement u varies linearly along the bar (length- L) then

u =a+bx : where a and b are constants.

Diagram showing a) one-dimensional pin jointed bar b) equivalent finite element

0We can now use this one dimension element, we will demonstrate its use in a typical application, stress analysis.

When a structure is loaded it attains some equilibrium state, the potential energy of the system is at minimum. The potential energy of the system can be defined by :-

For the single element we are considering, the strain energy stored in element following any deformation (assuming constant cross sectional area, A) is

Using this one dimensional linear element, we will demonstrate a one dimensional eigenvalue problem, that of a one dimensional string.

Remember the characteristic eigenvalue is

[K] is the global stiffness matrix and [m] the mass matrix.

Lets consider the ‘string’ to be made of 4 linear elements :-

Remembering that the stiffness matrix was written for nodes ‘i’ and ‘j’, the stiffness matrix for the element with nodes ‘j’ and ‘k’ would be

The global stiffness matrix, is the sum of all the element stiffness matrix, with respect to their nodes . For this case it would be

This assuming that the all elements have the same , of course they don’t have to,

generally the term would be multiplied through, before assembling the global stiffness, matrix.

So what about the [m], the simplest form this can take is lumped mass, where the mass is lumped at the nodes, this makes the [m] matrix diagonal. There are other schemes to produce mass matrices such as consistent mass where the shape function is used to define the distribution of mass across the element.

Now we can solve this equation

Assuming =2 and

For non-trivial solutions we demand that :

This piece analysis has been implemented in the matlab file, feaex.m It assembles the stiffness matrix as shown in the text, and solves the eigenvalue equation, the eigenvectors are normalised, but they so compare with the values calculated here.

Below are plots of the first 3 eigenmodes calculated above.

We are going to get a better approximation if we use more elements, lets run the program again, but with 8 elements.

As the matlab program shows the effect of using a greater number of elements and for interpolation function. The cost of using interpolation functions for greater accuracy, is the computation time needed. Sometimes using more elements gives results to the required accuracy.

Important things to notice from this is the way each element of the same type has the same elemental equation, and therefore its on local stiffness matrix k, the coordinates will be different for each element. In this case the material properties for each element in the same, however in the tabla model, the mass varies from one element to another.

The idea’s of the previous section are extended in this section. Since the area of finite element analysis is so big, only the aspects which are relevant to the specific problem and the Abaqus program will be covered.

As already discussed the eigenvalues and eigenvectors are very important in vibrational problems. They could provide away of implementing an efficient model, of the tabla.

Once the a finite elements, are defined in the input file to the Abaqus program, there are a number of analysis’s which can be done. We did a eigenvalue analysis, which is defined by :-

If we assume sinusoidal vibration

The standard form of the eigenvalue equation is

By comparing [Eq.16] with the eigenvalue equation :

The elements used in the model have been defined as membrane having membrane properties, that is forces that act along member axis and tangential to element surfaces, it has strength but no resistance to bending. A rubber balloon of such a material with membrane properties.

This makes the problem non-linear, something we avoided by simplification in the state space model. If the structure has an applied load, the stiffness of the structure will change due to the membrane forces, so we have to take this into account. The geometric stiffness matrix(stiffness due to the membrane forces) depends on the stresses of the structure so therefore is function displacement, thus making the problem non-linear.

The problem now takes the form.

To overcome this potential problem the geometric stiffness is calculated by recalculating the stiffness matrix after adjusting the nodal co-ordinates with the calculated displacements. This is usually done using the Newton-Raphson method.

If the displacement is going to be small the load is applied in one step, over several iterations, this can been seen in the input file for the model.

Including membrane forces, the eigenvalue equation, becomes

If the material has no bending stiffness at all k=0.

Once the local matrices are assembled into the global stiffness matrix, the boundary conditions have to be added just as in the example.

These are entries into the {d} matrix.

In the model, we used isoparametric triangular elements, and isoparametric Lagrange quadratic elements, the derivation of these elements will be show.

The next section defines the principle nodes. Once these are defined, the *nfill, *ngen, and ncopy commands, are used to generate a 2d system of nodes.

Lets look at the basic *node command

-where n-node number, x,y,z - co-ordinates.

Here the node command is used with the nset command, which simply labels a node or group of nodes with an identifier for reference latterly in the file.

Defines node number 1, which has an identifier N1 at coordinates 0,0,0 (spaces in between comma’s = zero’s). The second line defines node 9 (N9) at coordinates 30.0,0,0

The units of length used are up to the user, but must be kept consistent through out the file, this is very important when it comes to densities and Young’s modulas etc. which are inputted later in file.

So in the first few line nodes 1, 9, 10, 30, 36, and 38 are defined, and are labelled with the identifier’s Nx . Also nodes 40, and 1000 are defined. Lets look at these.

Notice that nodes 9 and 10 (N9 and N10), have the same coordinates, they are overlapping and not connected together in any way(yet).

The next procedure is to fill in the nodes between N1, N9/N10, N30, N36, and N38.

This is done for the ‘linear/one dimensional’ with the *NFILL command.

Lets take the first *NFILL command -

This fills in the nodes between N1 and N9, with 8 nodes, and the node numbers are incremented by 1 each time. The new set of nodes form the node set R1. Finally the nodes are placed with a bias of 1.1, so they do not have a uniform distance between them. The TWOSTEP tag, ensures that only the bias is applied to every second node, this means even numbered nodes, will be in the centre, of adjacent nodes. This is very important when second order elements are used, which require extra nodes at the mid-point of their vertices.

The second *NFILL command, fills in the nodes between N10 and N30, with 20 nodes, with a bias(which affects ever other node) of 1.025. These new nodes form part of the set R1.

This generates the nodes between 30 and 36, the new nodes form part of the set R1.

As you will see from the diagram, nodes 30 and 36 form the endpoints of a 90deg arc,(the edge of the drum). The LINE=C indicates that nodes are to be generated in a circle, the point of rotation is node 40, and the nodes numbers are in increments of one. Finally the last three numbers indicate which plane(x,y,z) the circle is, here its in the x plane.

This simply generates the nodes between 36 and 38 (37) with uniform distribution, the nodes numbers are incremented by 1, and they are added to the set R1.

This command takes our node set R1 and rotates (indicated by the ‘SHIFT’) and copies it, 35 times, the resulting set of nodes is called RS. Each time the node set R1 is rotated the node numbers have 40 added to them. The first line of numbers indicate the value of shift, in this case none, the first, the second line defines two points one at (0,0,-1) and (0,0,1), this essentially define the Z plane as that of rotation, the last number is the angle of rotation (10 degrees), through which R1 is rotated.

Now we have rotated the nodes we are left with a large system of nodes, what remains is to define the elements. First we define, some new node sets, the importance of these new node sets will become apparent later on.

The defines the node set CONS1, which starts from node 9 and ends at node 1409, in increments in 40’s. The GEN tag forces the set to be defined by existing node numbers.

This similarly defines the node set CONS2, which starts from node 10 and ends at node 1410, in increments in 40’s.

These sets each form a concentric circle, around the origin and have the same radius as the x coordinate’s of node9(CONS1) and node10(CONS2). These nodes overlap at the coordinates (30,0,0). As stated before although the nodes overlap, they are not in any way fixed to each other, when we rotated and copied the node set R1, the resulting node numbers 9+40n and 10+40n when n is an integer, although overlapping, are not fixed either.

This defines the node set pull, which are the nodes along the edge of the radius. Later there will be a displacement on them, to tension the drum.

Here we have used two different geometry’s of element , the m3d6 and m3d9,

these are membrane elements , with 6 and 9 nodes respectively.

This defines the master elements, they are all of the m3d6 type, and they form the element set, cent, these are the centre elements. The first number is the element number, these are consecutive. The next six numbers are the nodes numbers which the element is associated with. Notice how they are in the ordered, following the order in the top diagram.

This defines the master elements which are of the m3d9 type. The format is the same as before, the first number is the element number and subsequent numbers are the nodal coordinates.

This fills in the elements around the centre, using element 19, the four sided element, as a master element. 17 elements are to be defined, the node numbers are incremented by 80 each time. The element number is to be incremented by 1 each time. There are 3 rows (radially, in this case), and the nodes numbers are incremented by 2 (outwards), the element numbers are increment by 18 per row.

Also element 36 is copied outwards. These all form the element set cent, and represent the sythi.

The diagram below should clarify, the explanation.

This once again defines the rest of the elements which represent the lao, area of the puri. The format is the same as before. Using 101 as the master element 17 elements are defined in each radius, with a node number increment of 80, and an element number increment of 1. 14 rows are also defined, the node numbers increment by 2 from row to row, and the element numbers are incremented by 18. Then element 118, is used to generate elements, outwards. These form the set, out.

This makes all of the elements from 281 to 334 (in increments of 1) an element set. In physical terms these are the elements which are in contact with the rigid surface (not yet defined).

This defines the element set, 'CENT' as having membrane properties, and having the material properties of the material 'CC'. CC is just an identifier the actual properties are declared later in the program.

The membrane section, has thickness of 1.8mm.

This type of property is used to represent thin surfaces, which have strength but no bending stiffness.

Once again this define the element set 'OUT', as have membranes properties which has material properties as the material 'RR'. Its thickness is 1.8mm.

This defines the material properties of the material RR. For a membrane section the properties needed are the density, Poisson’s ratio, and the young’s modulas.

The density is defined as 4.6E-9 (kg/mm cubed), and the young’s modulas as 1E4, the Poisson ratio, is 0.45.

Once again this defines the material properties of material CC.

The density= 4.6E-9 kg/mm cubed, young’s modulus=1E4, and the Poisson’s ratio=0.45.

This defines a rigid surface, the type will be a volume of revolution, its identifier is MSURF, the reference node will be node number 1000. The SMOOTH=1 defines radius of curvature, over which to smooth discontinuities between adjoining straight line and circular arc elements.

The first row of numbers define the global coordinates,(from which subsequent local coordinates will be relative to) and the axis of rotation.

This defines the direction of the contact surface between the element set and whatever its touching(in this case the rigid surface). The element set is CNT, and the direction is SNEG - which means its a negative outward normal. The other possibility is SPOS (positive outward normal).

This defines the to surfaces ssurf and msurf to be in contact with each other. The interaction (friction) between are defined by the identifier MU.

This define the fricative constant between for a particular contact pair, in this case its the only contact pair, MU.

This defines a multi-point constraint, between the node sets CONS1 and CONS2

this effectively fixes the inner elements (type m3d6) with the outer elements(type m3d9). The 'TIE', makes all the degrees of freedom equal at each node.

This command allows a arbitrary time variations of the load displacements.

Here 'time', is step time, and the values given are relative magnitude, the displacement has an identifier 'ZDISP'.

This opens/creates a restart file, its updated at every increment of the analysis. All the results are stored here.

This is the first step in the analysis, the loading given by the amplitude command is applied, instantly at the beginning of the step. The NLGEOM tag, indicates that non-linear geometry should be taken into account, during this step. The INC=5 indicates that the maximum number of increments during the step is 5.

This indicates that the step as a static load step. A time scale is assigned to the analysis. over which the load are applied gradually in increments. The size of the increments can be automatic, or defined by the user.

Here the second line means that the initial time increment is 1, the time period of the step is 1, there is no minimum value of the of the time increment allowed, and the maximum size increment is one.

This describes the boundary conditions of the modal, the first number is the start dof, the second number is the end dof, and the third number is the displacement of all the dofs from the start to end numbers.

Thus the node set pull(the ones around the edge) have a displacement of 0 in the x and y dof's. The second line defines a displacement of -0.05 in the z direction (downwards). The third line sets all the displacements of all the dof's to 0, for node 1000.

This defines the output of the nodal variables to a file, here we output the displacements (U) and the reaction forces at the fixed boundaries.

The values are outputted every 200 increments.

This define the output of element variables, the values are the averages of values extrapolated to the nodes of the elements in the set. The output will be split between elements with different properties.

Here we are output the stresses in the elements, every 200 increments.

This define the nodal variables to be send to the results file, here we output the displacement, reaction forces and contact forces, every increment.

This define the element variables to send the results file, here we output the stresses(S) and section strain, curvature change, and twist components(SE).

This outputs the contact variables, which are contact pressure, frictional shear stresses, contact opening(clearance) and relative tangential(slip) motions. Here we have specified that we want the variables for both the master and slave contact surfaces, we could have one or the other, or any other contact interactions if there were anymore. The values are sent to the file every increment.

This allows a degree of freedom to be monitored during the execution of the analysis. Here we choose to watch the 3rd DOF (Z direction) of node 1, the output to the message and status file is every increment.

This marks the end of the first step.

This marks the beginning of a new step, once again the non-linear geometry is taken into account.

This does the eigenvalue extraction, the stiffness determined at the previous step, is used as the basis for the extraction process. This allows small vibrations of a preloaded structure to be modelled.

The first ten eigenvalues, are to be calculated, the maximum frequency of interest is 1000Hz. The shift point is 45, eigenvalues closest to this point, will be extracted.

This resets the boundary conditions, the OP=NEW, removes all of the boundary condition's, and replaces them with new conditions.

The node set pull, is restrained, in the x,y, and z directions, and node 1000 is restrained fixed in all dofs.

This suppresses all output to the output file regarding element variables.

This is a request to output the eigenvalues to the output file.

This ends the step and the program.