Engineering at Alberta Courses

Balance Equations: Momentum Balance

Learning Outcomes

Recall Newton’s equation of equilibrium applied to a particle.
Describe how Euler extended the equations of equilibrium to be applied to a continuum.
Identify that in static equilibrium, the stress distribution inside a continuum at equilibrium can be obtained by solving the three equations of equilibrium developed by Euler.
Compute the solution of static equilibrium in simple cases.

Newton and Euler Laws of Motion

Newton’s second law of motions states that the rate of change with respect to the time $t$ of the linear momentum of an object is equal to the net force acting on that object. If $F\in\mathbb{R}^3$ is the force vector acting on an object with mass $m\in\mathbb{R}$ moving with a velocity vector $v\in\mathbb{R}^3$ , then, by denoting the momentum by $p=mv$ , we have:

$F = \frac{dp}{dt} = \frac{dmv}{dt}$

If the mass $m$ is constant then:

$F = m\frac{dv}{dt} = ma$

where $a\in\mathbb{R}^3$ is the acceleration vector of the object.

Euler extended Newton’s second law of motion to a continuum object and introduced two laws. The first law is the balance of linear momentum while the second is the balance of angular momentum of a continuum object. These are as follows:

Euler’s Law of the Balance of Linear Momentum

The first equation is a direct extension of Newton’s first law to an arbitrary volume of a continuum body. Given an arbitrary volume of a continuum body represented by the set $\Omega\in\mathbb{R}^3$ , then the net force acting on the arbitrary volume is equal to the rate of change of the linear momentum of that volume. The net force is the resultant of the traction vectors acting on the boundary and the body forces acting on the continuum points. Thus, according to Euler’s Law:

$\frac{d}{dt}\int_{\Omega} \rho v \hspace{1mm} dx = \int_{\partial\Omega} t_n \hspace{1mm} ds+ \int_{\Omega} \rho b \hspace{1mm} dx$

where, $t$ is time, $\rho:\Omega\times t\rightarrow\mathbb{R}$ is the mass density distribution, $v:\Omega\times t\rightarrow\mathbb{R}^3$ is the velocity vector distribution, $b:\Omega\times t\rightarrow\mathbb{R}^3$ is the body force distribution, $t_n:\partial\Omega\times t\rightarrow\mathbb{R}^3$ is the traction vector distribution acting on the surface of $\Omega$ , and $\mathrm{d}x$ and $\mathrm{d}s$ are the differential volume and surface elements of the arbitrary volume $\Omega$ .

Euler’s Law of the Balance of Angular Momentum

Euler’s law of angular momentum preservation states that the rate of change of the angular momentum of an arbitrary volume of a continuum body around a fixed origin is equal to the torque (measured around the same fixed point) applied by the traction vectors on the boundary and the body forces acting on the continuum points. Thus, according to Euler’s Law:

$\frac{d}{dt}\int_{\Omega} r \times \rho v \hspace{1mm} dx = \int_{\partial\Omega} r \times t_n \hspace{1mm} ds+ \int_{\Omega} r \times \rho b \hspace{1mm} dx$

where, $r:\Omega\times t\rightarrow\mathbb{R}^3$ is the position vector relative to the fixed origin. It should be noted that the fixed origin can be arbitrary. If $r_2$ is the position vector relative to another fixed point such that $r=r_2+r_0$ , then, by replacing $r$ with $r_2+r_0$ and using the fact that $r_0$ is a fixed vector:

