## 1. Introduction

Computational Engineering Mathematics as applied to topics in Computational Geosciences and Computational Fluid Dynamics, among other themes and topics, is the foundation of computational modelling processes involved with Climate and Atmospheric processes. This chapter reviews the key mathematical underpinnings of computational mathematics used to represent climate, atmospheric and hydrologic related processes. The commonly used numerical and computational approaches of finite difference, finite element, and finite volume have common roots stemming from the mathematical analysis of partial differential equations, ordinary differential equations, vector differential and integral equations, among other topics including solution of large scale matrix systems involving tens of thousands of unknown variables, time‐stepped numerous times a modelling second, resulting in incredible numbers of computations. This chapter focuses on important aspects of matrix solutions and on the vector analysis towards use in development of finite volume and nodal domain numeric solutions of the governing Navier‐Stokes equations. Due to the continuing increasing computational capability forecast according to Moore’s Law, the use of Navier‐Stokes solvers in routine problems is becoming more prevalent, and will possibly be common practice in a few short years.

Methods for developing large scale matrix systems and their solution is a fundamental area of knowledge needed for the successful modelling of atmospheric processes and climate. This particular topic can be investigated in other references focused towards solving large scale matrix systems. Methods for efficiently constructing the matrix system for subsequent solution is also a subject for other references. Rather, in the current chapter, the focus is on the topics of matrix system solution and the existence and uniqueness of solutions. The reduction of the governing partial differential equations into numeric form or into finite element or finite volume form, and then further developed into numeric statements, is another subject dealt with in detail in other references.

In the current chapter, the mathematical underpinnings of vector analysis are reviewed as they are applied in the development of the ensemble of numeric statements for subsequent matrix solution. With the continued advances in computational power, there is increased interest in the field of atmospheric modelling to decrease the computational scale to a micro‐scale. This interest is partially motivated by the ability to solve large scale matrix systems in the number of occasions to enable a small‐scale time advancement to be approximated in a finite‐difference scheme. Solving entire large scale matrix systems several times a modelling second is now computationally feasible. Hence the motivation to increase computational detail by reducing modelling scale.

## 2. The Navier‐Stokes equations in computational engineering mathematics

The mathematical description of fluid processes are embodied in the well‐known Navier‐Stokes equations. These equations relate Newton’s second law of motion to a control volume of fluid. Their origins can be found in the work dated between 1827 and 1845, including the advances made by Saint‐Venant. Although all fluids possess viscous effects, in many problem situations the effects of viscosity can be neglected, reducing the governing equations considerably. The assumption of the fluid being incompressible in the range of temperatures and conditions under study further reduces the governing equations to being applicable to a fluid that is incompressible and having zero viscosity. That is, an ideal fluid. The analysis of ideal fluid flow is considerably simpler to undertake than the full Navier‐Stokes equations, and provides generally a good framework and basis for the analysis of the complete fluid flow mathematical regime. Ideal fluid flow ignores shear stresses in fluids and normal stresses (that is, thermodynamic pressure forces) are the internal force components carried forward into the analysis.

The resulting Euler equations are well known in fluid mechanics and are particularly amenable to solution when applied to flow on a streamline. For steady flow, the Euler equations indicate that flow velocity and pressure sum to as constant value along a horizontal streamline (or neglecting weight of the fluid control volume). Integration of the Euler equation along a streamline for steady flow conditions results in the Bernoulli equation, which is a commonly used relationship in the analysis of fluid flow. The Bernoulli equation is often misapplied by failure to adhere to the assumptions applied in the derivation of that equation. Namely, steady flow conditions, frictionless or shear‐free flow, incompressible fluid flow, and flow on a streamline. A common misapplication is to apply Bernoulli’s equation at two points within a fluid domain that do not lie on the same streamline. The sum of the Bernoulli equation components is called the Bernoulli constant, and has different values depending on the streamline under study. If the Bernoulli constant is constant throughout the entire fluid domain, regardless of streamline, then the flow is said to be irrotational. Irrotational flow can be assessed within a vector field by use of the vector operator Curl, defined in the development below. In particular, the streamline can be mathematically described in terms of vector notation with respect to usual measure of distance along the streamline,

