Solution of Differential Equations with Applications to Engineering Problems

Abstract

Over the last hundred years, many techniques have been developed for the solution of ordinary differential equations and partial differential equations. While quite a major portion of the techniques is only useful for academic purposes, there are some which are important in the solution of real problems arising from science and engineering. In this chapter, only very limited techniques for solving ordinary differential and partial differential equations are discussed, as it is impossible to cover all the available techniques even in a book form. The readers are then suggested to pursue further studies on this issue if necessary. After that, the readers are introduced to two major numerical methods commonly used by the engineers for the solution of real engineering problems.

1. Introduction

1.1. Classification of ordinary and partial equations

To begin with, a differential equation can be classified as an ordinary or partial differential equation which depends on whether only ordinary derivatives are involved or partial derivatives are involved. The differential equation can also be classified as linear or nonlinear. A differential equation is termed as linear if it exclusively involves linear terms (that is, terms to the power 1) of y, y′, y″ or higher order, and all the coefficients depend on only one variable x as shown in Eq. (1). In Eq. (1), if f(x) is 0, then we term this equation as homogeneous. The general solution of non-homogeneous ordinary differential equation (ODE) or partial differential equation (PDE) equals to the sum of the fundamental solution of the corresponding homogenous equation (i.e. with f(x) = 0) plus the particular solution of the non-homogeneous ODE or PDE. On the other hand, nonlinear differential equations involve nonlinear terms in any of y, y′, y″, or higher order term. A nonlinear differential equation is generally more difficult to solve than linear equations. It is common that nonlinear equation is approximated as linear equation (over acceptable solution domain) for many practical problems, either in an analytical or numerical form. The nonlinear nature of the problem is then approximated as series of linear differential equation by simple increment or with correction/deviation from the nonlinear behaviour. This approach is adopted for the solution of many non-linear engineering problems. Without such procedure, most of the non-linear differential equations cannot be solved. Differential equation can further be classified by the order of differential. In general, higher-order differential equations are difficult to solve, and analytical solutions are not available for many higher differential equations. A linear differential equation is generally governed by an equation form as Eq. (1).

$\frac{d^{n}y}{d{x}^n}+a_1\left(x\right)\frac{d^{n-1}y}{d{x}^{n-1}}+\dots +a_n\left(x\right)y=f\left(x\right)$ E1

“Non-linear” differential equation can generally be further classified as

1. Truly nonlinear in the sense that F is nonlinear in the derivative terms.

   $F\left(x_1,x_1,x_n,u,\frac{\partial u}{\partial x_1},\frac{\partial u}{\partial x_2},\frac{{\partial }^{2}u}{\partial x_1\partial x_2}\right)=0$ E2

2. Quasi-linear 1st PDE if nonlinearity in F only involves u but not its derivatives

   $A_1\left(x_1,x_2,u\right)\frac{\partial u}{\partial x_1}+A_2\left(x_1,x_2,u\right)\frac{\partial u}{\partial x_2}=B\left(x_1,x_2,u\right)$ E3

3. Quasi-linear 2nd PDE if nonlinearity in F only involves u and its first derivative but not its second-order derivatives

   $\begin{array}{l}{A}_{11}\left(x_1,x_2,u,\frac{\partial u}{\partial x_1},\frac{\partial u}{\partial x_2}\right)\frac{{\partial }^{2}u}{\partial x_1^2}+A_{12}\left(x_1,x_2,u,\frac{\partial u}{\partial x_1},\frac{\partial u}{\partial x_2}\right)\frac{{\partial }^{2}u}{\partial x_1\partial x_2}\\ +A_{22}\left(x_1,x_2,u,\frac{\partial u}{\partial x_1},\frac{\partial u}{\partial x_2}\right)\frac{{\partial }^{2}u}{\partial x_2^2}=F\left(x_1,x_2,u\frac{\partial u}{\partial x_1},\frac{\partial u}{\partial x_2}\right)\end{array}$ E4

Examples of differential equations:

1. $\frac{dy}{dx}=3x+2$ ; first-order ODE (linear)/nonhomogeneous

2. $\left(y-2x\right)dy-3ydx=0$ ; first-order ODE (nonlinear)/homogenous

3. $\frac{d^{2}y}{d{t}^2}+t^{2}y{\left(\frac{dy}{dt}\right)}^3+y=0$ ; second-order ODE (nonlinear)/homogenous

4. $\frac{d^{4}x}{d{t}^4}+5\frac{d^{2}x}{d{t}^2}+7x=sint$ ; fourth-order ODE (linear)/nonhomogeneous

5. $\frac{\partial z}{\partial x}+\frac{\partial z}{\partial y}=2z;$ first-order PDE (linear)/homogeneous

6. $\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+4x+3y-uz=0$ ; second-order PDE (linear)/nonhomogeneous

7. $x\frac{{\partial }^{2}u}{\partial x^2}+2u\frac{{\partial }^{2}u}{\partial y^2}+3u^2=0$ ; second-order PDE (linear)/homogeneous

8. $\frac{du}{dx}-\frac{dv}{dx}=6x$ ; 1st ODE (linear) for two unknowns/nonhomogeneous

1.2. Typical differential equations in engineering problems

Many engineering problems are governed by different types of partial differential equations, and some of the more important types are given below.

1. Tricomi equation: $y\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=0\{\begin{array}{c}y>0:elliptic\\ y<0:hyperbolic\end{array}$

2. Laplace equation (or variants): $\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial y^2}={\nabla }^2\phi =0$

3. Poisson’s equation: $\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial y^2}=f\left(x,y\right)$

4. Helmholtz equation: $\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial y^2}+c^2\phi =0$

5. Plate bending: ${\nabla }^2{\nabla }^{2}w={\nabla }^{4}w=\frac{q}{D}$

6. Wave equation (1D-3D): $\frac{{\partial }^{2}u}{\partial t^2}-c^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=0$

7. Fourier equation: $\frac{\partial T}{\partial t}=\alpha \left(\frac{{\partial }^{2}T}{\partial x^2}\right)$

There are many methods of solutions for different types of differential equations, but most of these methods are not commonly used for practical problems. In this chapter, the most important and basic methods for solving ordinary and partial differential equations will be discussed, which will then be followed by numerical methods such as finite difference and finite element methods (FEMs). For other numerical methods such as boundary element method, they are less commonly adopted by the engineers; hence, these methods will not be discussed in this chapter.

1.3. Separable differential equations

For equations which can be expressed in separable form as shown below, the solution can be obtained easily as

$\frac{dy}{dx}=F\left(x,y\right)\frac{dy}{\Phi \left(y\right)}=f\left(x\right)dx{\int }^{\text{}}\frac{dy}{\Phi \left(y\right)}={\int }^{\text{}}f\left(x\right)dx+c$ E5

$M\left(x,y\right)dx+N\left(x,y\right)dy=0M\left(x\right)dx=-N\left(y\right)dy$ E6

$\text{then}{\int }^{\text{}}M\left(x\right)dx=-{\int }^{\text{}}N\left(y\right)dy+c$ E7

Example:

$\begin{array}{l}\frac{dy}{dx}=x^3\left(y^2+1\right)\Rightarrow \frac{dy}{y^2+1}=x^{3}dx\\ {\int }^{\text{}}\frac{dy}{y^2+1}={\int }^{\text{}}{x}^{3}dx+c\Rightarrow ta{n}^{-1}y=\frac{1}{4}{x}^4+C\Rightarrow y=tan\left(\frac{1}{4}{x}^4+c\right)\end{array}$ E9000

Example:

$\frac{dy}{dx}=\frac{3x^2+4x+2}{2\left(y-1\right)}$ subject to $y\left(0\right)=-1$

Since this is a separable function, the problem can be solved as

$\begin{array}{l}2\left(y-1\right)dy=\left(3x^2+4x+2\right)dx\\ y^2-2y=x^3+2x^2+2x+c\end{array}$ E9008

Based on the boundary condition, c = 3, hence $y^2-2y=x^3+2x^2+2x+3$ .

This quadratic equation in y ² can be solved with two solutions by the quadratic equation as

$y=1-\sqrt{x^3+2x^2+2x+4}\text{ and }y=1+\sqrt{x^3+2x^2+2x+4}.$ E9001

Since the second solution does not satisfy the boundary condition, it will not be accepted; hence, the solution to this differential equation is obtained.

1.4. Variation of parameters

For the following equation form, it is possible to solve it by variations of parameters.

$\text{For }\frac{dy}{dx}=p\left(x\right)y+Q\left(x\right)$ E8

Put $y=c\left(x\right)e^{{\int }^{\text{}}p\left(x\right)dx}$ . By differentiating, it gives $\frac{dy}{dx}=\frac{dc\left(x\right)}{dx}{e}^{{\int }^{\text{}}p\left(x\right)dx}+{\displaystyle \underset{p\left(x\right)y}{{\displaystyle \underbrace{c\left(x\right)p\left(x\right)e^{{\int }^{\text{}}p\left(x\right)dx}}}}}$ . Substitute it to the original ODE $\frac{dc\left(x\right)}{dx}=Q\left(x\right)e^{-{\int }^{\text{}}p\left(x\right)dx}$ . Comparing the terms, it gives

