Linear Triangular Element

Related terms:

Finite Element Method, Boundary Condition, Degrees of Freedom, Discretization,

Shape Function, Stiffness Matrix, Nodal Displacement, Nodal Point

View all Topics

Two-Dimensional Interval Finite Ele-

ment
Sukanta Nayak, Snehashish Chakraverty, in Interval Finite Element Method with
MATLAB, 2018

6.4 Linear Triangular Element

A linear triangular element is a two-dimensional finite element that has three nodes
and three sides shown in Fig. 6.8. It has three vertices and the nodes have coordinates
(x1, y1) , (x2, y2) , (x3, y3) in global Cartesian coordinate system. Each linear triangular
element has six degrees of freedom where each node contributes two degrees of
freedom, viz. u and v, the translation along global X and Y axes, respectively.

Fig. 6.8. Linear triangular element.

We assume that the linear interpolation function for displacement field is

(6.42)

The nodal displacements are

(6.43)
Eq. (6.43) can be represented in the following matrix form

(6.44)

Now solving for c0 , c1 , c2, we have

(6.45)

where .

Here,

Substituting the values of c0 , c1 , c2 in Eq. (6.42) and rearranging, we get

(6.46)

Eq. (6.46) can be represented in the following compact form

(6.47)

where the shape functions Ni , i = 1 , 2 , 3 are

(6.48)

In structural mechanic problems (Chakraverty, 2008), each node of a linear triangular

element possesses two degrees of freedom, viz. u and v. Hence, the displacements
can be written in terms of the nodal displacements at any internal point. The
displacements are

(6.49)

Finally, using finite element notations (Nayak and Chakraverty, 2015), the displace-
ment field in terms of the nodal degrees of freedom at an interior point can be
written as

(6.50)

where u1 , v1 , u2 , v2 , u3 , v3 are the nodal displacements at nodes 1, 2, and 3,

respectively.

> Read full chapter

The Principle of Minimum Poten-
tial Energy for Two-Dimensional and
Three-Dimensional Elements
Dimitrios G. Pavlou PhD, in Essentials of the Finite Element Method, 2015

10.3.1 The Linear Triangular Element (or CST Element)

We recall the linear triangular element shown in Figure 10.4. As it is known from the
mechanics of solids, the matrix { } is correlated to the field displacements {u(x, y)
(x, y)} by the following equation:

(10.193)

Combining the above equation with Equation (10.57), the following expression can
be derived

(10.194)

where i,  j,  m,  i,  j,  m can be obtained from Equations (10.43)–(10.48). Equation
(10.194) can be written in the following abbreviated form:

(10.195)

where

(10.196)

However, according to Hooke’s law for an isotropic plate in plane stress conditions,
the stress { } is correlated to {   } by Equation (9.136)

(10.197)

where the matrix [E ] is already known from Equation (9.137). Combining Equations
(10.195) and (10.197), it can be written:

(10.198)

Therefore, using Equations (10.195) and (10.198), the strain energy U given by
Equation (9.154) now takes the form

(10.199)

Since (a) the surface forces acting on an element are usually simulated by concentrat-
ed forces acting on the nodes, and (b) we can neglect the inertia forces (body forces)
because the problem is static, the work of the external forces is given by Equation
(9.157). The last equation can be written in the following abbreviated form:

(10.200)

where

(10.201)

and

(10.202)

Therefore, combining Equations (9.151), (10.199), and (10.200), the total potential
energy can now be written as:

(10.203)

Applying the principle of the MPE, the above equation yields

(10.204)

or

(10.205)

For an element of constant thickness t, the above equation can be written as

(10.206)

From the above equation, it can be concluded that the stiffness matrix {k} of the
linear triangular element (or CST element) can be derived by the following equation:

(10.207)

However, the matrices [E] and [B] are independent of x or y. Therefore, the above
equation can now be written in the following form:

(10.208)

where A is the area of the triangle given by Equation (10.49).

> Read full chapter

Numerical Finite and Boundary Ele-

ment Methods
Martin H. Sadd, in Elasticity (Third Edition), 2014

16.4 FEM problem application