$\frac{d}{dt}\int_{\Omega} r_2 \times \rho v \hspace{1mm} dx = \int_{\partial\Omega} r_2 \times t_n \hspace{1mm} ds + \int_{\Omega} r_2 \times \rho b \hspace{1mm} dx + r_0 \times (-\frac{d}{dt}\int_{\Omega} \rho v \hspace{1mm} dx + \int_{\partial\Omega} t_n \hspace{1mm} ds+ \int_{\Omega} \rho b \hspace{1mm} dx$

By using Euler’s law of the balance of linear momentum we reach to the same statement of the balance of angular momentum with $r_2$ replacing $r$ :

$\frac{d}{dt}\int_{\Omega} r_2 \times \rho v \hspace{1mm} dx = \int_{\partial\Omega} r_2 \times t_n \hspace{1mm} ds + \int_{\Omega} r_2 \times \rho b \hspace{1mm} dx$

The Differential Equation of Equation of Equilibrium: Derivation of Using a Differential Volume

Balance of Linear Momentum

Newton’s laws of motion were originally formulated for systems of rigid particles mechanically interacting with each other. Newton’s laws can be extended to a continuum by considering a differential volume and the stresses acting on it. Let $\Omega\subset\mathbb{R}^3$ be a set representing a deformed configuration. Consider a differential cube of mass $\Delta m\in\mathbb{R}$ inside the body with dimensions $\Delta x_1$ , $\Delta x_2$ , $\Delta x_3$ and whose mass density is $\rho=\lim_{\Delta v\rightarrow 0}\frac{\Delta m}{\Delta v}$ . The differential volume of the cube is given by $\Delta v=\Delta x_1 \Delta x_2\Delta x_3$ . Let $F=\{F_1,F_2,F_3\}\in\mathbb{R}^3$ be the resultant force acting on the cube which is due to a body forces vector per unit mass $b=\{b_1,b_2,b_3\}\in\mathbb{R}^3$ acting on the cube and the stresses on the different faces of the cube. To write the differential form of the equilibrium equations, the stresses are assumed to be continuous functions. Such assumption allows the use of the first component in the Taylor expansion for the stresses. Thus, the resultant force component $F_1$ acting on the cube in the direction of the first basis vector can be evaluated as follows:

$F_1 = (-\sigma_{11} + \sigma_{11} + \frac{\partial\sigma_{11}}{\partial x_1}\Delta x_1)\Delta x_2 \Delta x_3 + (-\sigma_{21} + \sigma_{21} + \frac{\partial\sigma_{21}}{\partial x_2}\Delta x_2)\Delta x_1 \Delta x_3$

$+(-\sigma_{31} + \sigma_{31} + \frac{\partial\sigma_{31}}{\partial x_3}\Delta x_3)\Delta x_1 \Delta x_2 + \rho b_1 \Delta x_1 \Delta x_2 \Delta x_3$

$= (\frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{21}}{\partial x_2} + \frac{\partial \sigma_{31}}{\partial x_3} + \rho b_1)\Delta x_1 \Delta x_2 \Delta x_3$

Similarly,

$F_2 = (\frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{32}}{\partial x_3} + \rho b_2)\Delta x_1 \Delta x_2 \Delta x_3$

$F_2 = (\frac{\partial \sigma_{13}}{\partial x_1} + \frac{\partial \sigma_{23}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} + \rho b_3)\Delta x_1 \Delta x_2 \Delta x_3$

If $a=\{a_1,a_2,a_3\}\in\mathbb{R}^3$ is the acceleration vector of the differential cube, then the statement of Newton’s second law dictates that the resultant force is related to the acceleration via:

$F = \Delta m a = \rho \Delta v a = \rho a \Delta x_1 \Delta x_2 \Delta x_3$

Thus, in component form, the equilibrium equations can be written as:

$\frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{21}}{\partial x_2} + \frac{\partial \sigma_{31}}{\partial x_3} + \rho b_1 = \rho a_1$

(1) $\begin{equation*} \frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{32}}{\partial x_3} + \rho b_2 = \rho a_2 \end{equation*}$

$\frac{\partial \sigma_{13}}{\partial x_1} + \frac{\partial \sigma_{23}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} + \rho b_3 = \rho a_3$

The above equations are usually written in the following compact form with $1\leq i\leq 3$ :

$\sum_{j = 1}^{3} \frac{\partial \sigma_{ji}}{\partial x_j} + \rho b_i = \rho a_i$

If the divergencer operator is used, the equations of equilibrium can be written in the following vector form:

$div \sigma + \rho b = \rho a$

Balance of Angular Momentum