$c\left(x\right)={\int }^{\text{}}Q\left(x\right)e^{-{\int }^{\text{}}p\left(x\right)dx}dx+{\displaystyle \overline{c}}.$ E9

Example:

$\left(x+1\right)\frac{dy}{dx}-ny=e^x{\left(x+1\right)}^{n+1}$ E9003

This equation is now expressed as

$\begin{array}{l}\frac{dy}{dx}=p\left(x\right)y+Q\left(x\right)\\ \frac{dy}{dx}=\frac{n}{x+1}y+\underset{Q\left(x\right)}{\underbrace{e^x{\left(x+1\right)}^n}}\end{array}$ E9002

For x ≠ −1

Solving the homogeneous part of the ODE

$\frac{dy}{dx}=\frac{n}{x+1}y$ then $\frac{dy}{y}=\frac{n}{x+1}dx$

$\begin{array}{c}\hfill ln\left|y\right|=nln\left|x+1\right|+c_1\hfill \\ \hfill y=c{\left(x+1\right)}^n\hfill \end{array}$ E9005

Look for solution $y=c\left(x\right){\left(x+1\right)}^n$ , where c(x) is the variation of parameters. Substitute it to the ODE

$\begin{array}{c}\hfill \frac{dc\left(x\right)}{dx}{\left(x+1\right)}^n+nc\left(x\right){\left(x+1\right)}^{n-1}=nc\left(x\right){\left(x+1\right)}^{n-1}+e^x{\left(x+1\right)}^n\hfill \\ \hfill \frac{dy}{dx}=\frac{n}{x+1}y+e^x{\left(x+1\right)}^n\hfill \end{array}$ E9006

Comparison gives $\frac{dc\left(x\right)}{dx}=e^x$

Integration of this equation gives $c\left(x\right)=e^x+{\displaystyle \overline{C}}$

General solution is hence given by $y={\left(x+1\right)}^n\left(e^x+{\displaystyle \overline{C}}\right)$

The Bernoulli equation is an important equation type which can be solved in a similar way by variation of parameters. Consider the following form of equation

$\frac{dy}{dx}=p\left(x\right)y+Q\left(x\right)y^n$ E10

$Step 1:\text{ Put }z=y^{1-n}$ E11

$\begin{array}{c}Step 2:\text{ Then }\frac{dz}{dx}=\left(1-n\right)y^{-n}\frac{dx}{dy}\\ \frac{dz}{dx}=\left(1-n\right)P\left(x\right)z+\left(1-n\right)Q\left(x\right)\end{array}$ E12

The non-linear ODE now becomes linear ODE. It can be solved by formula

Step 3: n = −1, z = y ². Inverting z to get y

$\frac{dy}{dx}=\frac{y}{2x}+\frac{x^2}{2y}$ E13

$\frac{dz}{dx}=\frac{1}{x}z+x^2$ E14

$z=e^{\int \frac{1}{x}dx}\left({\displaystyle \int x^2{e}^{-\int \frac{1}{x}dx}}dx+c\right)=cx+\frac{1}{2}{x}^3$ E15

Back substitution of $z=y^2$

$y^2=cx+\frac{1}{2}{x}^3$ E16

1.5. Homogeneous equations

For equation of the following type, where all the coefficients are constant, it can be evaluated according to different conditions.

$\frac{dy}{dx}=\frac{a_{1}x+b_{1}y+c_1}{a_{2}x+b_{2}y+c_2}$ E17

Case 1: $c_1=c_2=0$

$\frac{dy}{dx}=\frac{a_{1}x+b_{1}y}{a_{2}x+b_{2}y}=\frac{a_1+b_1\frac{y}{x}}{a_2+b_2\frac{y}{x}}=g\left(\frac{y}{x}\right)$ E18

Step 1: Set $u=\frac{y}{x}$ , then $\frac{dy}{dx}=x\frac{du}{dx}+u$

Step 2: $\frac{du}{dx}=\frac{g\left(u\right)-u}{x}$ . The resulting non-linear ODE is hence separable and can be solved implicitly.

Step 3: Inverting u to get y.

Case 2: $\left|\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right|=0$

$a_1{b}_2-a_2{b}_1=0$ then $\frac{a_1}{a_2}=\frac{b_1}{b_2}=k$

$\frac{dy}{dx}=\frac{a_{1}x+b_{1}y+c_1}{a_{2}x+b_{2}y+c_2}=\frac{k\left(a_{2}x+b_{2}y\right)+c_1}{a_{2}x+b_{2}y+c_2}=f\left(a_{2}x+b_{2}y\right)$ E19

By change of variables as $u=a_{2}x+b_{2}y$

$\frac{du}{dx}=a_2+b_2\frac{dy}{dx}=a_2+b_{2}f\left(u\right)$ , then

$\frac{du}{a_2+b_{2}f\left(u\right)}=dx$ E20

The resulting non-linear ODE is now separable and can be solved.

Case 3: $\left|\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right|\ne 0$ $c_1\ne 0$ and $c_2\ne 0$