When dealing with fluid flow, particularly atmospheric fluid flow, thermodynamic pressure is a key component. Thermodynamic pressure is also called static pressure. The stagnation pressure is measured under an ideal frictionless process where the fluid is decelerated to zero speed (note that velocity is a vector whereas speed is a scalar term.) The dynamic pressure quantifies the magnitude of the flow velocity. The flow speed is determined given the stagnation pressure and the static pressure at a point. These relationships become impacted by the speed, particularly with higher Mach numbers greater than about 0.3.

Application of the Bernoulli equation finds good use in many problems of high interest involving fluid flow regimes involving contraction or expansion of flows. In all of these applications, however, care is required to preserve the fundamental assumptions employed in the derivation of the Bernoulli equation.

The above discussion focuses on the mass continuity and flow momentum equations embodied in the Navier‐Stokes equations. A similar resulting formulation arises using the energy equation. In this second situation, a control volume is determined that conforms to the stream lines of the flow regime, and steady incompressible flow is considered that is inviscid. The Bernoulli equation then arrives even though different concepts and boundary conditions are employed. Therefore, the Laws of Thermodynamics for steady, incompressible flow along a streamline simplifies to the well‐known Bernoulli equation which is also employed for analysis of mechanical energy fluid flow problems.

In civil engineering hydraulics problems, the energy equation formulation of the Bernoulli equation is typically shown graphically, as the Energy Grade Line or EGL displaying the sum of the three energy head components of pressure head plus velocity head plus gravity head, all in units of length. The sum of these three components is the Total Energy Head, which represents the height the fluid would rise up to in a tube open for the fluid to access. The Hydraulic Grade Line or HGL is plotted the velocity head below the EGL.

## 3. Irrotational flow and the vector operator curl

Another important vector calculus operator is divergence. This operator provides several measures of fluid flow and fluid properties including compressibility, and also a measure of the strength of source or sink conditions. A divergence value that is positive indicates outwards flux from the target point, whereas a negative value indicates flow trends indicative of a sink.

The vector operators of curl and divergence are key vector calculus tools used in the analysis of fluid flows of both compressible and incompressible flow regimes, and form the basis of continuity equations.

Modeling topics currently tend to divide along lines of numerical methods, particularly in leveraging parallelism versus increased processor performance. The partial differential equations describing the atmospheric processes involved still are heavily influenced by processes such as rainfall, convective processes, and others. Many models use grid‐based finite‐difference analogs, whereas other approaches are based upon finite element and finite volume schemes. Because many of these schemes involves grid density problems at the poles, further attention has been paid towards use of discretization's of the globe into different geometric finite volumes and tessellations. Finite volume methods, like finite element models, are attractive in conserving energy, mass, and momentum. Furthermore, the Divergence Theorem has a more direct application with use of finite volume or finite element methods.

In this chapter, the fundamental vector calculus principles are reviewed that are most relevant to fluid flow modeling such as involved with compressible and incompressible flow problems in water and the atmosphere. Details regarding these fields far surpass the scope of this chapter, but the fundamental elements are generally based in the computational mathematics involved in formulation of and solution of large scale matrix equations.

## 4. Vectors

### 4.1. Definition and representation

A vector is used to describe a quantity such as displacement, velocity, or force that has both magnitude and direction. In two dimensions or three dimensions, we can represent the magnitude and direction of a vector with an arrow (**Figure 1**) [1].

By convention, a vector is usually written in bold, or with an arrow, such as **Figure 2**).

It is important to note that a vector is not anchored to any particular place in the coordinate system. Rather, it represents only a *change* in terms of distance and direction. Thus, two vectors are equal if they have the same magnitude and direction (**Figure 3**).

Most often when dealing with vectors of more than two components, it is easier to represent them algebraically. We write a vector as a list of its components, which are equivalent to the change in each of the coordinate directions (**Figure 4**). This is fairly easy to represent in two or even three dimensions, but also extends to

A vector that denotes a specific point in the coordinate system is called a position vector (meaning a vector that, when its tail is placed at the origin, points to a specific point) Thus, a position vector for a point **Figure 5**)

### 4.2. Vector addition and scalar multiplication

Two basic operations that can be performed on vectors are addition and scalar multiplication. In other words, we can add and subtract vectors, as well as multiply them by a constant which will lengthen or shorten the vector.

**Definition:** If **v** are two vectors positioned so the initial point of **sum** **Figure 6**).

Algebraically, vector addition results in adding the vectors component‐wise:

In addition to representing a vector by its components, we can represent a vector by a linear combination of the basis vectors, a basis vector being a vector with a magnitude of one in the direction of only one of the coordinate axes.