The three angular momentum balance equations can be used to derive the symmetry of the stress tensor. Consider a differential rectangular volume oriented with the coordinate system with volume $\Delta v=\Delta x_1 \Delta x_2\Delta x_3$ . Let the stresses on the sides with the negative normals have stresses $\sigma_{ij}$ and the stresses on the sides with positive normals have stresses $\sigma_{ij}+\Delta \sigma_{ij}$ . The first two equations for the balance of linear momentum dictate the following:

$\frac{d}{dt}(\rho v_1 \Delta v) = \Delta \sigma_{11} \Delta x_2 \Delta x_3 + \Delta \sigma_{21} \Delta x_1 \Delta x_3 + \Delta \sigma_{31} \Delta x_1 \Delta x_2 + \rho b_1 \Delta v$

$\frac{d}{dt}(\rho v_2 \Delta v) = \Delta \sigma_{22} \Delta x_1 \Delta x_3 + \Delta \sigma_{12} \Delta x_2 \Delta x_3 + \Delta \sigma_{32} \Delta x_1 \Delta x_2 + \rho b_2 \Delta v$

The moment of the external forces (stresses on the boundary and body forces) acting on the rectangular volume around the axis $e_3$ have the following form:

$M_{external} = r_1 \rho b_2 \Delta v - r_2 \rho b_1 \Delta v + (\sigma_{12} - \sigma_{21})\Delta v$

$+ r_1 \Delta \sigma_{22} \Delta x_1 \Delta x_3 - r_2 \Delta \sigma_{11} \Delta x_2 \Delta x_3$

$+ (r_1 + \Delta /2)\Delta \sigma_{12} \Delta x_2 \Delta x_3 - (r_2 + \Delta x_2/2)\Delta \sigma_{21} \Delta x_1 \Delta x_3$

$+ r_1 \Delta \sigma_{32} \Delta x_1 \Delta x_2 - r_2 \Delta \sigma_{31} \Delta x_1 \Delta x_2$

The rate of change of angular momentum around the axis $e_3$ is denoted $T_3$ and has the form:

$T_3 = r_1 \frac{d \rho v_2 \Delta v}{dt} - r_2 \frac{d \rho v_1 \Delta v}{dt}$

By equating $M_{external}=T_3$ , using the equilibrium equations, and neglecting the smaller terms we get:

$(\sigma_{12} - \sigma_{21}) \Delta v = 0$

Therefore, $\sigma_{12}=\sigma_{21}$ . By repeating the procedure for the moment around other axes we reach:

$\sigma_{ij} = \sigma_{ji}$

The Differential Equation of Equilibrium: Derivation Using Integrals Over Arbitrary Volumes

The same differential equations of motion can be obtained using an equivalent but simpler derivation technique that relies on the integration over an arbitrary volume $\Omega\in\mathbb{R}^3$ .

Balance of Linear Momentum

Using Euler’s law of the balance of linear momentum we have:

$\frac{d}{dt}\int_{\Omega} \rho v \hspace{1mm} dx = \int_{\partial\Omega} t_n \hspace{1mm} ds+ \int_{\Omega} \rho b \hspace{1mm} dx$

The differential operator $\frac{\mathrm{d}}{\mathrm{d}t}$ on the left hand side cannot be interchanged with the integration operator unless the volume of integration is independent of time. Therefore, we can replace the differential volume $\mathrm{d}x$ with $J\mathrm{d}X$ where $J$ is the determinant of the deformation gradient and $\mathrm{d}X$ is the corresponding differential volume in the reference configuration, which is independent of time. We will also utilize the relationship between the traction vector and the stress tensor: $t_n=\sigma^Tn$ , and the divergence theorem which lead to:

$\frac{d}{dt} \int_{\Omega_0} \rho v J \hspace{1mm} dX = \int_{\partial \Omega} \sigma^T n \hspace{1mm} ds + \int_{\Omega} \rho b \hspace{1mm} dx$

$\int_{\Omega_0} \frac{d\rho v J}{dt} \hspace{1mm} dX = \int_{\Omega} div\sigma \hspace{1mm} dx + \int_{\Omega} \rho b \hspace{1mm} dx$

