A simple contour plotting program for finite element output

A simple contour plotting program for finite element output J. F. L Y N E S S * and L. A S Q U I T H "

Ulster Polytechnic, Shore Road, Newtownabbey, County Antrim BT37 OQB

The algorithm described is similar to that developed by D. Y. Kan in University College, Swansea, in the late 1960s. Piecewise linear contour plots are produced for twodimensional regions. The method differs from that of Kan in the drawing of a complete contour continuously rather than drawing contours element by element. The computer equipment used was a V A X l l / 7 8 0 , Tektronix4051 graphics display terminal and a Tektronix 4662 X-Y plotter.

A n algorithm to produce piecewise linear contour plots for finite e l e m e n t o u t p u t is described. Several

examples are shown to demonstrate the application of the algorithm.

INTRODUCTION The finite element method has proved itself a versatile technique for the solution of problems in structural mechanics, stress analysis, fluid mechanics and other field problems. The essence of the method is the sub-division of a domain, or continuum into a finite number of parts, or elements, the behavior of which is specified by a finite number of parameters. A matrix equation results which is solved to produce nodal values for the unknowns. The matrix equations are solved by computer and the output is usually a list of nodal values of the unknowns, functions or derivatives of the unknowns. Lists of output data are difficult to appraise quickly and the following simple contour plot algorithm can be added as a post-processor to give a visual representation of results.

FINITE ELEMENT DISCRETIZATION The various routes to finite dement discretization have been given by Zienkiewicz I and are shown in Fig. 1. An unknown function ~ is generally approximated within an element as a function of position and element nodal values. = approximation to unknown function ~(x, y) n

i=1 ~i i Nl(x, y) n

* Reader, School of Civil Engineering, Ulster Polytechnic. ~"Programmer,Computer Services,Ulster Polytechnic.

= = = =

value of unknown at node i node number shape function corresponding to node i total number of nodes around element

PHYSICAL PROBLEM

I

I

Variational principles

I

I

Weighted integrals of governing partial differential equations

I

I

Physical principles (e.g. minimization of strain energy)

Physical statements (e.g. virtual work)

I

Miscellaneous weighted residuals

Constrained Lagrangian forms

Collocation methods (point or subdomain)

Penalty function forms

Galerkin weighted residuals Least squares forms

t

Numerical approximation [ K ] { 0 } + { F } = 0

F~gure 1.

Routes to finite element approximation

041-1195/83/010023-09 $2.00 © 1983 CML Publications

Adv. Eng. Software, 1983, Vol. 5, No. 1

23

(a) Three node triangular element

Y 2

3 i=1

1

1

Ni(x, y) = ~

(ai + bix +ciY)

a i = x/y m -- XmY] bi = Y / - - Ym Ci : X m --X]

1 xi

v

2A=

×

1 xj

1

Figure 2.

Yt y!

I

Xm Ym

(b) Four node isoparametric quadrilateral element

3

x, y = global co-ordinates

Y

~, 7/= local (dement) co-ordinates 4 i=1

~o = ~ i

rio = rt,r/i

N] (~, n)--- ~($o + 1)(7/o + 1) 4

x= E lynx, i=1 4

Y = E NIYi '"'

l=l

v

×

F~ure 3. system

xi, Yi = co-ordinates o f element nodes in global

(c) Eight node isoparametric quadrilateral element:

Y

8

= Z NI ~t i--1

Go= ~ i

rio = ririi

Corner nodes: NI(~, r/) = ~(1 + ~o)(1 + rio)(~o + rio -- 1) Midside nodes:

%=--1 2

8

~l = 0

ii(~,

77i = 0

N/(~, rt) = ½(1 +/jo)(1 - ri2)

8

x=

E N]xt i=l

X

Figure 4.

24

Global to local mapping functions

Adv. Eng. Software, 1983, Vol. 5, No. 1

8

Ni Yi i=1

ri) = ½(1 -- ~2)(1 + rio)

-----BOUNDARY J •

~

•

/

REGIONDIVIDEDINTO TRIANGULARELEMENTS

18

'/