Applications using the linear triangular element discretize the domain into a con-
nected set of such elements (see, for example, Figure 16.1). The mesh geometry
establishes which elements are interconnected and identifies those on the boundary
of the domain. Using computer implementation, each element in the mesh is
mapped or transformed onto a master element in a local coordinate system where all
calculations are done. The overall problem is then modeled by assembling the entire
set of elements through a process of invoking equilibrium at each node in the mesh.
This procedure creates a global assembled matrix system equation of similar form
as (16.3.13). Boundary conditions are then incorporated into this global system to
reduce the problem to a solvable set of algebraic equations for the unknown nodal
displacements. We do not pursue the theoretical and operational details in these
procedures, but rather focus attention on a particular example to illustrate some of
the key steps in the process.

EXAMPLE 16.1: ELASTIC PLATE UNDER UNIFORM TENSION

Consider the plane stress problem of an isotropic elastic plate under uniform tension
with zero body forces, as shown in Figure 16.3. For convenience, the plate is taken
with unit dimensions and thickness and is discretized into two triangular elements
as shown. This simple problem is chosen in order to demonstrate some of the
basic FEM solution procedures previously presented. More complex examples are
discussed in the next section to illustrate the general power and utility of the
numerical technique.

FIGURE 16.3. FEM Analysis of an Elastic Plate Under Uniform Tension.

The element mesh is labeled as shown with local node numbers within each element
and global node numbers (1–4) for the entire problem. We start by developing the
equation for each element and then assemble the two elements to model the entire
plate. For element 1, the geometric parameters are 1 = −1, 2 = 1, 3 = 0, 1 =
0, 2 = −1, 3 = 1, and A1 = 1/2. For the isotropic plane stress case, the element
equation follows from our previous work(16.4.1)In similar fashion for element 2,
1 = 0, 2 = 1, 3 = −1, 1 = −1, 2 = 0, 3 = 1, A1 = 1/2, and the element equation

becomes(16.4.2)These individual element equations are to be assembled to model

the plate, and this is carried out using the global node numbering format by
enforcing x and y equilibrium at each node. The final result is given by the assembled
global system(16.4.3)where Ui and Vi are the global x and y nodal displacements, and
and are the local stiffness components for elements 1 and 2 as given in relations
(16.4.1) and (16.4.2).

The next step is to use the problem boundary conditions to reduce this global
system. Because the plate is fixed along its left edge, U1 = V1 = U4 = V4 = 0.
Using the scheme presented in equations (16.3.18)–(16.3.20), the tractions on the
right edge are modeled by choosing = T/2, = 0, = T/2, = 0. These conditions
reduce the global system to(16.4.4)This result can then be solved for the nodal
unknowns, and for the case of material with properties E = 207 GPa and =
0.25, the solution is found to be(16.4.5)Note that the FEM displacements are not
symmetric as expected from analytical theory. This is caused by the fact that our
simple two-element discretization eliminated the symmetry in the original problem.
If another symmetric mesh were used, the displacements at nodes 2 and 3 would
then be symmetric. As a postprocessing step, the forces at nodes 1 and 4 could
now be computed by back-substituting solution (16.4.5) into the general equation
(16.4.3). Many of the basic steps in an FEM solution are demonstrated in this
hand-calculation example. However, the importance of the numerical method lies
in its computer implementation, and examples of this are now discussed.

> Read full chapter

Control volume finite element method

(CVFEM)
Mohsen Sheikholeslami Kandelousi, Davood Domairry Ganji, in Hydrothermal
Analysis in Engineering Using Control Volume Finite Element Method, 2015

1.3 Element and Interpolation Shape Functions

A building block of discretization is the triangular element (Fig. 1.2). For linear
triangular elements the node points are placed at the vertices. In Fig. 1.2, the
nodes, moving in a counterclockwise direction, are labeled 1, 2, and 3. Values of
the dependent variable are calculated and stored at these node points.
Figure 1.2. An element indicating the areas used in shape function definitions [11].

In this way, values at an arbitrary point (x, y) within the element can be approximated
with linear interpolation

(1.1)

where the constant coefficients a, b, and c satisfy the nodal relationships