The mass balance equations can be used to simplify the left hand side:

$\int_{\Omega_0} \frac{d \rho J}{dt} v + \frac{dx}{dt} \rho J dX = \int_{\Omega} div \sigma dx + \int_{\Omega} \rho b dx$

$\int_{\Omega_0} \frac{dx}{dt} \rho J dX = \int_{\Omega} div \sigma dx + \int_{\Omega} \rho b dx$

$\int_{\Omega_0} \rho a dx = \int_{\Omega} div \sigma dx + \int_{\Omega} \rho b dx$

$0 = \int_{\Omega} div \sigma + \rho b - \rho a dx$

Since the integration domain is arbitrary, then the integrand is equal to zero. Therefore, the linear momentum balance equation is:

$div\sigma + \rho b = \rho a$

Balance of Angular Momentum

Using Euler’s law of the balance of angular momentum we have:

$\frac{d}{dt}\int_{\Omega} r \times \rho v \hspace{1mm} dx = \int_{\partial\Omega} r \times t_n \hspace{1mm} ds+ \int_{\Omega} r \times \rho b \hspace{1mm} dx$

$\frac{d}{dt} \int_{\Omega_0} r \times \rho v J dX = \int_{\partial \Omega} r \times t_n ds + \int_{\Omega} r \times \rho b dx$

$\int_{\Omega_0} \frac{dr \times \rho v J}{dt}dX = \int_{\partial \Omega} r \times t_n ds + \int_{\Omega} r \times \rho b dx$

$\int_{\Omega_0} r \times \frac{d \rho v J}{dt}dX = \int_{\partial \Omega} r \times t_n ds + \int_{\Omega} r \times \rho b dx$

$\int_{\Omega_0} r \times (\frac{d \rho J}{dt} v + \frac{dv}{dt}\rho J) dX = \int_{\partial \Omega} r \times t_n ds + \int_{\Omega} r \times \rho b dx$

$\int_{\Omega_0} r \times \frac{dv}{dt} \rho J dX = \int_{\partial \Omega} r\times t_n ds + \int_{\Omega} r \times \rho b dx$

Therefore, we have:

$\int_{\Omega} r \times (\rho a - \rho b) dx = \int_{\partial \Omega} r\times t_n ds$

Using the balance of linear momentum equation and replacing $t_n$ with $\sigma^T n$ we get:

$\int_{\Omega} r \times div \sigma dx = \int_{\partial \Omega} r \times \sigma^Tn ds$

The last equation can be written in component form to show the symmetry of the stress matrix. The $i^{th}$ component of the above vector equation can be written as:

$\int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} r_j \frac{\partial \sigma_{lk}}{\partial x_l}dx = \int_{\partial \Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} r_j \sigma_{lk}n_l ds$

Using the divergence theorem , the right hand side can be converted into a volume integral as follows:

$\int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} r_j \frac{\partial \sigma_{lk}}{\partial x_l}dx = \int_{\Omega} \sum_{j,k,l = 1}^3 \frac{\partial \varepsilon_{ijk} r_j \sigma_{lk}}{\partial x_l} ds$

The right hand side can be further simplified as follows:

