12

where b is the width of the beam and h the height of the beam and E the Young™s

modulus of the material. The ratio h/L of height over length equal to 0.1 is chosen.

Using the meshes depicted in Fig. 18.3, the results of Table 18.1 are obtained,

which presents the ratio of uL and the computed displacement uc at x = L.

It is clear that the displacement of the beam, using a mesh of linear trian-

gles is much too small. The poor performance of the linear triangle can easily

be understood; because of the linear interpolation of the displacement ¬eld u, the

associated strains computed from

Infinitesimal strain elasticity problems

324

Figure 18.3

Triangular and quadrilateral element distribution for the bar problem.

µ = Bu (18.69)

∼

∼

are constant per element. The bi-linear quadrilateral element, on the other hand,

is clearly enhanced. A typical shape function, of for example the ¬rst node in an

element, is given by (with respect to the local coordinate system)

1

N1 = (1 ’ ξ ) (1 ’ ·) . (18.70)

4

Hence

1

N1 = (1 ’ ξ ’ · + ξ ·) , (18.71)

4

which means that an additional non-linear term is present in the shape functions.

Therefore a linear variation of the stress ¬eld within an element is represented.

Two remarks have to be made at this point:

• The numerical analysis in this example was based on plane stress theory, while in the

chapter the equations were elaborated for a plane strain problem. How this elaboration

is done for a plane stress problem is discussed in Exercise 18.1.

• Strictly speaking, a concentrated force in one node of the mesh is not correct. Mesh

re¬nement in this case would lead to an in¬nite displacement of the node where the

force is acting. This can be avoided by applying a distributed load over a small part of

the beam.

Exercises

18.1 For Hooke™s law the Cauchy stress tensor σ is related to the in¬nitesimal

strain tensor µ via

σ = Ktr( µ) I + 2Gµ d .

(a) What is tr(µ) for the plane strain case?

(b) What is tr(µ) for the plane stress case?

325 Exercises

What are the non-zero components of µ and σ for the plane stress

(c)

case?

What are the non-zero components of µ and σ for the plane strain

(d)

case?

(e) Consider the plane stress case. Let

⎡ ¤ ⎡ ¤

σxx µxx

⎢ ⎥ ⎢ ⎥

σ = ⎣ σyy ¦ ∼ = ⎣ µyy ¦.

∼

σxy µxy

What is the related H matrix for this case, using

σ = H µ?

∼ ∼

(f) For the plane strain case the matrix H is given by

⎡ ¤ ⎡ ¤

4 ’2

110 0

⎢ ⎥ G⎢ ⎥

H = K ⎣ 1 1 0 ¦ + ⎣ ’2 4 0 ¦.

3

000 0 0 3

Rewrite this matrix in terms of Young™s modulus E and Poisson™s ratio

ν. Is the resulting matrix linearly dependent on E?

18.2 Consider the mesh given in the ¬gure below. The mesh consists of two

linear triangular elements and a linear elasticity formulation applies to this

element con¬guration. The solution u is given by

∼

uT = [ u1 w1 u2 w2 u3 w3 u4 w4 ] ,

∼

where u and v denote the displacements in the x- and y-direction, respec-

tively.

(a) What is a possible top array of this element con¬guration assuming

equal material and type identi¬ers for both elements?

(b) What is the pos array for this element con¬guration?

4 3

(2)

y

(1)

x

1 2

(c) What is the dest array?

(d) Based on the pos array the non-zero entries of the stiffness matrix

can be identi¬ed. Visualize the non-zero entries of the stiffness

matrix.

Infinitesimal strain elasticity problems

326

(e)

Suppose that the boundary nodes in the array usercurves are

stored as

usercurves=[

13

34

24

1 2];

The solution u is stored in the array sol. How are the displace-

∼

ments in the x-direction extracted from the sol array along the ¬rst

usercurve?

(f) Let the solution array sol and the global stiffness matrix q be given.

Suppose that both displacements at the nodes 1 and 2 are suppressed.

Compute the reaction forces in these nodes. Compute the total reac-

tion force in the y-direction along the boundary containing the nodes

1 and 2.

18.3 Consider the bi-linear element of the ¬gure below. It covers exactly the

domain ’1 ¤ x ¤ 1 and ’1 ¤ y ¤ 1. Assume that the plane strain

condition applies.

(a) Compute the strain displacement matrix B for this element in point

( x, y) = ( 1 , 1 ).

22

(b) Let the nodal displacements for this element be given by

uT = .

0 0 0 0 1 0 1 0

∼e

What are the strains in ( x, y) = ( 1 , 1 )?

22

y

(“1, 1) (1, 1)

4 3

x

1 2

(“1, “1) (1, “1)

If G = 1 and K = 2 what are the stress components σxx , σyy and σxy

(c)

at ( x, y) = ( 0, 0)? What is the value of σzz ?

18.4 Given the solution vector sol and the pos array,

(a) Extract the solution vector for a given element ielem.

(b) Compute the strain at an integration point.