(1.2)

Equation (1.1) can be more conveniently written in terms of the shape function N1,
N2, and N3, where

(1.3)

(1.4)

such that, over the element the continuous unknown field can be expressed as the
linear combination of the values at nodes i = 1,2,3:

(1.5)

With linear triangular elements a straightforward geometric derivation for the shape
functions can be obtained. With reference to Fig. 1.2, observe that the area of the
element is given by

(1.6)

and the area of the subelements with vertices at points (p, 2, 3), (p, 3, 1), and (p, 1, 2),
where p is an arbitrary and variable point in the element, are given by

(1.7)
Based on these definitions, it follows that the shape functions are given by

(1.8)

Note that, when point p coincides with node i (1, 2, or 3), the shape function Ni = 1,
and when point p is anywhere on the element side opposite node i, the associated
subelement area is zero, and, through Eq. (1.8), the shape function Ni = 0. Hence
the shape functions defined by Eq. (1.8) satisfy the required condition in Eq. (1.3).
Further, note that at any point p, the sum of the areas:

(1.9)

is such that the shape functions at (xp, yp) sum to unity. Hence the shape functions
defined by Eq. (1.8) also satisfy the conditions in Eq. (1.4). For future reference, it
is worthwhile to note that the following constants are the derivatives of the shape
functions in Eq. (1.8) over the element:

(1.10)

> Read full chapter

Detailed Explanation of Control Vol-

ume-based Finite Element Method
Mohsen Sheikholeslami, in Application of Control Volume Based Finite Element
Method (CVFEM) for Nanofluid Flow and Heat Transfer, 2019

1.3 The Element and the Interpolation Shape Functions

The building block of the discretization is the triangular element, Fig. 1.2. For
linear triangular elements the node points are placed at the vertices. In Fig. 1.2 the
nodes, moving in a counterclockwise direction, are labeled 1, 2, and 3. Values of the
dependent variable are calculated and stored at these node points.
Figure 1.2. An element indicating the areas used in shape function definitions.

In this way, values at an arbitrary point within the element can be approximated with
linear interpolation

(1.1)

where the constant coefficients a, b, and c satisfy the nodal relationships

(1.2)

Eq. (1.1) can be more conveniently written in terms of the shape function, , and ,
where

(1.3)

(1.4)

such that, over the element, the continuous unknown field can be expressed as the
linear combination of the values at nodes

(1.5)

With linear triangular elements a straightforward geometric derivation for the shape
functions can be obtained. With reference to Fig. 1.2, observe that the area of the
element is given by

(1.6)

and the area of the subelements with vertices at points and , where is an arbitrary
and variable point in the element, are given by

(1.7)

With these definitions it follows that the shape functions are given by
(1.8)

Note that, when point coincides with node , the shape function, and when point is
anywhere on the element side opposite node i, the associated subelement area is
zero, and through Eq. (1.8), the shape function. Hence the shape functions defined
by Eq. (1.8) satisfy the required condition in Eq. (1.3). Further, note that at any point
, the sum of the areas:

(1.9)

such that the shape functions at will sum to unity. Hence the shape functions defined
by Eq. (1.8) also satisfy the condition Eq. (1.4). For future reference it is worthwhile
to note that the derivatives of the shape functions in Eq. (1.8) over the element are
the following constants

(1.10)

> Read full chapter

Fluid Mechanics of Viscoelasticity

In Rheology Series, 1997

Linear and Quadratic Shape Functions

The number of nodes per element depends on the order of the interpolant and the
shape of the element. A linear triangular element has three nodes located at the
vertices of the element. Its nodal shape functions are given by

(C7.1)

The triangle defined by this interpolation has straight sides, due to the degree of
interpolation. To model a curved triangular element accurately, a quadratic interpo-
lation is employed, with the following nodal shape functions:

(C7.2)

Note that there are six nodes per element, the three extra nodes being located at the
midsides of the triangle.
One can also have a linear quadrilateral, i.e., with straight sides, having four nodes
located at the vertices with associated nodal shape functions:

(C7.3)

Quadrilaterals of curved sides can be modeled by a quadratic quadrilateral witheight