$\int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} r_j \frac{\partial \sigma_{lk}}{\partial x_l}dx = \int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} (\frac{\partial r_j}{\partial x_l}\sigma_{lk} + r_j\frac{\partial \sigma_{lk}}{\partial x_l}dx$

$\int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} r_j \frac{\partial \sigma_{lk}}{\partial x_l}dx = \int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} (\delta_{jl}\sigma_{lk} + r_j\frac{\partial \sigma_{lk}}{\partial x_l}dx$

$0 = \int_{\Omega} \sum_{j,k,l = 1}^3 \varepsilon_{ijk} \delta_{jl} \sigma_{lk}dx$

$0 = \int_{\Omega} \sum_{j,k = 1}^3 \varepsilon_{ijk} \sigma_{jk}dx$

Therefore, the integrand is equal to zero. Setting $i=1$ , $2$ and $3$ , the three components of the above vector equation have the form:

$0 = \varepsilon_{123}\sigma_{23} + \varepsilon_{132}\sigma_{32}$

$0 = \varepsilon_{213}\sigma_{13} + \varepsilon_{231}\sigma_{31}$

$0 = \varepsilon_{312}\sigma_{12} + \varepsilon_{321}\sigma_{21}$

Therefore:

$\sigma = \sigma^T$

It should be noted that the book by P. Chadwick has another proof showing the same result without relying on the component form of the vector equation.

Solution of the Equilibrium Equations in Static Problems

A continuum problem is solved if a stress field (the stress at every point inside the body) and an acceleration field are found such that the three equations of force equilibrium (conservation of linear momentum equations) are satisfied at every point. On the other hand, the restriction that $\sigma=\sigma^T$ ensures the conservation of angular momentum and no additional information can be obtained from considering the three equations of moment equilibrium. The unknown stress and velocity fields are required to also satisfy certain boundary conditions, which are traditionally given as either displacements or external traction forces on the exterior of the body and initial conditions for velocities. In static equilibrium, which is the purpose of many engineering applications, the acceleration vector is assumed to be zero and, thus, the equilibrium equation is simplified and the stress field becomes the only unknown. In that case, the equations become:

$\frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{21}}{\partial x_2} + \frac{\partial \sigma_{31}}{\partial x_3} + \rho b_1 = 0$

(2) $\begin{equation*} \frac{\partial \sigma_{12}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{32}}{\partial x_3} + \rho b_2 = 0 \end{equation*}$

$\frac{\partial \sigma_{13}}{\partial x_1} + \frac{\partial \sigma_{23}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} + \rho b_3 = 0$

and $\sigma_{ji}=\sigma_{ij}$ . The problem is mathematically formulated as follows:

Let $\Omega_0\subset \mathbb{R}^3$ be a set representing a reference configuration of a continuum body. Given a constant or variable body forces vector $b\in\mathbb{R}^3$ , find the distribution of the stresses inside the continuum body that would satisfy the equations of equilibrium (Eq. 2) above. The boundary conditions for the equations of equilibrium are given on two parts of the boundary of $\Omega_0$ . On the first part, $\Omega_n$ , the external traction vectors $t_n$ are known so we have the boundary conditions for $\sigma$ since $\sigma^Tn=t_n$ . On the second part, $\Omega_u$ , the displacement of the continuum is given. The boundary of $\Omega_0$ is $\partial \Omega_0=\Omega_n\cup \Omega_u$ .

Difficulties Associated with Obtaining a Solution to the Equilibrium Equations