Vector addition results in factoring each basis vector:

**Definition:** If

Algebraically, scalar multiplication results in multiplying each component of the vector by the scalar constant:

Vector subtraction results from applying both vector addition and scalar multiplication by

Given

(commutative property)

The length of a vector can be found by utilizing the Pythagorean Theorem. If a vector has the components **u** is found by:

### 4.3. The dot product

Multiplication between vectors is not the same as between scalars. With vectors we have two “products”, the dot product and the cross product. The dot product between vectors will play a vital part in vector calculus. We can think of the dot product as being a metric of how much one vector is influencing another. An application for the dot product is work, where work is defined as **Figure 7**).

**Definition:** If **number**

The dot product of two vectors is therefore a scalar (number), by definition, and not a vector. To find the dot product we multiply the corresponding components and add.

**Theorem:** if

**Important property:** Two non‐zero vectors

### 4.4. The cross product

The cross product is the other “product” between two vectors. The caveat is that the cross product is only defined for vectors in three dimensions. An application for the cross product is torque. This is similar to work which we saw with the dot product, but the cross product results in a vector orthogonal to the two vectors that created it.

**Definition:** If **vector**

The cross product of two vectors is therefore a vector (unlike the dot product). The cross product **orthogonal** to **Figure 8**).

**Theorem:** If

That is, the magnitude of the resultant vector from the cross‐product depends on the angle between them. The direction (while remaining orthogonal) depends on the order of the cross product,

## 5. Parametric equations

An important concept within vector calculus is parameterization. In this we will create a vector function that can build a line, a curve (called a space curve) or a surface in space. The idea is that each component of the vector is a function of a single independent variable called the parameter. The parameter dictates the value of the component for each value in the domain of the parameter. The terminal point of this vector will “draw out” the line, curve or surface as the parameter varies through its domain.

Example:

### 5.1. Equation of a line

To parameterize a line we will create a vector equation in which

The parameter is **Figure 9**.

As **r** “draws” the line. When we look at the resulting components from the vector equation, the parametric equations are:

If we wanted to parameterize a straight line between the points (0,1) to (2,1) then we could use

### 5.2. Equation of a plane

The vector equation for a plane is: **n** is a normal vector to the plane and

## 6. The gradient

An important calculation in vector calculus is the gradient. We will see that the gradient of a function will return a vector that points in the direction of greatest increase from any point in the function. First we will inspect where the gradient comes from.

### 6.1. Partial derivatives

As we move into three dimensions we see our functions represented as

The partial derivatives are found, largely, the same way as single derivatives. The only new aspect is the other variable, which we treat as a constant when taking the derivative.

Example:

### 6.2. Directional derivative

Now that we know how to find the rate of change in either the x or y direction, what about the other infinite directions we may be interested in? We will utilize a directional derivative to find the rate of change from a point in a particular direction.

**Theorem:** if f is a differentiable function of

The directional derivative appears to be a dot product between two vectors. In fact if we assume that it is and write the function as a dot product between two vectors the resulting equation reveals:

We immediately see our vector **u**, and the second vector of partial derivatives is known as the Gradient.

### 6.3. The gradient

When talking about the gradient we will introduce an operator called the dell operator, which is similar to an upside down triangle. The del operator is a vector.

The gradient is the product of the del operator and the function

Thus another way to annotate a directional derivative would be:[1]

It is important to note that the gradient of a function will point in the direction of greatest increase of a function at a specific point. The magnitude of the gradient at that point will be the greatest rate of change. When viewing a function as its level curves. The gradient will be orthogonal to the curve at any point.

Example:

Find the maximum rate of change of

### 6.4. Tangent planes

In a similar way to finding a line tangent to a curve at a point, we can find a plane tangent to a surface at a point. The idea of a tangent plane will be applied later in the chapter.

Suppose that we have a surface

If we think of the differences between the values of

## 7. Vector fields

We saw that a vector has both a magnitude and a direction. When we assign a vector to every point in space we build a vector field. A practical example in the real world would be wind speeds. It is easy to see that at every point the wind has both a velocity and a direction.

**Definition:** Let

### 7.1. Curl

Perhaps, the most important attribute of a vector field is the curl. The curl will determine if there is some sort of twist or spin within the vector field. This will become very important when we start to discuss work through a vector field. The curl is determined by taking the cross product between the del operator and the vector field.