nodes and nodal shape functions:

(C7.4)

with − 1 _ “1, “2 − 1. It should be noted that

(C7.5)

which is a property of the Lagrangian interpolation function.

> Read full chapter

Volume 1
Kiyoshi Yonekawa, Mutsuto Kawahara, in Computational Mechanics–New Frontiers
for the New Millennium, 2001

FINITE ELEMENT METHOD

Basic equation
Linear shallow water equation, which consist of momentum equation and continuity
equation are shown as:

(1)

(2)

where ui is velocity, is water elevation, h is water depth and g is the gravity

acceleration.

Boundary conditions
In the boundary , boundary conditions are given as:
(3)

(4)

where ˆ means given value and ni means direction cosine.

Finite element equations

To derive finite element equations from shallow water equation, it can be discretized
spatially and temporally. The Galerkin method with linear triangular element is
applied to the spatial discretization. The Balancing Tensor Diffusivity is applied to
the temporal discretization. The finite element equations are written as:

(5)

(6)

(7)

> Read full chapter

Computer Implementation of the CBS

Algorithm
O.C. Zienkiewicz, ... P. Nithiarasu, in The Finite Element Method for Fluid Dynamics
(Seventh Edition), 2014

15.1 Introduction
In this chapter we consider some essential steps in the computer implementation
of the CBS algorithm on structured or unstructured finite element grids. Only
linear triangular elements are used and the description given here is intended for
a two-dimensional version of the program. The program source and user manual
along with several solved problems are available at the website

http://www.zetacomp.com

free of charge. Readers also will find documentation of the code and other details
on this website including new results and corrections.

The program can be used to solve the following different categories of fluid mechan-
ics problems:
1. Compressible viscous and inviscid flow problems

2. Incompressible viscous flows

3. Incompressible flows with heat transfer

With further simple modifications, many other problems such as turbulent flows,
solidification, mass transfer, free surfaces, etc., can be solved. The procedures pre-
sented are largely based on the theory presented in Chapter 3. The program is written
in Fortran and it is assumed that the reader is familiar with coding and finite element
procedures discussed in this book [1,2].

We call the present program CBSflow since it is based on the CBS algorithm
discussed in Chapter 3 of this volume. We prefer to keep the compressible and
incompressible flow codes separate to avoid any confusion. However an experienced
programmer can incorporate both parts into a single code without much difficulty.
Each program listing is accompanied by some model problems which may be used
to validate an installation. In addition to the model inputs to programs, a complete
user manual is available to describe each part of the program in detail. Any error
reported by users will be corrected and the program will be continuously updated
by the authors.

The modules are (1) the data input module with preprocessing; (2) the solution
module; and (3) the output module. The CBSflow program contains the source
statements to solve transient Navier-Stokes (or Euler-Stokes) equations iteratively.
Here there are many possibilities such as fully explicit forms, semi-implicit forms,
quasi-implicit forms, and fully implicit forms as presented in Chapter 3. We concen-
trate mainly on the first two forms which require small memory and simple solution
procedures compared to other forms.

In both the compressible and incompressible flow codes, only nondimensional

equations are used. The reader is referred to the appropriate chapters (Chapters 3,
4, and 5Chapter 3Chapter 4Chapter 5) for different nondimensional parameters.

In Section 15.2 we describe the essential features of data input to the program. Here
either structured or unstructured meshes can be used to divide the problem domain
into finite elements. Section 15.3 explains how the steps of the CBS algorithm are
implemented. In that section, we briefly remark on the options available for shock
capturing, various methods of time stepping, and different procedures for equation
solving. In Section 15.4, the output generated by the program and postprocessing
procedures are considered.

> Read full chapter

Unstructured grid generation
M. Farrashkhalvat, J.P. Miles, in Basic Structured Grid Generation, 2003

8.4 Solving hosted equations using finite elements

Here we present a very brief introduction to the solution of field equations (the
hosted equations) in a domain which has been triangulated, using linear triangular
elements.

For further information, standard texts on finite element methods such as Hinton
and Owen (1979), Taylor and Hughes (1981), and George (1991) may be consulted.