(c) Compute the stress at an integration point.

327 Exercises

18.5 The m-¬le demo_bend in the directory twode analyses the pure bending

of a single element (for quadrilaterals) or two elements (for triangles). The

analysis is based on plane stress theory. The geometry is a simple square

domain of dimensions 1 — 1. Along the left edge the displacements in the

x-direction are set to zero, while at the lower left corner the displacement in

the y-direction is set to zero to prevent rigid body motions. The two nodes

at the right edge are loaded with a force F of opposite sign, to represent

a pure bending moment. Investigate the stress ¬eld for various elements:

linear triangle, bi-linear quadrilateral, and their quadratic equivalents.

Explain the observed differences. To make these choices use itype and

norder:

itype = 1 : quadrilateral element

itype = 20 : triangular element

norder = 1 :(bi-)linear element

norder = 2 :(bi-)quadratic element

18.6 In a shearing test a rectangular piece of material is clamped between a

top and bottom plate, as schematically represented in the ¬gure. This

experiment is generally set-up to represent the so-called ˜simple-shear™

con¬guration. In the simple-shear con¬guration the strain tensor is given as

µ = µxy ex ey + µxy ey ex

using the symmetry of the strain tensor. As a consequence, the stress-strain

relation according to Hooke™s law reduces to

σxy = 2Gµxy .

Hence, measuring the clamp forces and the shear displacement provides

h

a direct means to identify the shear modulus G. However, the ˜simple-

shear™ state is dif¬cult to realize experimentally, since the con¬guration

of the ¬gure does not exactly represent the simple-shear case. This may be

analysed using the m-¬le demo_shear.

(a) Analyse the shear and the simple-shear case using this m-¬le. What

is the difference in boundary conditions for these two cases?

Infinitesimal strain elasticity problems

328

(b) Why is the simple_shear=0 case not equal to the exact simple-

shear case?

What ratio /h is required to measure the modulus G within 10 %

(c)

accuracy? How many elements did you use to obtain this result?

Explain the way the m-¬le computes the modulus G.

References

[1] Adams, R. A. (2003) Calculus: A Complete Course (Addison, Wesley, Longman).

[2] Brekelmans, W. A. M., Poort, H. W. and Slooff, T. J. J. H. (1972) A new method to

analyse the mechanical behaviour of skeletal parts. Acta Orthop. Scandnav. 43,

301“17.

[3] Brooks, A. N. and Hughes, T. J. R. (1990) Streamline upwind/Petrov-Galerkin

formulations for convection dominated ¬‚ows with particular emphasis on the

incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics

and Engineering Archive “ Special edition, 199“259.

[4] Carslaw, H. S. and Jaeger, J. C. (1980) Conduction of Heat in Solids (Clarendon

Press).

[5] Cacciola, G. R. C. (1998) Design, simulation and manufacturing of ¬bre reinforced

polymer heart valves. Ph.D. thesis, Eindhoven University of Technology.

[6] Frijns, A. J. H., Huyghe, J. M. R. J. and Janssen, J. D. (1997) A validation of the

quadriphasic mixture theory for intervertebral disc tissue. Int. J. Engng. Sci.

35(15), 1419“29.

[7] Fung, Y. C. (1990) Biomechanics: Motion, Flow, Stress, and Growth

(Springer-Verlag).

[8] Fung, Y. C. (1993) Biomechanics: Mechanical Properties of Living Tissues, 2nd

edition (Springer-Verlag).

[9] Hill, A. V. (1938) The heat of shortening and the dynamic constants in muscle.

Proc. Roy. Soc. London 126, 136“65.

[10] Hughes, T. J. R. (1987) The Finite Element Method (Prentice Hall).

[11] Huyghe, J. M. R. J., Arts, T. and Campen, D. H. van (1992) Porous medium

¬nite element model of the beating left ventricle. Am.J.Physiol. 262

H1256“H1267.

[12] Huxley, A. F. (1957). Muscle structure and theory of contraction. Prog. Biochem.

Biophysic. Chem., 255“318.

[13] Mow, V. C., Kuei, S. C. and Lai, W. M. (1980) Biphasic creep and stress

relaxation of articular cartilage in compression. J. Biomech. Engng. 102 73“84.

[14] Oomens, C. W. J. (1985) A mixture approach to the mechanics of skin and

subcutis. Ph.D. thesis, Twente University of Technology.

References

330

[15] Oomens, C. W. J., Maenhout, M., Oijen, C. H. van, Drost, M. R. and

Baaijens, F. P. T. (2003) Finite element modelling of contracting skeletal muscle.

Phil. Trans. R. Soc. London B 358, 1453“60.

[16] Sengers, B. G., Oomens, C. W. J., Donkelaar, C. C. van, and Baaijens, F. P. T.

(2005) A computational study of culture conditions and nutrient supply in cartilage

tissue engineering. Biotechnology Progress 21, 1252“1261.

[17] Sengers, B. G. (2005) Modeling the development of tissue engineered cartilage.