+C. CONTOURVALUE((J/~i~¢E'~ J )'

EGIONDIVIDEDINTO 10 QUADRILATERALELEMENTS

Figure Z Determination of l c in terms of contour and nodal values

REGIONDIVIDEDINTO 8 CURVILINEARQUADRILATERAL ELEMENTS F~ure 5.

Division of region into elements

v~ Figure 8. Determination of Ic in terms of crossing point and nodal co-ordinates

CEONGT~EO

Figure 9.

CONTOURS

~

CLOSED

Unconnected contour segments

CONTOUR

/~~I

SEGMENTS~~~ 41

8

Figure 6.

Partition of quadrilateral elements

Figure 10. EDGE

7

9

10

11

Use of arrays MFLAG, MESHMAP, MESH-

Adv. Eng. Software, 1983, 17ol. 5, No. 1

25

MAIN P R O G R A M

START

/ l

)

TOPOLOGY FILENAME

/

ASSEMBLE TRIANGULAR MESH

NO --

YES

+

ASSEMBLE MESH OUTLINE IN SHAPE

ASSEMBLE FULL MESH I N. SHAPE

+ /G~T CO.TOUR/ INTERVAL /

T DRAW SHAPE / CALCULATE CONTOUR VA LUE

YES

J

ASSEMBLE CONTOUR IN SHAPE

I-

41

YES

YES

NO

/ F~gure 11.

26

Mainprogram

Adv. Eng. Software, 1983, Vol. 5, No. 1

SURFACE

(

~NO )

SUBROUTINE

CONTUR

t

INITIALISE FLAG ARRAY TO ZERO

q

IF ELEMENT IS CROSSED BY CONTOUR AND CONTOUR CUTS EDGE OFMESH,FLAG IS SET TO-1. ELSE FLAG IS SET TO 1

1

FLAG ALL ELEMENTS CROSSED BY CONTOUR INITIALISE ELEMENT INDEX TO ZERO

®

4

SEARCH FLAG ARRAY FOR FIRST-VE VALUE ELSE FIRST +VE VALUE,SET ELEMENT INDEX TO FLAG INDEX

I

~ RETURN 2

YES

[E. DEAL WITH ALL CONTOUR SEGMENTS WHICH ARE BROKEN BEFORE HANDLING CLOSED SEGMENTS A LOGICAL FLAG "CLOSED' IS USED TO DEFINE TOPOLOGY OF CURRENT SEGMENT

~ , ~NO

I

DETERMINE THE NODE NUMBERS OF THE SIDES OF THE ELEMENT CROSSED BY THE CONTOUR

~

ENT

NO

ENSURE START J AT EDGE OF MESHJ

ET FLAG ARRAY ELEMENTJ OF CURRENTELEMENT I

IS

INDEX TO ZERO AND I CALCULATE CO ORDINATES J OF CONTOUR INTERSECTION I

YES / I S - ~

SIDLE DEFINE D EDGE SIDE

YES

i

i

C~ENT

INITIALISE NEXT ELEMENT I z o. lNO , ELEMENT WITH SIDE COMMON SET NEXT ELEMENT FLAG TO THIS VALUE OF FLAG ARRAY I

I"

OTHER SIDE OF NEXT NO ELEMENT CUT BY CONTOUR

ISL~ FL~RO

YES jCALCULATE LAST POINT OF I ITHIS SEGMENT OF CONTOUR

F~ure 12.

Sub-routine CONTUR

Adv. Eng. Software, 1983, Vol. 5, No. 1

27

The matrix equation [K]{¢} + {F} = 0 is solved for {¢}. Output may consist of ¢ (e.g. displacement, temperature, velocity potential), its derivatives (e.g. strain, velocity) or some function of ¢ (e.g. stress, flow rate). A variety of elements are available for finite element modelling, the most commonly used two dimensional elements are described below. It can be seen that elements may have straight or curved sides. Examples of continua divided into finite elements are given in Fig. 5. In regions with steep gradients of the unknown function a timer mesh division is used (e.g. around stress concentrations). CONTOUR PLOTTING ALGORITHM

8

Input element topology and nodal co-ordinates Input vector of unknowns {~} Select contour interval Loop through all contours Loop through all elements containing contour Loop through element sides Determine where contour cuts sides Continue loop through sides Draw line joining points on two sides where contour crosses Continue loop through dements Quadrilateral elements are partitioned to form two triangles, as shown in Fig. 6. Curvilinear quadrilaterals are approximated with straight sides, this leads to some loss of information on nodal values. EXPRESSION OF CONTOUR CO-ORDINATES IN TERMS OF NODAL CO-ORDINATES AND NODAL VALUES

F~gure 13. Torsion o f square section

TOI~EON T ~ A R FJ.~S FOR CONTOUR PLOT PLOT ~ 1 CONTOUR "rNTERVAL O, E:I~eE÷81

Figure 14.

28

Triangular and quadrilateral finite element meshes

Adv. Eng. Software, 1983, Vol. 5, No. 1

Consider an element side i~ with nodal values ¢i, ¢1. Let l = length of side i], ¢c = contour value under considera-

TOP.$ZON OlJAD ELEMENTS FOR CONTOUR PLOT PLOT ~ R

I CONTOUR :]L'rlI"IERVAL e . r l 4 l l ~ 4 1 1

o

TORS'rON TR.1:ANGULAR ELEHENTS FOR CONTOUR PLOT PLOT NUHBF_.R I CONTOUR ]'NTERVAL 8.1ESeE÷81

F~gure 15.

TORSZON OUAD ELEMENTS FOR CONTOUR PLOT PLOT ~ R I CONTOUR ZNTERVAL 0 . S M + O 1

Stress function contours for triangular and quadrilateral meshes

=© =(23

X

o

Figure 16.

Torsion of cylinder

tion, lc= distance along side from node ] to the point where the ~c contour crosses side. From Fig. 7 by similar triangles:

By substituting for 1c from (1):

yc=yi+~)(yl--yi)

(3)

Similarly for the x co-ordinate:

.'. lc = l \ ~ !

(1)

--

(4)

The co-ordinates of the point (x~ Yc) can, therefore, be calculated from equations (3) and (4).

From Fig. 8 by similar triangles: Yj--Yf Yc--Yi t lc

le =(Ye--Yi] \YJ --YN

xc= x t + [ gi--~fl~xl

PROGRAMMING THE ALGORITHM

(2)

A given contour may consist of several unconnected segments which can be broken (i.e. cutting edge of region)

Adv. Eng. Software, 1983, Vol. 5, No. 1

29

T ~ ~OT

CY~R

C ~ V ~ ! ~

~

Figure 17.

/

~ Z~AL

CONTOUR PLOT O.I÷ll

~

PLOT 0 . 2 ~

/ HEAT FLOW RANDCbt PLOT NUPIBER

/

f

/

/

J

/

/

C ~ V ~ ~ I CONTOUR Z ~

Meshand stress function plots for torsion of cylinder

i

/

T ~ C Y ~ PLOT ~

_i

/

Figure 20. wall

ELEMENTS CONTOJR PLOT I CONTOUR I N T E R ' v ' & L 9 . 1 8 8 8 E + 0 0

Mesh and isotherm plot for heat flow through

i

x,~ t ~ ~ . ~ . ~ Figure 18.

CCNTOUR PLOT

II]l

~tLA7 FLO IA RArJDOM # L O T "~UP!BER

f"

ELE ,E ; .;

Finite element mesh

~mmmmmmii n|mmmmmml

nmmmmmnl

I~mmlmlii m~ummmumn

9~IIWfiIGi rlmnllW~miE~l InlilIPHKIII l ~ n m m n m n n

i~vnEnp~ CARI.ZI~F'ORD LOUGH AREA MAP REFS d l O l O , dlilP.I, d31J211,d31110 PLOT NUMBER ! CONTOUR ZNTERVAL O. S(IIlIE+03

CARLINGF'ORD LOUGH ARENA I'IAp RF..F$ dlOtO.dltl~O.d~l~k~l, d3010 PLOT ~ R I C..ONTOUR ZNTERVAL 8.16QQ~+EI3

Figure 19.

30

Stream function contours

Adv. Eng. Software, 1983, Vol. 5, No. 1

Figure 21.

Contour plots around Carlingford Lough

or dosed as shown in Fig. 9. After data input and partitioning of quadrilateral elements the program uses five arrays described below: NP = number of nodes in triangulated mesh NE = number of triangular elements COORD (NP, 2) LNODS (NE, 3)