Suppose that the field equation for a field quantity , subject to certain boundary
conditions, in a planar domain D (assumed simply connected here) is

(8.19)

where L is a second order partial differential operator. Suppose also that the domain
has been triangulated, with n triangles and m points (nodes).

In the Method of Weighted Residuals, we seek an approximate solution for from

the family of functions

(8.20)

where usually the function f0 is chosen to satisfy the given boundary conditions
and the fi s satisfy homogeneous boundary conditions in order to make satisfy the
boundary conditions for any choice of the constant coefficients i. In addition, we
choose a set of N weighting functions Wi(x, y) and impose the conditions

(8.21)

to make as close to zero as possible in some sense over the domain D. These
equations will generate a set of N equations for the N unknown constants i and
hence the best (in some sense) approximate solution from our original set.

In the special case of the Galerkin Method, the weighting functions are taken to be
identical to the approximating functions, which gives the set of equations

(8.22)

with given by eqn (8.20).

In the particular form of the finite-element method which is appropriate here for a
plane triangulated region, we look for approximate solutions in the form
(8.23)

where the i s are constant coefficients to be determined and the known Ri(x,
y) functions are called roof functions, which are continuous and piecewise-linear in
both x and y, such that Ri takes the value unity at the ith node and zero at all other
nodes. Each node of the triangulation is assigned a global number between 1 and
m. Moreover, a node constitutes one of three vertices of a triangular element, and
is assigned a local number between 1 and 3, these numbers being ordered in an
anti-clockwise sense. This assigning of numbers, together with the numbering of
triangular elements and the listing of cartesian co-ordinates of all nodes, is carried
out in the Delaunay grid generation program delaimay1.f listed below.

Roof functions may be built up from linear shape functions defined over each
triangular element. Here the suffix refers to the local number of a vertex;

takes the value unity at the ith vertex and zero at the other two. Over a particular
triangle with vertices numbered 1, 2, 3 in anti-clockwise manner, having known
cartesian co-ordinates (x1, y1), (x2, y2), (x3, y3), respectively, it is easy to show that a
shape function linear in x and y and taking the value 1 at the vertex 1 and zero at the
vertices 2 and 3 is given by

(8.24)

where a1 = x2 y3 – x3 y2, b1 = y2 – y3, C1 = x3 – x2, and Ae is the area of the triangle.

Expressions for similar shape functions and which take the value 1 at the vertices
2, 3, respectively, and zero at the other vertices may be immediately written down
by cyclic permutation of the suffixes in eqn (8.24). Over a typical triangular element
of the triangulation, we can then write the approximating solution as

(8.25)

where now the coefficients represent the values of e at the vertices.

As an example of a field equation, we take Poisson's equation

(8.26)

where is a constant and in heat transfer problems Q(x, y) could be a heat source
function. The boundary conditions are taken to be that is specified on a part C
of the boundary and the normal derivative ∂ /∂n (in the direction of the outward
normal) takes the value q (possibly a function of position) on the remaining part Cq.
This means that at nodes on that part of the boundary of the triangulation which
approximates C the value of the coefficients i in eqn (8.23) is effectively specified.

Suppose that N of the nodes are in the interior of the triangulation or on the
boundary Cq and that nodes with global vertex numbers N + 1, N + 2,…, m lie on
that part of the boundary which approximates C . Then instead of eqn (8.23) we can
put

(8.27)

where in the first term the coefficients i are known from the boundary conditions.
This representation of the family of approximating functions has the same form as
eqn (8.20), since the roof functions Ri(x, y) take the value zero at the boundary nodes
N + 1, N + 2,…, m. The Galerkin approach gives the integral form

(8.28)