To find the actual curl vector we will find the determinate of the matrix

If the curl results in a zero vector, then the vector field is determined to have zero curl and the vector field is conservative.

When we compare the plots of the two similar vector fields in **Figure 10**, we notice a subtle difference between them. In the vector field on the left there are vectors pointing towards the center as well as vectors pointing away from the center. This indicates curl. In the vector field on the right, all points are pointing away from the center which indicates that there is no curl in the vector field.

### 7.2. Divergence

The divergence of a vector field can be thought of as a quantification of the rate that the vector field is expanding or contracting per unit volume at any point in the field [3]. If the divergence is equal to zero the vector field is said to be incompressible. Divergence is calculated by taking the dot product between the del operator and the vector field [2].

If the divergence at a point is positive the point is considered a source. If the divergence is negative the point is a sink. If the divergence is zero, then the point is neither a sink nor a source.

## 8. Line integrals

If we remember back to the dot product, we defined work as **F** and the tangent vector to the curve at a point between them, **T**(the unit tangent) or *d***r** (**Figure 11**).

If

As

### 8.1. The fundamental theorem for line integrals

We talked earlier about conservative vector fields, where the curl was zero. There is an important implication to the curl being zero; there exists a function such that the gradient of that function is equal to the vector field. This function is called the potential function.

**Theorem:** Let

Clearly, if the path is closed, where the start point is the same as the end point, the work will be zero. Another implication to the conservative vector field is path independence. This means that the work done by a conservative vector field on an object moving between two points is the same regardless of the path.

## 9. Green's Theorem

In essence, Green's Theorem allows us to calculate the work done around a closed path through a non‐conservative vector field by looking at the area enclosed by the path rather than the path itself.

**Green's Theorem:** Let

We can recognize in this integral that we are translating the line integral into a double integral of the curl in the x‐y plane over the area enclosed by the curve. The positive orientation indicates that we are moving in a counter clockwise direction around the curve. Green's Theorem can also be applied to the

Example: [2]

Evaluate the given line integral around the path defined by the circles:

We are going to utilize the area between the two circles to solve this line integral.

By Green's Theorem:

Converting to polar coordinates:

## 10. Surface area

Previously we saw how to parameterize a line in which the *x*, *y* and *z* components were made to be functions of another variable (usually *u* and *v*). This will return a vector equation describing our surface as [5]:

Finding the surface area of a flat or geometric surface is relatively easy as there are known formulas for such cases. However, it is more difficult to find the surface area of a more complexly defined surface. In this case we will see that the relationship between surface area and surface integrals is similar to the relationship between arc length and line integrals.

If we have a space curve defined by the vector function

If the curve is in the x‐y plane where

In a similar sense to the arc length, where we utilized a secant line to approximate the length of a curve, we can utilize a tangent plane to a surface to approximate the surface area.

The primary idea is that to get as close to the true surface area as possible we must use an ever increasing number of planes approximating our surface. Then as the number of planes goes towards infinity the sum of their areas approaches the area of the surface.

It is important to remember how to define a plane. Two ways that will be important to surface area are finding two vectors within the plane and by finding a point with a normal vector to the plane. In order to find two vectors within the plane we can utilize the resulting vectors from the directional derivatives of our surface at our point within the surface. With the parameterization of the surface the vectors would become

In order to find the surface area of the plane we can scale each vector by the change of the parameter of the subinterval. We then take the cross product of the scaled vectors and by the rules for the cross product we can factor out the change in variables such that:

If we then sum the area of each plane, the approximation to the surface becomes:

Then the surface area for our surface S is: [4]

If the surface is defined in terms of

then the formula becomes: [4]

## 11. Surface integrals

In this section we will see that the surface integral is very similar to what a line integral is to curves in space.

It is clear that each coordinate is a function of our two parameters. The parameterized surface,

The first step is to assume the parameter domain, *D*, is an image of a function and has a rectangular shape which has been divided into rectangular sub sections identified as

Our surface is similarly divided into corresponding subsections, each defined by

As the area *m* and *n* go to infinity we have the definition of the surface integral as: [4]

Let's take a second to remember what the *dS* was:

If we remember back to how surface area was calculated utilizing an infinite number of tangent planes, then the change in the surface area *dS* [7].

Remember that *D*. If, for example, our surface is defined by an equation where the parameterization can be defined by