There are two difficulties associated with solving Eq. 2. First, these are three equations of static equilibrium. However, there are six unknowns (six stress variables). Thus, in their current form, many solutions could possibly satisfy the equilibrium equations! The second issue is that some of the boundary conditions contain expressions of displacements, while the equations themselves in this form do not have displacements as variables! These two major issues impel replacing the six unknown stress variables in the above equations with three unknown variables (usually the three displacements $u_1$ , $u_2$ , and $u_3$ . This can be performed by using a “constitutive equation” that describes the relationship between the stress and the strains inside the material. By replacing the stresses with the strains, and by using the relationship between the strains and the displacements described in the strain measures section, the problem becomes well posed, i.e., a solution can be obtained. The following few problems illustrate the applications and the solutions of the equilibrium equations given simplified assumptions allowing such solutions to exist.

The Differential Equation of Equilibrium: Lagrangian Formulation

In sections 6.3.2 and 6.3.3, the Eulerian form of the differential equation of equilibrium was derived. In the Eulerian form, the independent variables are the spatial position coordinates in the deformed configuration $\Omega$ . One difficulty arises from the fact that the set $\Omega$ representing the embedding of the a deforming object is in itself unknown. It is, therefore, desirable to have the differential equations written written in terms of the referential position coordinates represented by the set $\Omega_0$ . In the following two subsections, I will present the derivation of these equations using two approaches, a component form, and a vector form approaches.

Component Form Approach

In the interest of brevity of the derivation, I will be utilizing Einstein Summation Convention. The Eulerian form of the differential equation is given above as:

$\frac{\partial \sigma_{ji}}{\partial x_j} + \rho b_i = \rho a_i$

My intention now is to replace the variables $x_i$ with the variables $X_i$ . For the acceleration and the body forces vectors, that only means a simply change of the argument with the capital letters used for the mappings of the referential coordinates:

$a_i (x,t) = A_i(X,t), \hspace{5mm} b_i(x,t) = B_i(X,t)$

I will recall the density equation that relates the density in the reference configuration $\rho_0$ with the spatial density:

$\rho_0(X,t) = \rho(x,t)J$

I will also recall the equation that relates the First Piola-Kirchhoff stress tensor with the Cauchy stress tensor:

$\sigma_{ji} = P_{ik}F_{jk}\frac{1}{J}$

Substituting the above equations in the Eulerian form of the differential equation yields:

$\frac{\partial (P_{ik}F_{jk}\frac{1}{J})}{\partial x_j} + \rho_0 B_i \frac{1}{J} = \rho_0 A_i \frac{1}{J}$

I will decompose the first term, replace the derivation with respect to the components of $x$ with the derivation with respect to the components of $X$ , and will utilize the fourth statement in the matrix calculus section:

$\frac{\partial (P_{ik}F_{jk}\frac{1}{J})}{\partial x_j} = \frac{\partial P_{ik}}{\partial X_l} \frac{\partial X_l}{\partial x_j} F_{jk} \frac{1}{J} + P_{ik}\frac{\partial F_{jk} \frac{1}{J}}{\partial x_j}$

$= \frac{\partial P_{ik}}{\partial X_l}F_{lj}^{-1}F_{jk} \frac{1}{J} + P_{ik}\frac{1}{J} (\frac{\partial F_jk}{\partial x_j} - \frac{1}{J}F_{jk}\frac{\partial J}{\partial x_j})$

$= \frac{\partial P_{ik}}{\partial X_l}\frac{1}{J} + P_{ik}\frac{1}{J}(\frac{\partial F_{jk}}{\partial x_j} - F_{jk}\frac{\partial F_{il}}{\partial x_j}F_{li}^{-1})$

$= \frac{\partial P_{ik}}{\partial X_l}\frac{1}{J} + P_{ik}\frac{1}{J} (\frac{\partial F_{jk}}{\partial x_j} -\frac{\partial x_j}{\partial X_k}\frac{\partial \frac{\partial x_i}{\partial X_l}}{\partial x_j}\frac{\partial X_l}{\partial x_i})$

$= \frac{\partial P_{ik}}{\partial X_l}\frac{1}{J} + P_{ik}\frac{1}{J} (\frac{\partial F_{jk}}{\partial x_j} - \frac{\partial F_{ik}}{\partial x_i})$

$= \frac{\partial P_{ik}}{\partial X_l}\frac{1}{J}+0$

Therefore, the Lagrangian form of the equilibrium equation can be written in terms of the First Piola-Kirchhoff stress tensor as follows:

$\frac{\partial P_{ij}}{\partial X_j} + \rho_0B_i = \rho_0A_i$

The balance of angular moment imposes the following restriction on the components of $P$ :

$P_{ik}F_{jk} = P_{jk}F_{ik}$

Vector Form Approach

Euler’s law of the balance of linear moment that I showed above is repeated here:

$\frac{d}{dt} \int_{\Omega} \rho v dx = \int_{\partial \Omega} t_n ds + \int_{\Omega} \rho b dx$

By a simple change of variables, $x$ can be replaced by $X$ and the derivative operator $\mathrm{D}$ is used instead of $\mathrm{d}$ to indicate the Lagrangian formulation (material time derivatives with respect to the reference position components):

$\frac{D}{Dt} \int_{\Omega_0} \rho JV dX = \int_{\partial \Omega_0} T_n dS + \int_{\Omega_0} \rho JB dX$

I will use the mass balance equations to simplify the first expression and will use the First Piola-Kirchhoff stress tensor and the reference area vector $N$ for the traction vector: