Calculation of Direction Cosines

The main constituent of both the Neumann and Cauchy boundary conditions is the rate of change of the scalar field $\phi$ in the direction of the outward normal to the boundary. This has been denoted

$$\frac{\partial \phi}{\partial n} = ~~\parbox {55 mm}{rate change of $\phi$\ in the direction of the outward normal}$$ (4.13)

In many problems the value of this will be known in some form; in fact in many cases it will be a constant. However there will be occasions where the explicit value of $\partial \phi/\partial n$ is not known but is given as the rates of change of $\phi$ with respect to some set of global axes.

For example consider a two-dimensional problem formulated in a Cartesian coordinate system $(x,y)$. The boundary values of $\phi$ may well be given with respect to this system of axes. Thus on a boundary

$$\left. \begin{array}{c} \displaystyle\frac{\partial \phi}{\p... ...\partial y} = g(x,y) \end{array}\right\} \mbox{~~on~~} \Gamma$$ (4.14)

These are adequate provided the boundary is parallel to one of the coordinate directions.

Figure 4.1: A Simple Finite Element Mesh

In the example region shown in Figure 4.1, clearly along AD ($x=a$)

$$\frac{\partial \phi}{\partial n} = - \frac{\partial \phi}{\partial x}$$ (4.15)

and along DE ($y=d$)

$$\frac{\partial \phi}{\partial n} = + \frac{\partial \phi}{\partial y}$$ (4.16)

However, problems present themselves whenever the boundary is not parallel to the coordinate axes, for example along EF. The normal derivative must now be constructed from the known parts. Clearly this will be

$$\frac{\partial \phi}{\partial n} = l_x \frac{\partial \phi}{\partial x} + l_y \frac{\partial \phi}{\partial y}$$ (4.17)

where $l_x$ and $l_y$ are the directions cosines of the outward normal at some point along EF, given by

$$l_x = {\bf i}.{\bf n} \mbox{~~and~~} l_y = {\bf j}.{\bf n}$$ (4.18)

i and j being the unit vectors associated with the coordinate directions $x$ and $y$. (4.15) and (4.16) can be easily derived from (4.17) and (4.18) and the whole analysis extended to three-dimensional problems.

A straight boundary not parallel with the coordinate axes does not present real difficulty since $l_x$ and $l_y$ will be constants. However in the most general case, where the boundary could be any smooth curve, $l_x$ and $l_y$ will be functions of position. Within the context of finite element analysis the direction cosines will be required on sides of elements.

In general the element in global coordinates is generated from a parent element by means of some isoparametric transformation. The most used transformation is

$$x = \sum_i N_i(\xi ,\eta ,\zeta ) x_i$$ (4.19)

$$y = \sum_i N_i(\xi ,\eta ,\zeta ) y_i$$ (4.20)

$$z = \sum_i N_i( \xi ,\eta , \zeta) z_i$$ (4.21)

where the $N_i$ are a set of shape functions (see Chapter 2) and $(x_i,y_i,z_i)$ are the coordinates of the element nodes in global coordinates. The $N_i$ are defined in terms of a coordinate system natural to the parent element (most often called local coordinates).

As an example, consider the triangular elements in Figure 4.2. Figure 4.2(a) shows the elements in the global coordinate system and Figure 4.2(b) shows the parent element. The boundary $\Gamma$ is approximated by a series of element sides; for each side the outward normal is required.

Figure 4.2: Finite elements fitting a curved boundary

The construction of these direction cosines is based on two results (Rutherford, 1962). Firstly any element side can be presented as part of a set of equipotential surfaces in the $(x,y,z)$ coordinate system. These are characterised by

$$\Psi(x,y,z) = \mbox{~constant}$$ (4.22)

where $\Psi$ is some continuous function possessing continuous partial derivatives.

The second is a result from coordinate geometry:

If $\nabla \Psi= 0$ at any point P on the surface $\Psi = \mbox{~constant}$, then $\nabla \Psi$ is perpendicular to the surface, passing through P,

where

$$\nabla = {\bf i} \frac{\partial}{\partial x} + {\bf j} \frac{\partial}{\partial y} + {\bf k} \frac{\partial}{\partial z}$$ (4.23)

It is clear from this that the direction cosines of the outward normal are proportional to the components of $\nabla \Psi$:

$$l_x \propto \frac{\partial \Psi}{\partial x}, \, l_y \propto... ...}{\partial y}, \, l_z \propto \frac{\partial \Psi}{\partial z}$$ (4.24)

and just as this applies to the elements in the global system $(x,y,z)$ it will also apply to the parent element in the $(\xi ,\eta ,\zeta )$ coordinate system.

In the local coordinate system the direction cosines $(l_{\xi}, l_{\eta} \mbox{~and~} l_{\zeta})$ of the outward normal to any of the sides are simple to calculate. As in the global case these are proportional to the components of $\nabla \psi$, where $\psi$ is the equation of one of the element sides; hence

$$l_{\xi} \propto \frac{\partial \Psi}{\partial \xi}, \, l_{\e... ...ta}, \, l_{\zeta} \propto \frac{\partial \Psi}{\partial \zeta}$$ (4.25)

The chain rule governing partial differentiation reads:

$$\left[ \begin{array}{c} \displaystyle\frac{\partial}{\partia... ... \displaystyle\frac{\partial}{\partial z} \end{array} \right]$$ (4.26)

and using (4.20) the direction cosines of the outward normal to the parent element can be transformed into the outward normal to an element in the global coordinate system.

Only the direction of the normal is given so finally the vector is normalised. The components of this normalised vector will be the direction cosines $l_{\xi}, l_{\eta}$ and $l_{\zeta}$. Within the context of a finite element program (Segment 3.2 is the best example) the code to calculate the direction cosines is contained within the boundary integration loop since the normal direction will be required at the Gauss points on the boundary:

 
      CALL MATMUL(LDER, ILDER, JLDER, GEOM, IGEOM, JGEOM, JAC,
     *     IJAC, JJAC, DIMEN, NODEL, DIMEN, ITEST)
      CALL MATINV(JAC, IJAC, JJAC, JACIN, IJACIN, JJACIN, DIMEN,
     *     DET, ITEST)
      CALL DCSTRI(JACIN, IJACIN, JJACIN, SIDNUM, COSIN, ICOSIN,
     *     ITEST)

The routine DCSTRI uses the current transformation Jacobian and the side number ( SIDNUM) on which the direction cosines are required to make the calculation. The direction cosines are return in COSIN. The code above forms the Jacobian ( JAC) using (2.25) and then inverts it using MATINV. The routines DCSQUA and DCSBRK perform similar operations for rectangular and brick elements.