Then the resulting surface integral becomes: [7]

Similar formulas can be easily derived should the surface be defined in a way other than by

## 12. Oriented surfaces

It is important to talk about the “orientation” of a surface. It is easy to see that there are always two vectors that are normal to any surface as they point in exactly opposite directions. We can think about a “positive” orientation being one defined by a normal vector pointing away from the origin and a “negative” orientation being defined by a normal vector pointing towards the origin. In vector calculus, the only surfaces that are used are those that are orientable. A type of surface that is not orientable, such as the Mobius Strip, is not applicable within vector calculus.

## 13. Surface integrals and flux

Suppose that we have an oriented surface, *S*, with a unit normal vector, **n** (remember, that a unit normal vector is a vector orthogonal to the surface with a length of 1 unit. If we think of this surface being within a vector field describing, perhaps, a fluid with a density function of *S*, then the amount of the “fluid” passing through our surface at one point would be the dot product between the vector field at the point and the normal vector to the surface at that same point.

The rate of flow (mass per unit time) per unit area is **n** per unit time by:

If we add this mass for each sub section as the area of each subsection goes to 0 it results in an infinite sum or integral which calculates the rate of flow through *S*:

This type of surface integral occurs frequently in physics. If the density function is constant then we are only left with the velocity vector field.

More simply, if **F** is a continuous vector field defined on an oriented surface *S* with a unit normal vector **n**, then the surface integral of **F** over *S* is: [4]

We are then able to translate this integral into our previously shown double integral over a region such that: [4]

With regards to the actual mechanics of how this integral is calculated it is important to have an understanding of the relationship between the parts. It is also important to pay attention to how the surface is defined. Let's assume

The differential for surface area is:

The resulting integral then becomes:

The denominator within the dot product and the radical within the surface differential will cancel which leaves us with the final integral: [4]

## 14. Stokes Theorem

Stokes’ Theorem is very similar to Green's Theorem. In fact it can be thought of that Green's Theorem is a special case of Stokes’ Theorem, or that Stokes’ Theorem is a higher‐dimensional version of Green's Theorem. What Stokes is going to do is relate a line integral with a surface integral. Stokes’ Theorem will find the work done around a closed path by computing a surface integral of the curl of a vector field over a “surface cap [8]” on top of the path.

**Stokes’ Theorem:** Let *S* be a smooth oriented surface in *C* whose orientation is consistent with that of S. Assume that *S*. Then [9]

In other words the circulation of a vector field around the boundary of an oriented surface in space in the direction counter clockwise with respect to the surface's unit normal vector equals the integral of the curl for vector field over the normal component of the surface. The surface must also be piecewise smooth.

We can think of this integral as relating a line integral around a closed path to a surface integral over a “capping surface” of the path. The interesting idea is that it doesn't matter what the surface is over the path, once the surface integral is calculated it reaches the same value!

While this result is difficult to prove in the general case, an easier intuition is gained through inspection of Green's Theorem demonstrated earlier. This can then be extended to Stokes, which is the same concept in three dimensions.

Example:

Evaluate

First check the curl to see if the vector field is conservative.

We are going to utilize Stokes’ Theorem to solve this problem. Remembering that we can define our surface over the path as anything, let's use the plane

If we project the surface into the x‐y plane it will form a circle. By making the substitution 5 for

Further substituting 5 for

## 15. Divergence Theorem

In this final section we move away from line integrals and again visit flux. In this case we will be looking at how a vector field moves through a solid shape in space.

The divergence theorem states: Let *S* be the boundary surface of **F** be a vector field whose component functions have continuous partial derivatives on an open region that contains *E*. Then [10]

In other words the Divergence Theorem states that the flux of **F** across the boundary surface of **F** over

The Divergence Theorem is to surface integrals what Green's Theorem is to line integrals in that it allows us to convert a surface integral over a closed surface into a triple integral over a closed region.

While the equation looks to be a direct equivalence between the double and triple integrals, it is not quite that simple. If we were to try and utilize the surface integral definition for flux on a closed surface, we would have to compute the flux through each sub surface individually then add the result. In contrast, the divergence theorem allows us to find the flux through the closed surface by utilizing the divergence through the volume.

Example:

Find the outward flux of **F** through the closed surface of the cylinder

If we project our region into the x‐y plane, we see that it is a quarter circle with a radius of 2. In this problem we will utilize cylindrical coordinates to solve the integral.