but to be able to deal with approximating functions with discontinuous slopes across
triangular edges, we need to transform the double integral by integration by parts
(the

Divergence Theorem, or Green's formula) to

(8.29)

where the contour integral is taken anti-clockwise around the boundary C of D. If

we replace by q on Cq, and Ri by zero on C , we obtain

(8.30)

Using eqn (8.27), this becomes, for i = i,…, N,

(8.31)

The contribution of a triangular element to the terms

in eqn (8.31) is then, using eqn (8.25),

(8.32)

where is clearly symmetric and can be regarded as the element stiffness matrix. By
eqn (8.24) we can write, in the present case,

(8.33)

For a given triangulation, reducing eqns (8.31) to a global matrix equation

(8.34)

involves first the evaluation of the stiffness matrices for all elements and then
assembling and adding them into a large N × N matrix Kij. This is usually a straight
forward task once a connectivity table has been established listing the global and local
vertex numbers for each element in the triangulation. This enables global numbers
to be made to correspond to the row (and column) numbers of Kij so that the 3 ×
3 entries can be added in the correct places. The ‘force’ terms Fi can be calculated
from the RHS of eqn (8.31) using the values of Q and q at appropriate nodes. Finally
a

numerical scheme must be used to solve the global matrix equation for the unknown
coefficients i.

Similar procedures can be used to solve non-linear hosted equations, such as

the Navier-Stokes equations, using iterative methods based on an initial approxi-
mate solution. For example, see Rao, S.S. (1982) and Lewis, Morgan, Thomas, and
Seetharamu (1996).

> Read full chapter

Coupled Finite Element–Agent-Based

Models for the Simulation of Vascular
Growth and Remodeling
D. Nolan, C. Lally, in Numerical Methods and Advanced Simulation in Biomechanics
and Biological Processes, 2018

16.3.2.2.4 Element Selection

The choice of element type is important not only for modeling incompressibility,
but it will determine the accuracy of the solution. Fig. 16.3A shows a schematic for
a linear triangular element with its integration point. This element is also known
as a constant strain triangle as the single integration point means that only one
stress/strain state exists within its boundary. These elements are useful because
they can mesh complex geometries; however, because of the single integration
point, relatively more elements must be used when modeling a heterogeneous
stress field compared with quadrilateral elements. Fig. 16.3B shows a schematic of a
quadrilateral element and its integration points for both full integration and reduced
integration formulations. In the full integration formulation the four internal inte-
gration points permit a complex stress gradient within an element, which generally
means that fewer elements are required. In problems where there is a high degree of
bending the full integration elements are susceptible to a numerical problem known
as shear locking, where a bending deformation is misinterpreted as a shear deforma-
tion at the integration points. To prevent this occurring, a reduced integration linear
quadrilateral element may be used instead. However, reduced integration elements
do not permit internal stress gradients because of their single integration point and
suffer from their own numerical problem known as hourglassing.

Figure 16.3. (A) Linear triangular element showing the nodes and integration point.
(B) Linear quadrilateral element showing the integration points for a full integration
and reduced integration formulation.

Interestingly, reduced integration elements are often the default setting in com-
mercial FE codes, such as Abaqus and Ansys, because they have a relatively low
computational cost, the hourglassing problem can be numerically managed, and
they behave well when an element is highly distorted. The novice user cannot make
serious unconscious errors by using this default. In the authors' opinion, however,
full integration elements should be used where appropriate. In the balloon angio-
plasty and stent implantation simulations, there are no bending loading modes.
However, in the stent implantation simulation, there is a high gradient of stress at
the site of implantation, which warrants the use of full integration elements.

Fig. 16.4 demonstrates the importance of conducting a mesh convergence study,

particularly in regions with a heterogeneous stress field such as that found around
the stent strut. Here, a model using a low-density mesh fails to capture the stress
field behind the stent strut.

Figure 16.4. Contour plot of von Mises stress in the region around the implanted
stent strut for (A) a low mesh density (element length 7 μm and (B) a high mesh
density (element length 3 μm). Although the stresses are similar in locations
furthest from the stent, the stress immediately behind the strut is very different
for (A) and (B) demonstrating the importance of having a sufficiently dense
 mesh in regions of highly heterogeneous stress. Therefore, approximately 60 plane
strain elements along the contact edge of the stent are required in the undeformed
configuration.

> Read full chapter

ScienceDirect is Elsevier’s leading information solution for researchers.

Copyright © 2018 Elsevier B.V. or its licensors or contributors. ScienceDirect ® is a registered trademark of Elsevier B.V. Terms and conditions apply.