Set $\{\begin{array}{c}{a}_{1}x+b_{1}y+c_1=0\\ a_{2}x+b_{2}y+c_2=0\end{array}$ . Intersecting point of these two lines on xy - plane and (α, β) ≠ 0

$xy-\text{plane and }\left(\alpha ,\beta \right)\ne \left(0,0\right)$ E21

Apply change of variables

$\{\begin{array}{c}X=x-\alpha \\ Y=y-\beta \end{array}\{\begin{array}{c}x=X+\alpha \\ y=Y+\beta \end{array}$ E22

$\begin{array}{c}{a}_{1}x+b_{1}y+c_1=a_1\left(X+\alpha \right)+b_1\left(Y+\beta \right)+c_1=a_{1}X+b_{1}Y+\left(a_1\alpha +b_1\beta +c_1\right)\\ a_{2}x+b_{2}y+c_2=a_2\left(X+\alpha \right)+b_2\left(Y+\beta \right)+c_2=a_{2}X+b_{2}Y+\left(a_2\alpha +b_2\beta +c_2\right)\end{array}$ E23

The original ODE will now become $\frac{dY}{dX}=\frac{a_{1}X+b_{1}Y}{a_{2}X+b_{2}Y}$ which is homogeneous and separable!

Example: $\frac{dy}{dx}=\frac{x+y-1}{x-y+3}$

Solve for $\{\begin{array}{c}x+y-1=0\\ x-y+3=0\end{array}$ we have $\alpha =-1,\beta =2$

Change of variables X = x + 1, Y = y $-$ 2

Then, $\frac{dY}{dX}=\frac{dy}{dx}=\frac{x+y-1}{x-y+3}=\frac{X+Y}{X-Y}=\frac{1+\frac{Y}{X}}{1-\frac{Y}{X}}$

Use a change of variable $u=\frac{Y}{X} X\frac{du}{dX}=\frac{1+u^2}{1-u}\frac{\left(1-u\right)du}{1+u^2}=\frac{dX}{X}$

$\begin{array}{l}\Rightarrow ta{n}^{-1}u-\frac{1}{2}ln\left(1+u^2\right)=ln\left|X\right|+c\\ \Rightarrow ta{n}^{-1}u=ln\left[\sqrt{1+u^2}X\right]+c=ln\left[\sqrt{\left(X^2+Y^2\right)}\right]+c\\ \Rightarrow ta{n}^{-1}\left(\frac{y-2}{x+1}\right)=ln\sqrt{{\left(x+1\right)}^2+{\left(y-2\right)}^2}+c\end{array}$ E9007

There are various tricks to solve the differential equations, like integration factors and other techniques. A very good coverage has been given by Polyanin and Nazaikinskii [29] and will not be repeated here. The purpose of this section is just for illustration that various tricks have been developed for the solution of simple differential equations in homogeneous medium, that is, the coefficients are constants inside a continuous solution domain. The readers are also suggested to read the works of Greenberg [14], Soare et al. [34], Nagle et al. [28], Polyanin et al. [30], Bronson and Costa [4], Holzner [18], and many other published books. There are many elegant tricks that have been developed for the solution of different forms of differential equations, but only very few techniques are actually used for the solution of real life problems.

1.6. Partial differential equations

In many engineering or science problems, such as heat transfer, elasticity, quantum mechanics, water flow and others, the problems are governed by partial differential equations. By nature, this type of problem is much more complicated than the previous ordinary differential equations. There are several major methods for the solution of PDE, including separation of variables, method of characteristic, integral transform, superposition principle, change of variables, Lie group method, semianalytical methods as well as various numerical methods. Although the existence and uniqueness of solutions for ordinary differential equation is well established with the Picard-Lindelöf theorem, but that is not the case for many partial differential equations. In fact, analytical solutions are not available for many partial differential equations, which is a well-known fact, particularly when the solution domain is nonregular or homogeneous, or the material properties change with the solution steps.

1.6.1. Classification of second-order PDE

Refer to the following general second-order partial differential equation:

$A\frac{{\partial }^{2}u}{\partial x^2}+B\frac{{\partial }^{2}u}{\partial x\partial y}+C\frac{{\partial }^{2}u}{\partial y^2}+D\frac{\partial u}{\partial x}+E\frac{\partial u}{\partial y}+Fu+G=0$ E24

To begin with, let us consider a review of conic curves (ellipse, parabola and hyperbola)

$A{x}^2+Bxy+C{y}^2+Dx+Ey+F=0$ E25

The conic curve can be classified with the following criterion.

$B^2-4AC=\{\begin{array}{c}>0 \text{hyperbola}\\ =0 \text{parabola}\\ <0 \text{ellipse}\end{array}$ E26

Following the conic curves, the general partial differential is also classified according to similar criterion as

$\text{Classification}\{\begin{array}{c}{B}^2-4AC>0:\text{elliptic}\\ B^2-4AC=0:\text{parabolic}\\ B^2-4AC<0:\text{hyperbolic}\end{array}$ E27

This classification was proposed by Du Bois-Reymond [41] in 1839. In this section, only some of the more common techniques are discussed, and the readers are suggested to read the works of Hillen et al. [16], Salsa [33], Polyanin and Zaitsev [31], Bertanz [2], Haberman [15] and many other published texts.

1.7. Parabolic type: heat conduction/soil consolidation/diffuse equation

The following equation form is commonly found in many engineering applications.

${\alpha }^2\frac{{\partial }^{2}u}{\partial x^2}=\frac{\partial u}{\partial t},0<x<L,t>o$ E28

Initial condition: $u\left(x,0\right)=f\left(x\right),0\le x\le L$

Boundary condition: $u\left(0,t\right)=0,u\left(L,t\right)=0,t>0$

α ² is a constant known as the thermal diffusivity or coefficient of consolidation. For soil consolidation problem, the governing conditions are given by

Initial excess pore pressure

$\begin{array}{l}{u}_e\left(z,0\right)=u_i\left(z\right),0\le z\le 2d\\ u_e\left(0,t\right)=0,u_e\left(2d,t\right)=0,t>0\end{array}$ E29

Drained boundary

$\begin{array}{l}{\alpha }^2{u}_{xx}=u_t,0<x<L,t>0\\ u\left(0,t\right)=0,u\left(L,t\right)=0,t>0\\ u\left(x,0\right)=f(x),0\le x\le L\end{array}$ E30

Assuming variable u(x, t) can be separated, using separation of variables

$u\left(x,t\right)=X\left(t\right)T\left(t\right)$ E31

$\begin{array}{l}{\alpha }^2{X}^{″}T=X{T}^{\prime }\hfill \\ \frac{X^{″}}{X}=\frac{1}{{\alpha }^2}\frac{T^{\prime }}{T}\hfill \\ \frac{X^{″}}{X}=\frac{1}{{\alpha }^2}\frac{T^{\prime }}{T}=-\lambda \to \{\begin{array}{c}{X}^{″}+\lambda X=0\\ T^{\prime }+{\alpha }^2\lambda T=0\end{array}\hfill \end{array}$ E32

A PDE now becomes two ODE which can be solved readily. Based on the boundary condition $u\left(0,t\right)=0,u\left(L,t\right)=0,t>0$

$\begin{array}{r}u\left(0,t\right)=X\left(0\right),T\left(t\right)=0\\ X\left(0\right)=0,X\left(L\right)=0\\ X^{″}+\lambda X=0,X(0)=0,X(L)=0\end{array}$ E33

This is an eigenvalue problem which has solution only for certain λ. The eigenvalues are given by

${\lambda }_n=\frac{n^2{\pi }^2}{L^2},n=1,2,3,\dots$ E34

Hence the eigenfunctions are expressed as

$X_n\left(x\right)=sin\left(\frac{n\pi x}{L}\right),n=1,2,3\dots$ E35

For the time-dependent function T,

$T^{\text{'}}+{\alpha }^2\lambda T=0$ E36

$\begin{array}{l}\frac{dT}{T}=-{\alpha }^2\lambda dt\\ ln\left|T\right|=\frac{-{\alpha }^2{n}^2{\pi }^{2}t}{L^2}+C\end{array}$ E37

hence $T_n=k_n{e}^{-{\left(n\pi \alpha /L\right)}^{2}t}, k_n$ constant. The fundamental solutions are then expressed as

$u\left(x,t\right)=e^{-{\left(n\pi \alpha /L\right)}^{2}t}sin\left(\frac{n\pi x}{L}\right),n=1,2,3\dots$ E38

The Fourier series expansion in x is given by

$u\left(0,t\right)=f\left(x\right),0\le x\le L$ E39

$u\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{u}_n\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{e}^{-{\left(n\pi \alpha /L\right)}^{2}t}sin\left(\frac{n\pi x}{L}\right)$ E40

Initial condition is given as

$u\left(x,0\right)=f\left(x\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_{n}sin\left(\frac{n\pi x}{L}\right)$ E41

$\begin{array}{l}{\displaystyle {{\displaystyle \int }}_0^L}f\left(x\right)sin\left(\frac{m\pi x}{L}\right)dx={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{\displaystyle {{\displaystyle \int }}_0^L}sin\left(\frac{m\pi x}{L}\right)sin\left(\frac{n\pi x}{L}\right)dx\\ {\displaystyle {{\displaystyle \int }}_0^L}f\left(x\right)sin\left(\frac{n\pi x}{L}\right)dx=c_n{\displaystyle {{\displaystyle \int }}_0^L}si{n}^2\left(\frac{m\pi x}{L}\right)dx=c_n\frac{L}{2}\end{array}$ E8000

Solution of the soil consolidation equation is hence given by

$u\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{e}^{-{\left(n\pi \alpha /L\right)}^{2}t}sin\left(\frac{n\pi x}{L}\right)$ E42

$c_n=\frac{2}{L}{\displaystyle {{\displaystyle \int }}_0^L}f\left(x\right)sin\left(\frac{n\pi x}{L}\right)dx(\text{EulerFourier formulas})$ E43

1.8. One-dimensional wave equation

One-dimensional (1D) wave equation appears in many physical and engineering problems. For example, a vibrating string or pile driving process is given by this type of differential equation. This problem is also commonly solved by the method of separation of variables

$\begin{array}{l}{a}^2{u}_{xx}=u_{tt},0<x<L,t>0\\ u\left(0,t\right)=0,u\left(L,t\right)=0,t\ge 0\\ u\left(x,0\right)=f(x),u\left(x,0\right)=0,0\le x\le L\end{array}$ E44

Consider u(x, t) is given by X(x)T(t). The wave equation will give

$\frac{X^{″}}{X}=\frac{1}{a^2}\frac{T^{\prime }}{T}=-\lambda \to \{\begin{array}{c}{X}^{″}+\lambda x=0\hfill \\ T^{\prime }+a^2\lambda t=0\end{array}$ E45

The partial differential equation will then be given by two equivalent ODEs.

$\begin{array}{l}{u}_t\left(x,0\right)=X\left(x\right)T^{\text{'}}\left(0\right)=0,0\le x\le L\to T^{\text{'}}\left(0\right)=0\\ u\left(0,t\right)=X\left(0\right)T\left(t\right)=0,u\left(L,t\right)=X\left(L\right)T\left(t\right) 0\le x\le L\to T^{\text{'}}\left(0\right)=0\end{array}$ E46

$X^{″}+\lambda X=0,X(0)=X(L)=0$ E47

$X_n\left(x\right)=sin\left(\frac{n\pi x}{L}\right),n=1,2,3,\dots$ E48

${\lambda }_n=\frac{n^2{\pi }^2}{L^2},n=1,2,3,\dots$ E49

For the time-dependent function T,

$T^{\text{'}}+a^2\lambda T=0$ E50

$\begin{array}{c}{T}^{\text{'}}\left(0\right)=0 {\lambda }_n=n\pi /L\\ \text{Then }T\left(t\right)=k_{1}cos\left(n\pi at/L\right)-k_{2}sin\left(n\pi at/L\right)\end{array}$ E51

Since $T^{\text{'}}\left(0\right)=0k_2=0$

Therefore, $T\left(t\right)=k_{1}cos\left(n\pi at/L\right)$

Fundamental solution is given by

$u_n\left(x,t\right)=sin\left(\frac{n\pi x}{L}\right)cos\left(\frac{n\pi at}{L}\right),n=1,2,3\dots$ E52

The general solution is then given by

$u\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{u}_n\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_{n}sin\left(\frac{n\pi x}{L}\right)cos\left(\frac{n\pi at}{L}\right)$ E53

Applying the boundary condition

$\begin{array}{l}u\left(x,0\right)=f\left(x\right),0\le x\le L\\ u\left(x,0\right)=f\left(x\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_{n}sin\left(\frac{n\pi x}{L}\right)\to c_n=\frac{2}{L}{\displaystyle {{\displaystyle \int }}_0^L}f\left(x\right)sin\left(\frac{n\pi x}{L}\right)dx\end{array}$ E54

The final solution is then given by

$u\left(x,t\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_{n}sin\left(\frac{n\pi x}{L}\right)cos\left(\frac{n\pi at}{L}\right)$ E55

$c_n=\frac{2}{L}{\displaystyle {{\displaystyle \int }}_0^L}f\left(x\right)sin\left(\frac{n\pi x}{L}\right)dx$ E56

1.9. Laplace equation

Laplace equation forms an important governing condition for many types of problems. Some of the more common forms are given by

1. three-dimensional Laplace equation $u_{xx}+u_{yy}+u_{zz}=0$

2. two-dimensional heat conduction ${\alpha }^2\left(u_{xx}+u_{yy}\right)=u_t$

3. two-dimensional seepage problem $\left(k_x{u}_{xx}+k_y{u}_{yy}\right)=0$

There are two major types of boundary conditions to this problem:

1. Dirichlet problem: boundary conditions prescribed as u

2. Neumann problem: normal derivative u _x or u _y are usually prescribed on the boundary for many mathematical problems. This case can be solved by the use of complex analysis or series method for which many analytical solutions are available in the literature. In many anisotropic seepage problems, however, the normal of a derived quantity at any arbitrary direction (seepage flow normal to an impermeable surface) is 0 instead of u _x or u _y being zero. For such cases, it is very difficult to obtain the analytical solution if the solution domain is nonhomogeneous, and the use of numerical method such as the finite element method appears to be indispensable.

Consider the given Laplace equation, using separation of variables for the analysis.

$\begin{array}{l}{u}_{xx}+u_{yy}=0,0<x<a,0<y<b\\ u\left(x,0\right)=0,u\left(x,b\right)=0,0<x<a\\ u\left(0,y\right)=0,u\left(a,y\right)=f\left(y\right),0<y\le b\end{array}$ E57

Using separation of variables, $u\left(x,t\right)=X\left(x\right)Y\left(y\right)$

$\begin{array}{l}{X}^{″}Y+X{Y}^{″}=0\\ \frac{X^{″}}{X}=-\frac{Y^{″}}{Y}=\lambda \to \begin{array}{c}{X}^{″}-\lambda X=0\\ Y^{\prime }+\lambda Y=0\end{array}\end{array}$ E58

$u_{xx}+u_{yy}=0,0<x<a,0<y<b$ E59

$\begin{array}{l}u\left(x,0\right)=0,u\left(x,b\right)=0,0<x<a\\ u\left(0,y\right)=0,u\left(a,y\right)=f\left(y\right),0<y\le b\end{array}$ E60

$\begin{array}{l}u\left(0,y\right)=X\left(0\right)Y\left(y\right)=0,0<y<b\to X\left(0\right)=0,\\ u\left(x,0\right)=X\left(x\right)Y\left(0\right)=0,0<x<a\to Y\left(0\right)=0,\\ u\left(x,b\right)=X\left(x\right)Y\left(b\right)=0,0<x<a\to Y\left(b\right)=0,\end{array}$ E61

$\begin{array}{c}{X}^{″}-\lambda X=0,X(0)=0\\ Y^{″}+\lambda Y=0,Y(0)=0,Y(b)=0\end{array}$ E62

${\lambda }_n=\frac{n^2{\pi }^2}{b^2},Y_n\left(y\right)=sin\left(\frac{n\pi y}{b}\right),n=1,2,3,\dots$ E63

$X^{″}-\lambda X=0$ , hence $X\left(x\right)=k_{1}cosh\left(n\pi x/b\right)-k_{2}sin\left(n\pi x/b\right)$

Since X(0) = 0, k ₁ = 0

$X\left(x\right)=k_{2}sinh\left(\frac{n\pi x}{b}\right)$ E64

$u_n\left(x,y\right)=sinh\left(\frac{n\pi x}{b}\right)sin\left(\frac{n\pi y}{b}\right)n=1,2,3\dots$ E65

$\begin{array}{l}u\left(a,y\right)=f\left(y\right),0\le y\le b\\ u\left(x,y\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_n{u}_n\left(x,y\right)={\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{c}_{n}sin\left(\frac{n\pi x}{b}\right)cos\left(\frac{n\pi y}{b}\right)\end{array}$ E66

Based on the Fourier expansion as given by

$\begin{array}{l}{\displaystyle {{\displaystyle \int }}_0^b}f\left(y\right)sin\left(\frac{m\pi y}{b}\right)dy={\displaystyle {{\displaystyle \sum }}_{n=1}^{\infty }}{c}_{n}sinh\left(\frac{n\pi a}{b}\right){\displaystyle {{\displaystyle \int }}_0^b}sin\left(\frac{m\pi y}{b}\right)sin\left(\frac{n\pi y}{b}\right)dy\\ {\displaystyle {{\displaystyle \int }}_0^b}f\left(x\right)sin\left(\frac{n\pi x}{b}\right)dx=sinh\frac{m\pi a}{b}{c}_n{\displaystyle {{\displaystyle \int }}_0^b}si{n}^2\left(\frac{n\pi x}{b}\right)dx=sinh\frac{m\pi a}{b}{c}_n\frac{b}{2}\\ u\left(a,y\right)=f\left(y\right)={\displaystyle {{\displaystyle \sum }}_{n=1}^{\infty }}{c}_{n}sinh\left(\frac{n\pi a}{b}\right)sin\left(\frac{n\pi y}{b}\right)\end{array}$ E67

$\begin{array}{l}{c}_{n}sinh\left(\frac{n\pi a}{b}\right)=\frac{2}{b}{\displaystyle {{\displaystyle \int }}_0^b}f\left(y\right)sin\left(\frac{n\pi y}{b}\right)dy\\ c_n=\frac{2}{b}sinh{\left(\frac{n\pi a}{b}\right)}^{-1}{\displaystyle {{\displaystyle \int }}_0^b}f\left(y\right)sin\left(\frac{n\pi y}{b}\right)dy\end{array}$ E9009

1.10. Introduction to numerical methods

In general, analytical solutions are not available for most of the practical differential equations, as regular solution domain and homogeneous conditions may not be present for practical problems. Moreover, the solution domain may be indeterminate (free surface seepage flow), the displacement is large so that the solution may deform under motion, or in an extreme case part of the material may tear off from the main body with continuous formation and removal of contacts. Many engineering problems fall into such category by nature, and the use of numerical methods will be necessary. Currently, there are several major numerical methods commonly used by the engineers: finite difference method, finite element method, boundary element method and distinct element. There are also other less common numerical methods available for practical problems, and many researchers also try to combine two or even more fundamental numerical methods so as to achieve greater efficiency in the analysis. In general, the solution domain is discretized into series of subdomains with many degrees of freedom. The number of variables or degrees of freedom may even exceed millions for large-scale problems, and sometimes very special material properties are encountered so that the system is highly sensitive to the method of discretization and the method of solution. Similar to the ODE and PDE, it is impossible to discuss the details of all the numerical methods and the author choose to discuss the finite element method due to the wide acceptance of the method and this method is more suitable for general complicated methods.

Except for some simple problems with regular geometry and loading, it is very difficult to solve most of the boundary value problems with the yield of analytical solutions. Towards this, the use of numerical method seems indispensable, and the finite element is one of the most popular methods used by the engineers [32, 38]. There are two fundamental approaches to FEM, which are the weighted residual method (WRM) and variational principle, but there are also other less popular principles which may be more effective under certain special cases. In finite element analysis of an elastic problem, solution is obtained from the weak form of the equivalent integration for the differential equations by WRM as an approximation. Alternatively, different approximate approaches (e.g. collocation method, least square method and Galerkin method) for solving differential equations can be obtained by choosing different weights based on the WRM and the Galerkin method appears to be the most popular approach in general.

Specifically, in elasticity for instance, the principle of virtual work (including both principle of virtual displacement and virtual stress) is considered to be the weak form of the equivalent integration for the governing equilibrium equations. Furthermore, the aforementioned weak form of equivalent integration on the basis of the Galerkin method can also be evolved to a variation of a functional if the differential equations have some specific properties such as linearity and selfadjointness. Principles of minimum potential energy and complementary energy are two variational approaches equivalent to the fundamental equations of elasticity.

Since displacement is usually the basic unknown quantity in FEM, only the principle of virtual displacement and minimum potential energy will be introduced in the following section. In this case, the FEM introduced herein is also called displacement finite element method (DFEM). There are other ways to form the basis of FEM with advantages in some cases, but these approaches are less general and will not be discussed here.

1.11. Principle of virtual displacement

The principle of virtual displacement is the weak form of the equivalent integration for equilibrium equations and force boundary conditions. Given the equilibrium equations and force boundary conditions in index notation,

${\sigma }_{ij,j}+f_i=0,\left(\text{in domain} V\right)$ E68

${\sigma }_{ij}{n}_j-T_i=0,\left(\text{on domain boundary }{S}_{\sigma }\right)$ E69

In WRM, without loss of generality, the variation of true displacement $\delta u_i$ and its boundary value (i.e. $-\delta u_i$ ) can be selected as the weight functions in the equivalent integration

${\displaystyle {\int }_V\delta u_i\left({\sigma }_{ij,j}+f_i\right)dV-{\displaystyle {\int }_{S_{\sigma }}\delta u_i\left({\sigma }_{ij}{n}_j-T_i\right)dS=0}}$ E70

The weak form of Eq. (70) is given as

${\displaystyle {\int }_V\left(-\delta {\epsilon }_{ij}{\sigma }_{ij}+\delta u_i{f}_i\right)dV+{\displaystyle {\int }_{S_{\sigma }}\delta u_i{T}_i}}dS=0$ E71

It can be seen clearly from Eq. (71) that the first item in the volume integral indicates the work done by the stresses under the virtual strain (i.e. internal virtual work), while the remaining items indicate the work done by the body force and surface force under the virtual displacement (i.e. external virtual work). In other words, the summation of the internal and external virtual works is equal to 0, which is called the principle of virtual displacement. Under this case, we can conclude that a force system will satisfy the equilibrium equations if the summation of the work done by it under any virtual displacement and strain is equal to 0.

1.12. Principle of minimum potential energy (PMPE)

Based on Eq. (71), we can deduce that

${\displaystyle {{\displaystyle \int }}_V}\left(\delta {\epsilon }_{ij}{D}_{ijkl}{\epsilon }_{kl}-\delta u_i{f}_i\right)dV+{\displaystyle {\int }_{S_{\sigma }}}\delta u_i{T}_{i}dS=0$ E72

Due to the symmetry of the constitutive matrix $D_{ijkl},$ we can further obtain

$\left(\delta {\epsilon }_{ij}\right)D_{ijkl}{\epsilon }_{kl}=\delta \left(\frac{1}{2}{D}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}\right)=\delta U\left({\epsilon }_{mn}\right)$ E73

where $U\left({\epsilon }_{mn}\right)$ is the unit volume strain energy. Given the assumptions in linear elasticity

$-\delta \varphi \left(u_i\right)=f_i\delta u_i,-\delta \psi \left(u_i\right)=T_i\delta u_i$ E74

Eq. (72) is further simplified to

$\delta {\Pi }_P=0$ E75

${\Pi }_P$ is the total potential energy of the system, which is equal to the summation of the potential energy of deformation and external force and can be expressed as

${\Pi }_P={\Pi }_P\left(u_i\right)={\displaystyle {{\displaystyle \int }}_V}\left(\frac{1}{2}{D}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}-f_i{u}_i\right)dV-{\displaystyle {{\displaystyle \int }}_{S_{\sigma }}}{T}_i{u}_{i}dS$ E76

Eq. (75) shows that, among all the potential displacements, the total potential energy of system will take stationary value at the real displacement, and it can be further verified that this stationary value is exactly the minimum value which is the principle of minimum potential energy.

1.13. General expressions and implementation procedure of FEM

The solution of a general continuum problem by FEM always follows an orderly step-by-step process which is easy to be programmed and used by the engineers. For illustration, a three-node triangular element for plane problems is taken as an example to illustrate the general expressions and implementation procedures of FEM.

1.13.1. Discretization of domain

The first step in the finite element method is to divide the structure or solution region into subdivisions or elements. Hence, the structure is to be modelled with suitable finite elements. In general, the number, type, size, and arrangement of the elements are critical towards good performance of the numerical analysis. A typical discretization with three-node triangular element is shown schematically in Figure 1 .

Figure 1.

Discretization of a two-dimensional domain with three-node triangular element.

Mesh generation can be a difficult process for a general irregular domain. If only triangular element is to be generated, this is a relatively simple work, and many commercial programs can perform well in this respect. There are also some public domain codes (EasyMesh or Triangle written in C) which are sufficient for normal purposes. For quadrilateral or higher elements, mesh generation is not that simple, and it is preferable to rely on the use of commercial programs for such purposes.

1.13.2. Interpolation or displacement model

As can be seen from Figure 1(b) , the nodal number of a typical three-node triangular element is coded in anticlockwise order (i.e. in the order of i, j and m), and each node has two degrees of freedom (DOFs) or two displacement components which is stored in a column vector in index notation as follows:

$a_i=\left[\begin{array}{c}{u}_i\\ v_i\end{array}\right]\left(i,j,m\right)$ E77

Totally, each element has six nodal displacements, i.e. six DOFs. Putting all the displacements in a column vector, we can obtain the element nodal displacement column matrix as

$a^e=\left[\begin{array}{c}{a}_i\\ a_j\\ a_m\end{array}\right]={\left[\begin{array}{cccccc}{u}_i& v_i& u_j& v_j& u_m& v_m\end{array}\right]}^T$ E78

In FEM, a nodal displacement is chosen as the basic unknowns, so interpolation at any arbitrary point is based on the three nodal displacements of each element, which is called a displacement mode. For a three-node triangular element, linear polynomial is utilized, and the element displacement in both x -direction and y-direction are

$u={\beta }_1+{\beta }_{2}x+{\beta }_{3}y$ E79

$v={\beta }_4+{\beta }_{5}x+{\beta }_{6}y$ E80

Obviously, displacements of all the three nodes should satisfy Eqs. (79) and (80). By substituting the six nodal displacement components into these equations, it is easy to obtain another form of displacement mode as

$u=N_i{u}_i+N_j{u}_j+N_m{u}_m$ E81

$v=N_i{v}_i+N_j{v}_j+N_m{v}_m$ E82

where

$N_i=\frac{1}{2A}\left(a_i+b_{i}x+c_{i}y\right)\left(i,j,m\right).$ E83

In Eq. (81), $N_i,N_j \text{and} N_m$ denote the interpolation function or shape function for the three nodes, respectively. A is the area of the element, and $a_i,b_i,c_i\cdots ,c_m$ are constants related to the coordinates of the three nodes. Similarly, Eqs. (81) and (82) can also be expressed in the form of matrix as

$u=\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{ccc}\begin{array}{cc}{N}_i& 0\end{array}& \begin{array}{cc}{N}_j& 0\end{array}& \begin{array}{cc}{N}_m& 0\end{array}\\ \begin{array}{cc}0& N_i\end{array}& \begin{array}{cc}0& N_j\end{array}& \begin{array}{cc}0& N_m\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{u}_i\\ v_i\end{array}\\ \begin{array}{c}{u}_j\\ v_j\end{array}\\ \begin{array}{c}{u}_m\\ v_m\end{array}\end{array}\right]=N{a}^e$ E84

where $N$ is the shape function matrix and $a^e$ is the element nodal displacement vector. For the geometric equations, element strains are

$\begin{array}{c}\epsilon =\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right]=Lu=LN{a}^e=L\left[\begin{array}{ccc}{N}_i& N_j& N_m\end{array}\right]a^e\\ =\left[\begin{array}{ccc}{B}_i& B_j& B_m\end{array}\right]a^e=B{a}^e\end{array}$ E85

where $L$ is the differential operator and $B$ is the element strain displacement matrix which can be given as

$B_i=L{N}_i=\left[\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]\left[\begin{array}{cc}{N}_i& 0\\ 0& N_i\end{array}\right]=\left[\begin{array}{cc}\frac{\partial N_i}{\partial x}& 0\\ 0& \frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial y}& \frac{\partial N_i}{\partial x}\end{array}\right]\left(i,j,m\right)$ E86

Substitute Eq. (85) by the stress-strain relation,

$\sigma =\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]=D\epsilon =DB{a}^e=S{a}^e$ E87

where

$S=DB=D\left[\begin{array}{ccc}{B}_i& B_j& B_m\end{array}\right]=\left[\begin{array}{ccc}{S}_i& S_j& S_m\end{array}\right]$ E88

S is called the element stress matrix. It should be noted that both the strain and stress matrices are constant for each element, because in a three-node triangular element, the displacement mode is a first-order function, and differentiating this function will give a constant function.

1.13.3. Stiffness equilibrium equation (SEE) of FEM derived from PMPE

For elastic plane problems, the total potential energy ${\Pi }_P$ in Eq. (76) can be expressed in matrix formulation as follows:

${\displaystyle {\prod }_P=}{\displaystyle {\int }_{\Omega }\frac{1}{2}{\epsilon }^{T}D\epsilon t}dxdy-{\displaystyle {\int }_{\Omega }{u}^{T}ft}dxdy-{\displaystyle {\int }_{S_{\sigma }}{u}^{T}Tt}dS$ E89

where t, f, and T denote the thickness, body force and surface force, respectively. For an FEM problem, the total potential energy is the summation of that from all the elements. Therefore, substituting Eqs. (84) and (85) into Eq. (89) gives

$\begin{array}{l}{\Pi }_P={\displaystyle \sum }_e{\Pi }_P^e={\displaystyle \sum }_e\left(a^e{}^T\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}\frac{1}{2}{B}^{T}DBtdxdy{a}^e\right)\\ -{\displaystyle \sum }_e\left(a^e{}^T\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{N}^{T}ftdxdy\right)-{\displaystyle \sum }_e\left(a^e{}^T\underset{S_{\sigma }{}^e}{\overset{}{{\displaystyle \int }}}{N}^{T}Ttdxdy\right)\end{array}$ E90

Eq. (90) can be viewed as

$\begin{array}{l}{K}^e={\displaystyle {{\displaystyle \int }}_{{\Omega }_e}}{B}^{T}DBtdxdy,P_f^e={\displaystyle {{\displaystyle \int }}_{{\Omega }_e}}{N}^{T}ftdxdy\\ P_S^e={\displaystyle {{\displaystyle \int }}_{S_{\sigma }{}^e}}{N}^{T}Ttdxdy,P^e=P_f^e+P_S^e\end{array}$ E91

where $K^e$ and $P^e$ are named as the element stiffness matrix and equivalent element nodal load matrix, respectively. Substitute Eq. (91) to Eq. (90), the total potential energy of the structure can be simplified as

${\Pi }_P=a^T\frac{1}{2}{\displaystyle {\displaystyle \sum }_e}\left(K^e\right)a-a^T{\displaystyle {\displaystyle \sum }_e}\left(P^e\right)$ E92

Given

$K={\displaystyle {\displaystyle \sum }_e}{K}^e,P={\displaystyle {\displaystyle \sum }_e}{P}^e$ E93

Eq. (92) is further simplified as

${\Pi }_P=\frac{1}{2}{a}^{T}Ka-a^{T}Pa$ E93a

where $K$ and $P$ are global stiffness matrix and global nodal load matrix, respectively.

For PMPE, the variation of ${\Pi }_P \text{is}$ equal to 0 and the unknown variable is $a$ , thus Eq. (75) gives

$\frac{\partial {\Pi }_P}{\partial a}=0$ E94

which finally comes to the SEE of FEM as

$Ka=P$ E95

From Eq. (93), we know that the global stiffness matrix and the global load matrix are the assemblage of the element stiffness matrices and equivalent element nodal load matrices, respectively. Specifically, in order to solve Eq. (93), element stiffness matrix, element equivalent nodal load vector, global stiffness matrices and global nodal load vector are all determined together with some given displacement boundary conditions. Without the provision of adequate boundary condition, the system is singular as rigid body motion will produce no stress in the system and such mode will be present in the SEE.

1.13.4. Derivation of element stiffness matrices (ESM)

For a three-node triangular element, the element strain matrix B is constant, thus Eq. (91) gives

$K^e=B^{T}DBtA=\left[\begin{array}{ccc}{K}_{ii}& K_{ij}& K_{im}\\ K_{ji}& K_{jj}& K_{jm}\\ K_{mi}& K_{mj}& K_{mm}\end{array}\right]$ E96

of which the submatrix

$K_{ij}=\left[\begin{array}{cc}{k}_{ij}{}^{xx}& k_{ij}{}^{xy}\\ k_{ij}{}^{yx}& k_{ij}{}^{yy}\end{array}\right]$ E97

$K_{ij}$ indicates the ith nodal force along the x- and y-directions in the Cartesian coordinate system when the displacement of the jth node is unit along the x- and y-directions, which can be easily obtained. Moreover, the element stiffness matrix is symmetric, and the computational memory required in an FEM program can be reduced by using this property.

It should be noted that for a higher order triangular element (e.g. six-node triangular element) or quadrilateral element for which higher order terms are involved, the strain matrix B is not constant any more so that the element stiffness matrix needs to be evaluated by numerical integration (direct integration is seldom adopted). Towards this, numerical integration methods such as the Gaussian integration or the Newton-Cotes integration can be utilized.

1.13.5. Assembling of ESMs and ENLMs

For an FEM process, we need to solve Eq. (95) which is the global equilibrium equation. Most of the elements in the matrix $K$ are 0 simply because each node is only shared by a few surrounding elements. In view of that, a rectangular matrix can represent the global stiffness matrix (which is a square matrix), and the half bandwidth D can be defined as

$D=\left(1+NDIF\right)\times NDOF$ E98

where NDIF denotes the largest absolute difference between the element node numbers among all the elements in the finite element mesh.

In conclusion, the properties of the global stiffness matrix can be summarized as: symmetric, banded distribution, singularity and sparsity. Among all the properties, singularity will vanish by introducing appropriate boundary conditions to Eq. (95) to eliminate the rigid body motion. Also, other properties like banded distribution should be fully taken into consideration to reduce the computational memory and enhance the computation efficiency.

1.13.6. Isoparametric element and numerical integration

Most of the engineering structure is not regular in shape, and some of them even have very complicated boundary shapes. Although the use of triangular element can always fit a complicated boundary, the accuracy of this element is low in general. To cope with the irregular boundary shape with a higher accuracy in analysis, one of the most common approaches is the use of higher-order element, and the isoparametric formulation is the most commonly used at present. Consider an arbitrary four-node quadrilateral element as an example which is schematically shown in Figure 2 . If we can find the transformation from Figure 2(a) to (b), then it will become easier to carry numerical integration with complicated shapes for an arbitrary element. In Figure 2(a), we define the Cartesian coordinate system, while in Figure 2(b) , we define the local coordinate system (or natural coordinate system) within a specific domain (i.e. $\xi ,\eta \in \left(-1,1\right)$ ). The relation between these two kinds of coordinate system can be described as

Figure 2.

Isoparametric transition.

$\left\{\begin{array}{c}x\\ y\end{array}\right\}=f\left\{\begin{array}{c}\xi \\ \eta \end{array}\right\}$ E99

which can be further modified by the interpolation function at nodes in the local coordinate system as follows:

$x={\Sigma }_{i=1}^m{N}_i^{\text{'}}{x}_i, y={\displaystyle {\displaystyle \sum }_{i=1}^m}{N^{\prime }}_i{y}_i$ E100

where $\left(x_i,y_i\right)$ are coordinates in the Cartesian coordinate system corresponding to the ith node in local coordinate system, $N_i^{\text{'}}$ is interpolation function of the ith node in local coordinate system and m is the number of nodes chosen to transform the coordinates. Therefore, the regular element in the natural coordinate system can be transformed to the irregular element in the Cartesian coordinate system. The former element is called the parent element, while the latter is called the subelement. Specifically, Eq. (101) can be further expanded as

$\left\{\begin{array}{c}x\\ y\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{N}_1& 0\end{array}& \begin{array}{cc}{N}_2& 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{N}_3& 0\end{array}& \begin{array}{cc}{N}_4& 0\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{cc}0& N_1\end{array}& \begin{array}{cc}0& N_2\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}0& N_3\end{array}& \begin{array}{cc}0& N_4\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_1\\ y_1\end{array}\\ \begin{array}{c}{x}_2\\ y_2\end{array}\end{array}\\ \begin{array}{c}\begin{array}{c}{x}_3\\ y_3\end{array}\\ \begin{array}{c}{x}_4\\ y_4\end{array}\end{array}\end{array}\right\}$ E101

Using the same interpolation functions, the element displacement model can be written as

$\left\{\begin{array}{c}u\\ v\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{N}_1& 0\end{array}& \begin{array}{cc}{N}_2& 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{N}_3& 0\end{array}& \begin{array}{cc}{N}_4& 0\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{cc}0& N_1\end{array}& \begin{array}{cc}0& N_2\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}0& N_3\end{array}& \begin{array}{cc}0& N_4\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_1\\ v_1\end{array}\\ \begin{array}{c}{u}_2\\ v_2\end{array}\end{array}\\ \begin{array}{c}\begin{array}{c}{u}_3\\ v_3\end{array}\\ \begin{array}{c}{u}_4\\ v_4\end{array}\end{array}\end{array}\right\}$ E102

where $J$ denotes the Jacobi matrix while the interpolation functions are given by

$\begin{array}{l}{N}_1=\frac{1}{4}\left(1+\xi \right)\left(1+\eta \right),N_2=\frac{1}{4}\left(1-\xi \right)\left(1+\eta \right)\\ N_3=\frac{1}{4}\left(1+\xi \right)\left(1+\eta \right),N_2=\frac{1}{4}\left(1-\xi \right)\left(1+\eta \right)\end{array}$ E103

As mentioned before, during the derivation of the element stiffness matrix and the equivalent load vector, the derivative of the shape function and the integration in element surface or volume in the Cartesian coordinate system are required. Since the shape functions adopted herein are expressed in natural coordinates, therefore, derivative and integration transformation relationships are essential when isoparametric element is used.

1.13.7. Derivative and integral transformation

According to the law of partial differential,

$\begin{array}{l}\frac{\partial N_i}{\partial \xi }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \xi }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \xi },\\ \frac{\partial N_i}{\partial \eta }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \eta }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \eta },\end{array}$ E104

or in matrix form

$\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\end{array}\right\}=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}=J\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}$ E105

Inverse of Eq. (105) gives

$\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}=J^{-1}\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\end{array}\right\}$ E106

where

$\begin{array}{c}J=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]=\left[\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{i=1}^4}\frac{\partial N_i}{\partial \xi }{x}_i& {\displaystyle {\displaystyle \sum }_{i=1}^4}\frac{\partial N_i}{\partial \xi }{y}_i\\ {\displaystyle {\displaystyle \sum }_{i=1}^4}\frac{\partial N_i}{\partial \eta }{x}_i& {\displaystyle {\displaystyle \sum }_{i=1}^4}\frac{\partial N_i}{\partial \eta }{y}_i\end{array}\right]\\ =\left[\begin{array}{cc}\frac{\partial N_1}{\partial \xi }& \begin{array}{ccc}\frac{\partial N_2}{\partial \xi }& \frac{\partial N_3}{\partial \xi }& \frac{\partial N_4}{\partial \xi }\end{array}\\ \frac{\partial N_1}{\partial \eta }& \begin{array}{ccc}\frac{\partial N_2}{\partial \eta }& \frac{\partial N_3}{\partial \eta }& \frac{\partial N_4}{\partial \eta }\end{array}\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\\ x_4& y_4\end{array}\right]\end{array}$ E107

For an infinitely small element, the area under the Cartesian coordinate system and the natural coordinate system are related by

$ds=dxdy=\left|J\right|d\xi d\eta ,$ E108

where $\left|J\right|$ is the determinant of the Jacobian matrix J . Therefore, element stiffness matrix and equivalent nodal load matrix in Eq. (91) can be transformed to

$\begin{array}{l}{K}^e={\displaystyle {{\displaystyle \int }}_{{\Omega }_e}}{B}^{T}DBt\left|J\right|d\xi d\eta ,P_f^e={\displaystyle {{\displaystyle \int }}_{{\Omega }_e}}{N}^{T}f\left|J\right|d\xi d\eta \\ P_S^e={\displaystyle {{\displaystyle \int }}_{S_{\sigma }{}^e}}{N}^{T}T\left|J\right|d\xi d\eta \end{array}$ E109

For solving the integral equation, usually the Gaussian integration method is employed. In practice, both two and three integration points along each direction of integration are commonly used. Since the discretized system is usually overstiff, it is commonly observed that the use of two integration points along each direction of integration will slightly reduce the stiffness of the matrix and give better results as compared with the use of three integration points. The use of exact integration is possible for some elements, but such approaches are usually tedious and are seldom adopted. The advantage in using the exact integration is that the integration is not affected by the shape of the element while the transformation as shown in Eq. (109) may be affected if the poor shape of the element is poor. The author has developed many finite element programs for teaching and research purposes which can be obtained at [e-mail]. The programs available include plane stress/strain problem, thin/thick plate bending problem, consolidation in 1D and 2D (Biot), seepage problem, slope stability problem, pile foundation problems and others.

1.14. Distinct element method

In practical applications, a limit equilibrium method based on the method of slices or method of columns and strength reduction method based on the finite element method or finite difference method are used for many types of stability problems. These two major analysis methods take the advantage that the in situ stress field which is usually not known with good accuracy is not required in the analysis. The uncertainties associated with the stress-strain relation can also be avoided by a simple concept of factor of safety or the determination of the ultimate limit state. In general, this approach is sufficient for engineering analysis and design. If the condition of the system after failure has initiated is required to be assessed, these two methods will not be applicable. Even if the in situ stress field and the stress-strain relation can be defined, the post-failure collapse is difficult to be assessed using the conventional continuum-based numerical method, as sliding, rotation and collapse of the slope involve very large displacement or even separation without the requirement of continuity.

The most commonly used numerical methods for continuous systems are the FDM, the FEM and the boundary element method (BEM). The basic assumption adopted in these numerical methods is that the materials concerned are continuous throughout the physical processes. This assumption of continuity requires that, at all points in a problem domain, the material cannot be torn open or broken into pieces. All material points originally in the neighbourhood of a certain point in the problem domain remain in the same neighbourhood throughout the whole physical process. Some special algorithms have been developed to deal with material fractures in continuum mechanics-based methods, such as the special joint elements by Goodman [13] and the displacement discontinuity technique in BEM by Crouch and Starfield [5]. However, these methods can only be applied with limitations [21]:

1. large-scale slip and opening of fracture elements are prevented in order to maintain the macroscopic material continuity;

2. the amount of fracture elements must be kept to relatively small so that the global stiffness matrix can be maintained well-posed, without causing severe numerical instabilities; and

3. complete detachment and rotation of elements or groups of elements as a consequence of deformation are either not allowed or treated with special algorithms.

Before a slope starts to collapse, the factor of safety serves as an important index in both the LEM and SRM to assess the stability of the slope. The movement and growth after failure have launched which is also important in many cases that cannot be simulated on the continuum model, and this should be analyzed by the distinct element method (DEM).

In continuum description of soil material, the well-established macro-constitutive equations whose parameters can be measured experimentally are used. On the other hand, a discrete element approach will consider that the material is composed of distinct grains or particles that interact with each other. The commonly used distinct element method is an explicit method based on the finite difference principles which is originated in the early 1970s by a landmark work on the progressive movements of rock masses as 2D rigid block assemblages [6]. Later, the works by Cundall are developed to the early versions of the UDEC and 3DEC codes [9, 10, 12]. The method has also been developed for simulating the mechanical behaviour of granular materials [8], with a typical early code BALL [7], which later evolved into the codes of the PFC group for 2D and 3D problems of particle systems (Itasca, 1995). Through continuous developments and extensive applications over the last three decades, there has accumulated a great body of knowledge and a rich field of literature about the distinct element method. The main trend in the development and application of the method in rock engineering is represented by the history and results of the code groups UDEC/3DEC. Currently, there are many open source (Oval, LIGGGHTS, ESyS, Yade, ppohDEM, Lammps) as well as commercial DEM programs, but in general, this method is still limited to basic research instead of practical application as there are many limitations which include: (1) difficult to define and determine the microparameters; (2) there are still many drawbacks in the use of matching with the macro response to determine the microparameters; (3) not easy to set up a computer model; (4) not easy to include structural element or water pressure; (5) extremely time consuming to perform an analysis; and (6) postprocessing is not easy or trivial. It should also be noted that DEM can be formulated by an energy-based implicit integration scheme which is the discontinuous deformation analysis (DDA) method. This method is similar in many respect to the force-based explicit integration scheme as mentioned previously.

In DEM, the packing of granular material can be defined from statistical distributions of grain size and porosity, and the particles are assigned normal and shear stiffness and friction coefficients in the contact relation. Two types of bonds can be represented either individually or simultaneously; these bonds are referred to the contact and parallel bonds, respectively (Itasca, 1995). Although the individual particles are solid, these particles are only partially connected at the contact points which will change at different time step. Under low normal stresses, the strength of the tangential bonds of most granular materials will be weak and the material may flow like a fluid under very small shear stresses. Therefore, the behaviour of granular material in motion can be studied as a fluid-mechanical phenomenon of particle flow where individual particles may be treated as “molecules” of the flowing granular material. In many particle models for geological materials in practice, the number of particles contained in a typical domain of interest will be very large, similar to the large numbers of molecules.

One of the primary objectives of the particle model is the establishment of the relations between microscopic and macroscopic variables/parameters of the particle systems, mainly through micromechanical constitutive relations at the contacts. Compared with a continuum, particles have an additional degree of freedom of rotation which enables them to transmit couple stresses, besides forces through their translational degrees of freedom. At certain moment, the positions and velocities of the particles can be obtained by translational and rotational movement equations and any special physical phenomenon can be traced back from every single particle interactions. Therefore, it is possible for DEM to analyze large deformation problems and a flow process which will occur after slope failure has initiated. The main limitation of DEM is that there is great difficulty in relating the microscopic and macroscopic variables/parameters; hence, DEM is mainly tailored towards qualitative instead of quantitative analysis.

DEM runs according to a time-difference scheme in which calculation includes the repeated application of the law of motion to each particle, a force-displacement law to each contact, and a contact updating scheme. Generally, there are two types of contact in the program which are the ball-wall contact and the ball-ball contact. In each cycle, the set of contacts is updated from the known particles and known wall positions. Force-displacement law is firstly applied on each contact, and new contact force is then calculated according to the relative motion and constitutive relation. Law of motion is then applied to each particle to update the velocity, the direction of travel based on the resultant force, and the moment and contact acting on the particles. Although every particle is assumed as a rigid material, the behaviour of the contacts is characterized using a soft contact approach in which finite normal stiffness is taken to represent the stiffness which exists at the contact. The soft contact approach allows small overlap between the particles which can be easily observed. Stress on particles is then determined from this overlapping through the particle interface.

1.15. General formulation of DEM

The PFC runs according to a time-difference scheme in which calculation includes the repeated application of the law of motion to each particle, a force-displacement law to each contact, and a contact updating a wall position. Generally, there are two types of contact exist in the program which are ball-to-wall contact and ball-to-ball contact. In each cycle, the set of contacts is updated from the known particle and the known wall position. The force-displacement law is first applied on each contact. New contact force is calculated and replaces the old contact force. The force calculations are based on preset parameters such as normal stiffness, density, and friction. Next, a law of motion is applied to each particle to update its velocity, direction of travel based on the resultant force, moment and contact acting on particle. The force-displacement law is then applied to continue the circulation.

1.16. The force-displacement law

The force-displacement law is described for both the ball-ball and ball-wall contacts. The contact arises from contact occurring at a point. For the ball-ball contact, the normal vector is directed along the line between the ball centres. For the ball-wall contact, the normal vector is directed along the line defining the shortest distance between the ball centre and the wall. The contact force vector F _i is composed of normal and shear component in a single plane surface

$F_i=F_{ij}^n\left(t\right)+F_{ij}^s\left(t+\Delta t\right)$ E110

The force acting on particle i in contact with particle j at time t is given by

$F_{ij}^n\left(t\right)=k_n\left(r_i+r_j-l_{ij}\left(t\right)\right)$ E111

where r _j and r _i stand for particle i and particle j radii, l _ij(t) is the vector joining both centres of the particles and k _n represents the normal stiffness at the contact. The shear force acting on particle i during a contact with particle j is determined by

$F_{ij}^s\left(t+\text{Δ}t\right)=\pm \text{min}(F_{ij}^s\left(t\right)+k_s\text{Δ}{s}_{ij},f|F_{ij}^n\left(t+\text{Δ}t\right)|)$ E112

where f is the particle friction coefficient, k _s represents the tangent shear stiffness at the contact. The new shear contact force is found by summing the old shear force (min F _ij (t)) with the shear elastic force. Δs _ij stands for the shear contact displacement-increment occurring over a time step Δt.

$\Delta s_{ij}=v_{ij}^s\Delta t$ E113

where $V_{ij}^s$ is the shear component of the relative velocity at contact between particles i and j over the time step Δt.

1.17. Law of motion

The motion of the particle is determined by the resultant force and moment acting on it. The motion induced by resultant force is called translational motion. The motion induced by resulting moment is rotational motion. The equations of motion are written in vector form as follows:

1. (Translational motion)

   ${\displaystyle {\displaystyle \sum }_j}{F}_{ij}+m_{i}g+F_i^d=m_i{x}_i^n$ E114

2. (Rotational motion)

   ${\displaystyle {\displaystyle \sum }_j}{r}_i{F}_{ij}+M_i^d=I_r{\theta }_i^n$ E115

where ${x^{″}}_i$ and ${{\theta }^{″}}_i$ stand for the translational acceleration and rotational acceleration of particles i. I _r stands for moment of inertia. $F_i^d$ and $M_i^d$ stand for the damping force and damping moment. Unlike finite element formulation, there are now three degree of freedom for 2D problem and six degree of freedom for 3D problems. In Cundall and Strack’s explicit integration distinct element approach, solution of the global system of equation is avoided by considering the dynamic equilibrium of the individual particles rather than solving the entire system simultaneously. That means, Newton’s law of motion is applied directly. This approach also avoids the generation and storage of the large global stiffness matrix that will appear in finite element analysis. On the other hand, the implicit DDA approach will generate a global stiffness matrix which is even larger than that in finite element analysis, as the rotation is involved directly in the stiffness matrix.

In a typical DEM simulation, if there is no yield by contact separation or frictional sliding, the particles will vibrate constantly and the equilibrium is difficult to be achieved. To avoid this phenomenon which is physically incorrect, numerical or artificial damping is usually adopted in many DEM codes, and the two most common approaches to damping are the mass damping and non-viscous damping. For mass damping, the amount of damping that each particle “feels” is proportional to its mass, and the proportionality constant depends on the eigenvalues of the stiffness matrix. This damping is usually applied equally to all the nodes. As this form of damping introduces body forces, which may not be appropriate in flowing regions, it may influence the mode of failure. Alternatively, Cundall [11] proposed an alternative method where the damping force at each node is proportional to the magnitude of the out-of-balance-force, with a sign to ensure that the vibrational modes are damped rather than the steady motion. This form of damping has the advantage that only accelerating motion is damped and no erroneous damping forces will arise from steady-state motion. The damping constant is also non-dimensional and the damping is frequency independent. As suggested by Itasca [20], an advantage of this approach is that it is similar to the hysteretic damping, as the energy loss per cycle is independent of the rate at which the cycle is executed. While damping is one way to overcome the non-physical nature of the contact constitutive models in DEM simulations, it is quite difficult to select an appropriate and physically meaningful value for the damping. For many DEM simulations, particles are moving around each other and the dominant form of energy dissipation is for frictional sliding and contact breakages. The choice of damping may affect the results of computations. Currently, most of the DEM codes allow the use of automatic damping or manually prescribed the damping if necessary.

To capture the inherent non-linearity behaviour of the problem (with generation and removal of contacts, non-linear contact response and stress-strain behaviour and others), the displacement and contact forces in a given time step must be small enough so that in a single time step, the disturbances cannot propagate from a particle further than its nearest neighbours. For most of the DEM programs, this can be achieved automatically and the default setting is usually good enough for normal cases. It is, however, sometimes necessary to manually adjust the time step in some special cases when the input parameters are unreasonably high or low. Most of the DEM codes use the central difference time integration algorithm which is a second-order scheme in time step.

1.18. Measuring logic

If the local results in DEM are analyzed, it is found that there will be large fluctuations with respect to both locations and time. Such results are not surprising, as the results are highly sensitive to the interaction between particles and hence the time step under which the results are monitored. It can be viewed that such local results can be meaningless unless the results are monitored over a long time span or region. A number of quantities in a DEM model are defined with respect to a specified measurement circle. These quantities include coordinate number, porosity, sliding fraction, stress and strain rate. The coordination number and stress are defined as the average number of contacts per particle. Only particles with centroids that are contained within the measurement circle are considered in computation. In order to account for the additional area of particles that is being neglected, a corrector factor based on the porosity is applied to the computed value of stress.

Since measurement circle is used, stress in particle is described as the two in-plane force acting on each particle per volume of particle. Average stress is defined as the total stress in particle divided by the volume of measurement circle. Thus, shape of particle is regardless of the average stress measurement because the reported stress is easily scaled by volume unity. The reported stress is interpreted as the stress per volume of measurement circle.

1.19. Discussion and conclusion

There are also various publications on the numerical solutions of differential equations, and the readers are suggested to the works of Lee and Schiesser [24], Jovanoic and Suli [22], Veiga et al. [37], Sewell [35], Morton and Mayers [27], Logg et al. [25], Holmes [17], Lui [26], Lapids and Pinder [23] and Iserles [19]. It is impossible for the author to cover every available analytical or numerical method; hence, the author has chosen some methods that are actually used for teaching and research. The readers are strongly encouraged to consult the numerous resources available in various books and publications. There are still new developments available for the solutions of specific differential equations in large-scale problems, and this is also the current trend in the development of differential equation solution.

Due to the importance of the solution of differential equations, there are other important numerical methods that are used by different researchers but are not discussed here, which include the finite difference and boundary element methods (computer codes for learning can also be obtained from the author). Differential equations rely on the Taylor’s series, and the derivatives in the differential equation can be replaced with finite difference approximations on a discretized domain. This will result in a system of algebraic equations that can be solved implicitly or explicitly. There are various ways to form the derivatives, and the most common methods are the forward difference, backward difference and the central difference schemes. While the finite difference methods may be more suitable for different types of differential equations, this method is less convenient to deal with irregular boundary conditions as compared with the finite element method. For highly irregular domain where it is not easy to form a nice discretization, the finite element method will also be much easier and natural to deal with for such condition. In this respect, it is not surprising that many engineering programs are written by the use of the finite element method than the finite difference method.

The boundary element method (BEM) is another numerical method for solving linear partial differential equations which can be formulated as integral equations. The boundary element method uses the given boundary conditions to fit boundary values into the integral equation. In the post-processing stage, the integral equation will be used to calculate the solution directly at any given point inside the solution domain numerically. BEM is applied to problems for which Green’s functions can be calculated, thus this method is initially designed for problems in linear homogeneous media. The dimension of the problem will then be reduced by one. For example, two-dimensional problem will be effectively reduced to one-dimensional problem along the boundary, and this will greatly improve the efficiency of computation. The requirement from the boundary element method imposes considerable restrictions on the range and generality of problems to which the boundary element method can usefully be applied. There are some new developments to the boundary element method so that it can be used for non-linear problem or problems with several major materials (problems with random distribution of material properties are still not applicable). The fundamental solutions are often difficult to integrate. One important property of boundary element analysis is the solution of a fully populated matrix as compared with that in the finite element/difference method. For complicated problems, the boundary element will lose its advantage as compared with other numerical methods. Due to the various limitations, there are only limited boundary element programs available to the researchers. Interested readers can consult the works of Banerjee [1], Brebbia et al. [3] and Trevelyan [36]. It appears that there are less interest in the use and development of the boundary element method in the recent years, due to the various limitations of this method in general non-linear non-homogeneous problem.

In history, various techniques have been developed for ordinary differential equations and partial differential equations under different boundary conditions. While these tricks appear to be elegant, they are not readily adopted for normal engineering use due to various limitations. Being an engineer, the author seldom adopted the methods as outlined in this chapter in actual applications (but do adopt for teaching), except the numerical methods as outlined in this chapter. At present, there are many proprietary or open source finite elements or distinct element codes being used for many complicated real problems. The computer codes (usually in Fortran or C) are usually difficult to be read (if available), and the computer codes for all the partial differential forms (including some extended formats) that have been discussed in this chapter can be readily available from the author for learning purposes. There are also very powerful and general finite element tools or differential equations solver such as FreeFem++, Comsol, Matlab, Mathematica, Maple and Maxima which are used by many scientists and engineers [39, 40]. The use of parallel computing is also strongly influenced by the needs to solve complicated partial differential equations over large solution domain.