nodal co-ordinates triangular element topology, list of nodes around element, consistently clockwise/anticlockwise for mesh MESHEDGE (NE, 3) ordered list of adjacent elements for every element in the mesh MFLAG (NE) flag array showing elements crossed by a particular contour = 0 for element not crossed by contour = 1 for element crossed by contour = -- 1 for element crossed by contour which also crosses a mesh edge side. The mesh edge sides are found by selecting pairs of node numbers on each element and comparing the pair with sides of all other elements. If the pair is not common to any other element pair then they form part of the mesh boundary and the node numbers of the pair are put into array MESHEDGE which is then ordered into clockwise/ anticlockwise pairs around the mesh boundary. For a particular contour, elements are first searched to find those which are crossed by the contour. The array MFLAG has been initialized to zero, all elements crossed by the contour are flagged with 1. Using array MESHEDGE those elements on the edge of the mesh crossed by the contour are found and flagged in array MFLAG with --1. An example of the use of the arrays is shown below for Fig. 10. Element No. 1

2 3 4 5 6 7 8 9 10 11

MESHMAP II III II III IV I VI III viii v x

VI 0 IV V X VI VIII VII x IX 0

0

0 VIII 0 0 0 0 IX 0 xI 0

MFLAG

MESHEDGE

1

1

3

- 1 1 1 1 -1 0 1 1 1 0

2 3 7 11 10 9 8 4

3 7 11 10 9 8 4 1

When array MFLAG has been found for a contour, calculation of points on piecewise linear segments may start. MFLAG is searched for the f i r s t - 1, contour segments starting from the boundary of the mesh are processed first, the co-ordinates of the start of the contour segment are calculated using equations (3) and (4). Adjacent elements are found from MESHMAP, contour co-ordinates calculated and flags set equal to 0 in MFLAG, array SHAPE contains co-ordinates for the current contour

segment which is drawn piecewise linearly. After all contour segments crossing the mesh edge have been drawn closed contour segments are drawn. Array MFLAG is zeroed after all contour segments for a particular contour value have been drawn. Contour intervals are checked to ensure drawing clarity and for contours passing through nodes. If a contour passes through a node the nodal value is incremented by a factor of 0.0001 to ensure that the contour crosses between nodes. Flow charts for the program and subroutine CONTUR are shown in Figs. 11 and 12. CONTOUR PLOTFING EXAMPLES (1) Torsion of square section b a r - stress function plot. One quarter of the section of the bar was analysed using 32 triangular elements and 16 linear quadrilaterals with boundary conditions shown in Fig. 13. Mesh and contour plots are shown in Figs. 14 and 15. (2) Torsion of cylindrical s e c t i o n - stress function plot. Twenty curvilinear elements were used with boundary conditions as shown in Fig. 16. Figure 17 shows mesh and contour plots obtained using the program. (3) Streamline flow around cylinder. Eighty-one triangular elements were used as shown in Figs. 18 and 19. (4) Test using random array of elements - i s o t h e r m plot. This example was run to check generated contours, the example illustrated in Fig. 20 is for heat flow through a wall with prescribed face temperatures. (5) Contour map plot. This example shows the program used to plot contours other than those obtained from finite element analysis. Figure 21 shows 500 ft and 150 ft contours around Carlingford Lough, Northern Ireland. The area was divided into quadrilaterals as shown. CONCLUSIONS It has been shown that the piecewise linear contours produced by the algorithm give an effective representation of finite element output. The piecewise linear representation was used rather than curve fitting techniques to prevent 'meandering' of contours in coarse mesh regions. It can be seen from Figs. 14 and 15, for the stress function contours at the comer of the square section, that mesh design does have some influence on contour shapes. Refinement of the finite element mesh is used in regions of high gradient producing contours of greater accuracy than those found in regions of low gradient. Further refinement of the contours, of the eight-noded element, could be made by taking an average of the element nodal values at the centroid and dividing these elements into eight rather than two triangles. REFERENCE

1 Zienkiewicz, O. C. The Finite Element Method, McGraw-Hill, London, 1977

Adv. Eng. Software, 1983, Vol. 5, No. 1

31