Introduction

The main sources of theories, hypothesis and equations are both books:

  • Analysis and Design of Elastic Beams, from Walter D. Pilkey

  • The Boundary Element Method for Engineers and Scientists, from John T. Katsikadelis

The used theory can be divided in three parts:

The geometric properties, like computing momentum of inertias, are easy developed. But solving the linear PDE, to use in torsion and shear, is more challenging and more described through this page.

While Pilkey’s book uses Finite Element Method (FEM for short), their equations must be adapted to Boundary Element Method (BEM for short). The transformations of equations are developed throughout in this text. Since the main sources are from these two books, the repeted citations are ommited.

Geometric properties

Cross-section area

The cross sectional are is computed by

\[A := \int_{\Omega} \ dx \ dy\]

First moment of area

The first moment of area are computed by

\[Q_y := \int_{\Omega} x \ dx \ dy\]
\[Q_x := \int_{\Omega} y \ dx \ dy\]

Note

The index \(x\) and \(y\) are switched and they doesn’t represent the internal function

Second Moment of Area

The second moment of inertia are

\[I_{yy} := \int_{\Omega} x^2 \ dx \ dy\]
\[I_{xy} := \int_{\Omega} xy \ dx \ dy\]
\[I_{xx} := \int_{\Omega} y^2 \ dx \ dy\]

They are also known as global second moment of inertia

They can be arranged in a tensor form:

\[\begin{split}\mathbb{I}_{global} := \begin{bmatrix}I_{xx} & I_{xy} \\ I_{xy} & I_{yy}\end{bmatrix}\end{split}\]

Third Moment of Area

The third moment of inertia is computed as:

\[I_{yyy} := \int_{\Omega} x^3 \ dx \ dy\]
\[I_{xyy} := \int_{\Omega} x^2y \ dx \ dy\]
\[I_{xxy} := \int_{\Omega} xy^2 \ dx \ dy\]
\[I_{xxx} := \int_{\Omega} y^3 \ dx \ dy\]

They are used only in Shear center

Geometric center

We denote the geometric centroid by \(\boldsymbol{G}\)

\[\boldsymbol{G} := \left(x_{gc}, \ y_{gc}\right)\]
\[x_{gc} := \dfrac{Q_y}{A}\]
\[y_{gc} := \dfrac{Q_x}{A}\]

Local Second Moment of Area

The local second moment of inertia are computed with respect to the Geometric center \(\boldsymbol{G}\)

\[I_{\overline{yy}} := \int_{\Omega} (x-x_{gc})^2 \ dx \ dy = I_{yy} - \dfrac{Q_{y}^2}{A}\]
\[I_{\overline{xy}} := \int_{\Omega} (x-x_{gc})(y-y_{gc}) \ dx \ dy= I_{xy} - \dfrac{Q_{x}Q_{y}}{A}\]
\[I_{\overline{xx}} := \int_{\Omega} (y-y_{gc})^2 \ dx \ dy= I_{xx} - \dfrac{Q_{y}^2}{A}\]

They can be arranged in a tensor form:

\[\begin{split}\mathbb{I}_{local} := \begin{bmatrix}I_{\overline{xx}} & I_{\overline{xy}} \\ I_{\overline{xy}} & I_{\overline{yy}}\end{bmatrix}\end{split}\]

Radius of Gyration

The radius of gyration is one mesure of spread the body is.

\[r_{x} := \sqrt{\dfrac{I_{xx}}{A}}\]
\[r_{y} := \sqrt{\dfrac{I_{yy}}{A}}\]

For a ring, the radius of gyration matches its radius

Principal Axis Properties

The principals moment of inertia are the eigenvalues of the tensor \(\mathbb{I}_{local}\), from the Local Second Moment of Area.

\[\begin{split}\mathbb{I}_{local} := \begin{bmatrix}I_{\overline{xx}} & I_{\overline{xy}} \\ I_{\overline{xy}} & I_{\overline{yy}}\end{bmatrix}\end{split}\]

For a 2D matrix, \(I_{11}\) and \(I_{22}\) can be easily calculated

\[\Delta = \sqrt{\left(\dfrac{I_{\overline{xx}}-I_{\overline{yy}}}{2}\right)^2+I_{\overline{xy}}^2}\]
\[I_{11} = \dfrac{I_{\overline{xx}}+I_{\overline{yy}}}{2} + \Delta\]
\[I_{22} = \dfrac{I_{\overline{xx}}+I_{\overline{yy}}}{2} - \Delta\]

The direction principal moment of inertia is the eigenvector related to the higher eigenvalue.

It’s also computed as

\[\phi = \arg\left(I_{\overline{xy}} + i \cdot \left(I_{\overline{xx}}-I_{11}\right)\right) = \text{arctan}\left(\dfrac{I_{\overline{xx}}-I_{11}}{I_{\overline{xy}}}\right)\]

Bending Center

The bending center \(\mathbf{B}\) is defined as the point such, when any bending moment is applied, the axial stress is zero.

\[\mathbf{B} := \left(x_{bc}, \ y_{bc}\right)\]

For non-composite sections, it’s the same as the Geometric center \(\mathbf{G}\)

\[\mathbf{B} = \mathbf{G}\]

Torsion Properties

Warping Function

From Saint-venant theory, the warping function \(\omega(x, \ y)\) is fundamental to compute torsion properties.

From \(\omega\), it’s possible to find the Torsion constant, Torsion center and stresses/strains due to Torsion Moment.

\[\begin{split}\begin{align*} \nabla^2 \omega & = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{on} \ \Omega \\ \left\langle \nabla \omega, \ \mathbf{n}\right\rangle & = \mathbf{n} \times \mathbf{p} \ \ \ \ \ \ \text{on} \ \partial \Omega \end{align*}\end{split}\]

The \(\mathbf{n}\) is the normal vector at a point \(\mathbf{p} = (x, \ y)\) on the boundary.

This warping function is found by Boundary Element Method.

Torsion constant

The torsion constant can be computed

\[J := I_{xx} + I_{yy} - \mathbb{J}_{\omega}\]

With \(I_{xx}\) and \(I_{yy}\) found in Second Moment of Area and

\[\mathbb{J}_{\omega} := \int_{\Omega} y \dfrac{\partial \omega}{\partial x} - x \dfrac{\partial \omega}{\partial y} \ dx \ dy\]

Torsion center

The torsion center \(\mathbf{T}\) is defined as the point such, when a torsion moment applied, the shear stresses is zero.

\[\mathbf{T} := \left(x_{tc}, \ y_{tc}\right)\]

The quantities \(x_{tc}\), \(y_{tc}\) and \(c_{tc}\) by solving the linear system.

\[\begin{split}\left(\int_{\Omega} \begin{bmatrix}1 & x & y \\ x & x^2 & xy \\ y & xy & y^2 \end{bmatrix} \ d\Omega\right) \begin{bmatrix}c_{tc} \\ y_{tc} \\ -x_{tc}\end{bmatrix} = \int_{\Omega} \omega\begin{bmatrix}1 \\ x \\ y\end{bmatrix} \ d\Omega\end{split}\]

The matrix on the left side is already computed in

while the values on the right side are computed

\[Q_{\omega} := \int_{\Omega} \omega \ dx \ dy\]
\[I_{x\omega} := \int_{\Omega} x \omega \ dx \ dy\]
\[I_{y\omega} := \int_{\Omega} y \omega \ dx \ dy\]

Note

This torsion center is ignored by Saint-Venant’s theory and comes from elasticity theory by minimizing the strain energy produced by axial normal warping stresses.

Shear properties

Shearing Functions

From Saint-venant theory, the shearing functions \(\Psi\) and \(\Phi\) are fundamental to compute shear properties. They satisfy the Poisson’s equation

\[\nabla^2 \Psi = 2\left(- xI_{\overline{xx}} + y I_{\overline{xy}} \right)\]
\[\nabla^2 \Phi = 2\left(xI_{\overline{xy}} - yI_{\overline{yy}}\right)\]

And boundary conditions

\[\begin{split}\begin{bmatrix}\nabla \Psi \\ \nabla \Phi\end{bmatrix} \cdot \mathbf{n} = \mathbb{H} \cdot \mathbf{n}\end{split}\]
\[\begin{split}\mathbb{H} = \dfrac{\nu}{2}\left((x^2-y^2)\cdot\begin{bmatrix}I_{xx} & I_{xy} \\ -I_{xy} & -I_{yy}\end{bmatrix} + 2xy \cdot \begin{bmatrix}-I_{xy} & I_{xx} \\ I_{yy} & -I_{xy}\end{bmatrix}\right)\end{split}\]

We find them by using the Boundary Element Method.

Shear center

The shear center \(\boldsymbol{S}\) is defined as the point such, when the resultant shear force pass through this point, the generated torsion moment is zero.

\[\boldsymbol{S} := \left(x_{sc}, \ y_{sc}\right)\]
\[\begin{split}\boldsymbol{S} = \dfrac{\nu}{2\Delta}\begin{bmatrix}I_{yy} & I_{xy} \\ I_{xy} & I_{xx}\end{bmatrix}\begin{bmatrix}I_{yyy}+I_{xxy} \\ I_{xyy}+I_{xxx} \end{bmatrix} - \dfrac{1}{\Delta}\int \begin{bmatrix}\Psi \\ \Phi\end{bmatrix} \left\langle \mathbf{p}, \ \mathbf{p}'\right\rangle \ dt\end{split}\]

Which values on the left are the Second Moment of Area and Third Moment of Area and

\[\Delta = 2(1+\nu)(I_{xx}I_{yy}-I_{xy})\]

TODO

Stress and Strain

Introduction

The stress \(\boldsymbol{\sigma}\) and strain \(\boldsymbol{\varepsilon}\) are one of the fundamental quantities to evaluate. They arrive from 4 different phenomenums:

Here we develop expressions to compute stress and strain for any point \(\mathbf{p}\) inside the section. The stress and strain tensor in a beam are given by

\[\begin{split}\boldsymbol{\sigma} = \begin{bmatrix}0 & 0 & \sigma_{xz} \\ 0 & 0 & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz}\end{bmatrix} \ \ \ \ \ \ \ \ \ \boldsymbol{\varepsilon} = \begin{bmatrix}\varepsilon_{xx} & 0 & \varepsilon_{xz} \\ 0 & \varepsilon_{yy} & \varepsilon_{yz} \\ \varepsilon_{xz} & \varepsilon_{yz} & \varepsilon_{zz} \end{bmatrix}\end{split}\]

The elasticity law relates both tensors

\[\boldsymbol{\sigma} = \lambda \cdot \text{trace}\left(\boldsymbol{\varepsilon}\right) \cdot \mathbf{I} + 2\mu \cdot \boldsymbol{\varepsilon}\]
\[\boldsymbol{\varepsilon} = \dfrac{1+\nu}{E} \cdot \boldsymbol{\sigma} - \dfrac{\nu}{E} \cdot \text{trace}\left(\boldsymbol{\sigma}\right) \cdot \mathbf{I}\]

With \(\lambda\) and \(\mu\) being Lamé Parameters, \(E\) beging Young Modulus and \(\nu\) the Poisson’s coefficient.

\[\lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)} \ \ \ \ \ \ \ \ \ \ \ \mu = \dfrac{E}{2(1+\nu)}\]
\[E = \dfrac{\mu\left(3\lambda+2\mu\right)}{\lambda+\mu} \ \ \ \ \ \ \ \ \ \ \ \nu = \dfrac{\lambda}{2(\lambda+\mu)}\]

To clear the equations, sometimes we use the pair \(\left(\lambda, \ \mu\right)\), other times we use \(\left(E, \ \nu\right)\)

Axial Force

The axial force only leads to axial stress. Meaning, in pure axial force case, the stress tensor and strain are given by

\[\begin{split}\boldsymbol{\varepsilon} = \begin{bmatrix}\varepsilon_{xx} & 0 & 0 \\ 0 & \varepsilon_{yy} & 0 \\ 0 & 0 & \varepsilon_{zz}\end{bmatrix} \ \ \ \ \ \ \ \ \ \ \ \sigma = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \sigma_{zz}\end{bmatrix}\end{split}\]

The axial stress is constant when an axial force \(\mathrm{F}_{z}\) is given by

\[\sigma_{zz} = \dfrac{\mathrm{F}_{z}}{A}\]

Where \(A\) is the Cross-section area.

Hence, the strain is given by elasticity law:

\[\varepsilon_{xx} = \varepsilon_{yy} = -\nu \cdot \dfrac{\mathrm{F}_{z}}{EA}\]
\[\varepsilon_{zz} = \dfrac{\mathrm{F}_{z}}{EA}\]
\[\begin{split}\boldsymbol{\varepsilon} = \dfrac{\mathrm{F}_{z}}{EA}\begin{bmatrix}-\nu & 0 & 0 \\ 0 & -\nu & 0 \\ 0 & 0 & 1\end{bmatrix}\end{split}\]

Bending Moments

The bending moments \(\mathrm{M}_{x}\) and \(\mathrm{M}_{y}\) causes only axial stresses. The tensors are

\[\begin{split}\boldsymbol{\varepsilon} = \begin{bmatrix}\varepsilon_{xx} & 0 & 0 \\ 0 & \varepsilon_{yy} & 0 \\ 0 & 0 & \varepsilon_{zz}\end{bmatrix} \ \ \ \ \ \ \ \ \ \ \ \sigma = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \sigma_{zz}\end{bmatrix}\end{split}\]

The expression of \(\sigma_{zz}\) depends on the position of the point \(\mathbf{p}\) in the section. In the Bending Center \(\boldsymbol{B}\) the stress and the strain are zero while they increase/decrease depending on the distance to the bending center.

Let \(\bar{x}=x-x_{bc}\) and \(\bar{y}=y-y_{bc}\), the function \(\sigma_{zz}(x, \ y)\) satisfy

\[\int_{\Omega} \sigma_{zz} \ d\Omega = 0\]
\[\begin{split}\int_{\Omega} \sigma_{zz} \cdot \begin{bmatrix}\bar{y} \\ -\bar{x}\end{bmatrix} \ d\Omega = \begin{bmatrix}M_{x} \\ M_{y}\end{bmatrix}\end{split}\]

Add the hypothesis that \(\sigma_{zz}\) is linear with respect to \(x\) and \(y\), then

\[\begin{split}\sigma_{zz}(x, \ y) & = \dfrac{1}{\det \left(\mathbb{I}_{local}\right)} \begin{bmatrix}\bar{y} & \bar{x}\end{bmatrix} \left[\mathbb{I}_{b}\right] \begin{bmatrix}M_{x} \\ M_{y}\end{bmatrix} \\ & = -\left(\dfrac{I_{\overline{xy}}\mathrm{M}_{x} + I_{\overline{xx}}\mathrm{M}_{y}}{I_{\overline{xx}}I_{\overline{yy}}-I_{\overline{xy}}^2}\right) \cdot \bar{x} + \left(\dfrac{I_{\overline{yy}}\mathrm{M}_{x} + I_{\overline{xy}}\mathrm{M}_{y}}{I_{\overline{xx}}I_{\overline{yy}}-I_{\overline{xy}}^2}\right) \cdot \bar{y}\end{split}\]

With constants given in Local Second Moment of Area

The neutral line is the set of pairs \((x, \ y)\) such \(\sigma_{zz}(x, \ y) = 0\). That means the neutral line is the line that pass thought \(\boldsymbol{B}\) and it’s parallel to the vector \(\left[\mathbb{I}_{b}\right] \cdot \left(\mathrm{M}_{x}, \ \mathrm{M}_{y}\right)\)

It’s possible to obtain strain values from elasticity law:

\[\varepsilon_{xx} = \varepsilon_{yy} = -\nu \cdot \dfrac{\sigma_{zz}}{E}\]
\[\varepsilon_{zz} = \dfrac{\sigma_{zz}}{E}\]
\[\begin{split}\boldsymbol{\varepsilon} = \dfrac{\sigma_{zz}}{E} \cdot \begin{bmatrix}-\nu & 0 & 0 \\ 0 & -\nu & 0 \\ 0 & 0 & 1\end{bmatrix}\end{split}\]

Torsion Moment

The torsion moment \(\mathrm{M}_{z}\) causes only shear stresses. The tensors are

\[\begin{split}\boldsymbol{\varepsilon} = \begin{bmatrix}0 & 0 & \varepsilon_{xz} \\ 0 & 0 & \varepsilon_{yz} \\ \varepsilon_{xz} & \varepsilon_{yz} & 0\end{bmatrix} \ \ \ \ \ \ \ \ \ \ \ \boldsymbol{\sigma} = \begin{bmatrix}0 & 0 & \sigma_{xz} \\ 0 & 0 & \sigma_{yz} \\ \sigma_{xz} & \sigma_{xz} & 0\end{bmatrix}\end{split}\]

The Warping Function \(\omega\) is used to compute them

\[\sigma_{xz}(x, \ y) = \dfrac{\mathrm{M}_{z}}{J} \cdot \left(\dfrac{\partial \omega}{\partial x} - y\right)\]
\[\sigma_{yz}(x, \ y) = \dfrac{\mathrm{M}_{z}}{J} \cdot \left(\dfrac{\partial \omega}{\partial y} + x\right)\]
\[\varepsilon_{xz}(x, \ y) = \dfrac{1}{2\mu} \cdot \sigma_{xz}\]
\[\varepsilon_{yz}(x, \ y) = \dfrac{1}{2\mu} \cdot \sigma_{yz}\]

Which \(J\) is the Torsion constant and \(\mu\) is the second Lamé Parameter.

To compute the partial derivatives, two approaches are used:

  • For a point \(\mathbf{p}\) on the boundary

    \[\begin{split}\nabla \omega & = \dfrac{\partial \omega}{\partial t} \cdot \mathbf{t} + \dfrac{\partial \omega}{\partial n} \cdot \mathbf{n} \\ & = \left\langle \mathbf{p}, \ \mathbf{t}\right\rangle \cdot \mathbf{n} + \mathbf{t} \cdot \sum_{j=0}^{n-1} \varphi_{j}'(t) \cdot W_{j}\end{split}\]

    The derivatives by themselves don’t matter, but the evaluation of \(\sigma_{xz}\) and \(\sigma_{yz}\), which are rewritten as

    \[\begin{split}\begin{bmatrix}\sigma_{xz} \\ \sigma_{yz}\end{bmatrix} = \dfrac{\mathrm{M}_z}{J} \cdot \left[\left\langle\mathbf{p}, \ \mathbf{n}\right\rangle + \sum_{j=0}^{n-1}\varphi_{j}'(t) \cdot W_{j}\right] \cdot \mathbf{t}\end{split}\]

  • For interior points, \(\mathbf{p} \in \text{open}\left(\Omega\right)\)

Shear Forces

The shear forces \(\mathrm{F}_{x}\) and \(\mathrm{F}_{y}\) causes only shear stresses. The tensors are

\[\begin{split}\boldsymbol{\varepsilon} = \begin{bmatrix}0 & 0 & \varepsilon_{xz} \\ 0 & 0 & \varepsilon_{yz} \\ \varepsilon_{xz} & \varepsilon_{yz} & 0\end{bmatrix} \ \ \ \ \ \ \ \ \ \ \ \boldsymbol{\sigma} = \begin{bmatrix}0 & 0 & \sigma_{xz} \\ 0 & 0 & \sigma_{yz} \\ \sigma_{xz} & \sigma_{xz} & 0\end{bmatrix}\end{split}\]

Depending on the application of the shear force, as seen at Shear center, it may cause a torsion moment.

For computations of these shear stresses, we suppose both shear forces pass through the shear center and therefore no torsion is caused.

TODO

Scalar transformations

To evaluate if a material will fail under a certain load condition, criterias are used to transform transform the tensor of stress \(\mathbb{\sigma}\) into a scalar value \(\sigma_{eq}\).

\[\begin{split}\mathbb{\sigma} = \begin{bmatrix}0 & 0 & \sigma_{xz} \\ 0 & 0 & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz} \end{bmatrix} \Rightarrow \sigma_{eq}\end{split}\]

A material will not fail if the criteria is satisfied, for the yield strength \(S_{y}\) obtained from a tensile test.

\[\sigma_{eq} < S_{y}\]

The most used criterias are Von Mises and Tresca:

  • Von Mises

\[\sigma_{eq} = \sqrt{\sigma_{zz}^2 + 3\left(\sigma_{xz}^2 + \sigma_{yz}^2\right)}\]

  • Tresca

\[\sigma_{eq} = \sqrt{\dfrac{\sigma_{zz}^2}{4} + \sigma_{xz}^2 + \sigma_{yz}^2}\]

Boundary Element Method

Introduction

The Boundary Element Method (BEM for short) is a method that solves a linear PDE by transforming the problem in a boundary problem. Once the problem is solved, all the informations on the boundary are known and then the interior informations are easy computed after that.

In our case, BEM is used to solve the laplace’s equation

(1)\[\nabla^2 u = 0\]

BEM transforms Eq. (1) into a boundary version Eq. (2)

(2)\[\alpha\left(\mathbf{s}\right) \cdot u\left(\mathbf{s}\right) = \int_{\Gamma} u \cdot \dfrac{\partial v}{\partial n} \ d\Gamma - \int_{\Gamma} \dfrac{\partial u}{\partial n} \cdot v \ d\Gamma\]

Which \(\mathbf{s}\) is the source point of the Green function \(v\) and \(\alpha(\mathbf{s})\) is the angle at the point \(\mathbf{s}\).

(3)\[v(\mathbf{p}, \ \mathbf{s}) = \ln r = \ln \|\mathbf{r}\| = \ln \|\mathbf{p} - \mathbf{s}\|\]

Since all the PDEs used in this package have only Neumann’s boundary conditions, all values of \(\dfrac{\partial u}{\partial n}\) are known and the objective is finding all the values of \(u\) at the boundary.

Once \(u\) and \(\dfrac{\partial u}{\partial n}\) are known at the boundary, it’s possible to compute \(u(x, y)\) and its derivatives at any point inside by using Eq. (2).

Solution at the boundary

Parametrize the curve \(\Gamma\) by \(\mathbf{p}(t)\)

(4)\[\mathbf{p}(t) = \sum_{j=0}^{m-1} \phi_{j}(t) \cdot \mathbf{P}_{j} = \langle \mathbf{\phi}(t), \ \mathbf{P}\rangle\]

Set \(u(t)\) as a linear combination of \(n\) basis functions \(\varphi(t)\) and weights \(\mathbf{U}\).

(5)\[u(t) = \sum_{j=0}^{n-1} \varphi_j(t) \cdot U_j = \langle \mathbf{\varphi}(t), \ \mathbf{U}\rangle\]

Fix the source point \(\mathbf{s}_i = \mathbf{p}(\hat{t}_i)\) at the boundary and expand Eq. (2) by using Eq. (5) to get Eq. (6)

(6)\[\sum_{j=0}^{n-1} A_{ij} \cdot U_{j} = \sum_{j=0}^{n-1} M_{ij} \cdot U_{j} - F_{i}\]

With the auxiliar values which depends only on the geometry, the source point and the basis functions

\[A_{ij} = \alpha\left(\mathbf{s}_i\right) \cdot \varphi_j\left(\hat{t}_i\right)\]
\[M_{ij} = \int_{\Gamma} \varphi_j \cdot \dfrac{\partial v_i}{\partial n} \ d\Gamma\]
\[F_{i} = \int_{\Gamma} \dfrac{\partial u}{\partial n} \cdot v_i \ d\Gamma\]

Applying for \(n\) different source points \(\mathbf{s}_i\) at boundary, we get the matrices \(\mathbb{A}\), \(\mathbb{M}\) and \(\mathbf{F}\) such

(7)\[\left(\mathbb{M}-\mathbb{A}\right) \cdot \mathbf{U} = \mathbf{F}\]

Finding the values of \(\mathbf{U}\) means solving the linear system Eq. (7). The objective then is computing these matrices to solve Eq. (7).

Matrix \(\mathbb{A}\)

The angle \(\alpha\) is the mesure for a given point with respect to its position to the domain \(\Omega\).

\[\begin{split}\alpha\left(\mathbf{s}\right) = \begin{cases}\in \left(0, \ 2\pi\right) \ \ \ \ \text{if} \ \mathbf{s} \in \partial \Omega \\ 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{if} \ \mathbf{s} \notin \text{closed}\left(\Omega\right) \\ 2\pi \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{if} \ \mathbf{s} \in \text{open}\left(\Omega\right) \end{cases}\end{split}\]

When \(\mathbf{s}_i \in \partial \Omega\), there is a value \(\hat{t}_i\) such \(\mathbf{p}(\hat{t}_i) = \mathbf{s}_i\) and the angle \(\alpha\) is computed by

\[\mathbf{v}_0 = -\lim_{\delta \to 0^{+}} \mathbf{p}'\left(\hat{t}_i - \delta\right)\]
\[\mathbf{v}_1 = \lim_{\delta \to 0^{+}} \mathbf{p}'\left(\hat{t}_i + \delta\right)\]
\[\alpha = \arg\left(\langle\mathbf{v_0}, \ \mathbf{v_1} \rangle + i \cdot \left(\mathbf{v_0} \times \mathbf{v_1}\right)\right)\]

For smooth regions, \(\mathbf{p}'\) is continuous and therefore then \(\alpha = \pi\).

Matrix \(\mathbb{M}\)

We use

\[\dfrac{\partial v}{\partial n} ds = \dfrac{\mathbf{r} \times \mathbf{p}'}{\left\langle\mathbf{r}, \ \mathbf{r}\right\rangle} \ dt\]

to write

\[M_{ij} = \int_{t_{min}}^{t_{max}} \varphi_{j}(t) \cdot \dfrac{\mathbf{r} \times \mathbf{p}'}{\left\langle\mathbf{r}, \ \mathbf{r}\right\rangle} \ dt\]

Vector \(\mathbf{F}\) for warping

For the warping function

\[\dfrac{\partial u}{\partial n} = \mathbf{n} \times \mathbf{p} = \dfrac{\langle \mathbf{p}, \ \mathbf{p}'\rangle}{\|\mathbf{p}'\|}\]
\[F_i = \int_{t_{min}}^{t_{max}} \left\langle \mathbf{p}, \ \mathbf{p}'\right\rangle \cdot \ln \|\mathbf{r}_i\| \ dt\]

Vector \(\mathbf{F}\) for shear

TODO

Evaluating matrices

The matrices highly depend on the geometry and the basis functions \(\varphi\).

To compute the coefficients \(M_{ij}\) and \(F_{i}\), it’s used numerical integration, like Gaussian-Quadrature. Unfortunatelly, when \(r = 0\) at some point, the integrants are singular and special techniques are used.

The main idea to compute them is decompose the integral in intervals and use

For polygonal domains the Inside integration is not required cause it can be done analiticaly. But for higher degrees, it’s indeed necessary

Constraint solution

Although the matrix \(\mathbb{K}=\mathbb{M}-\mathbb{A}\) is not singular, all the PDEs have Neumann’s boundary conditions and has no unique solution. If \(u^{\star}\) is found as solution, then \(\left(u^{\star} + \text{const}\right)\) also is a solution.

Although both functions give the same properties cause it envolves only the derivatives of \(u\), we restrict the solution by solving the system with Lagrange Multiplier.

\[\begin{split}\begin{bmatrix}K & \mathbf{C} \\ \mathbf{C}^T & 0\end{bmatrix} \begin{bmatrix}\mathbf{U} \\ \lambda \end{bmatrix} = \begin{bmatrix}\mathbf{F} \\ 0\end{bmatrix}\end{split}\]

The matrix \(\mathbf{C}\) is a matrix obtained from the Torsion center and Shear center. It comes from minimizing the total strain energy.

Lagrange Multiplier is also used to restrict the solution when there are two boundary curves that touch each other and the expected solution must be equal for all points that both curves touch each other as described in Multiple curves.

Numerical Integration

Regular integrals

The numerical integral are computated by using quadrature schemas, rewriting

\[\int_{0}^{1} f(x) \ dx = \sum_{i=0}^{n-1} w_i \cdot f(x_i)\]

With specific position nodes \(x_i\) and weights \(w_i\). \(n\) is the number of integration points

Depending of the nodes and weights, we get different approximations. Although the error is unknown, it’s still possible to estimate a value called residual.

\[\left| \int_{0}^{1} f(x) \ dx - \sum_{i=0}^{n-1} w_i \cdot f(x_i) \right| \le R(n, f)\]

There are nodes and weights for different quadratures.

Closed Newton Cotes Quadrature

\(n\)

\(x_i\)

\(w_i\)

2

0

1/2

1

1/2

3

0

1/6

0.5

4/6

1

1/6
Open Newton cotes Quadrature

\(n\)

\(x_i\)

\(w_i\)

1

1/2

1

2

0

1/2

1

1/2

3

1/4

2/3

2/4

-1/3

3/4

2/3
Gaussian Quadrature

\(n\)

\(x_i\)

\(w_i\)

1

1/2

1
Chebyshev Quadrature

\(n\)

\(x_i\)

\(w_i\)

1

1/2

1

Singular integrals

Singular integrals are used when the integrating function is not defined in the entire interval due to singularities. We decompose the integrating function in two functions:

  • The weight function \(g(x)\), such contains known singularities

  • The integrable function \(f(x)\), which is a unknown function defined in all interval

Therefore, we compute

\[\int_{0}^{1} f(x) \cdot g(x) \ dx \approx \sum_{i=0}^{n-1} w_i \cdot f(x_i)\]

With \(n\) specific position nodes \(x_i\) and weights \(w_i\), computed depending on the fonction \(g(x)\) and the position of the singularities.

For our specific case, there are only two types of singular integrals developed in Boundary Element Method:

\[\int_{0}^{1} f(x) \cdot \ln x \ dx\]
\[\int_{-1}^{1} f(x) \cdot \dfrac{1}{x} \ dx\]

Note

The current implementation allows only polygonal domains. Hence, singular integrals are evaluated analiticaly as shown in Boundary Element Method

Logarithm singularity

We are interested in computing the integral

\[I = \int_{0}^{1} f(x) \ \cdot \ln x \ dx\]

Describing the function \(f(x)\) by taylor series

\[f(x) = \sum_{i=0}^{\infty} a_i \cdot x^{i}\]

The integral is well defined

\[I = - \sum_{i=0}^{\infty} \dfrac{a_i}{\left(1+i\right)^2}\]

Although it’s well defined, in general the \(a\) coefficients are unknown.

A logarithm quadrature was created by Stroud and Sladek with given values in table bellow

\[\int_{0}^{1} f(x)\ln x \ dx \approx -\sum_{i=0}^{n-1} w_{i} \cdot f(x_{i})\]
Nodes and Weights for Logarithm Quadrature

\(n\)

\(x_i\)

\(w_i\)

2

0.112008806166976

0.718539319030384

0.602276908118738

0.281460680969615

3

0.0638907930873254

0.513404552232363

0.368997063715618

0.391980041201487

0.766880303938941

0.0946154065661491

Odd singularity

We are interested in computing the integral

\[\int_{-1}^{1} \dfrac{1}{x} \cdot f(x) \ dx\]

The given integral is computed as the Cauchy Principal Value, which symbol is further ommited

\[PV\int_{-1}^{1} \dfrac{f(x)}{x} \ dx = \lim_{\varepsilon \to 0^{+}} \int_{-1}^{-\varepsilon} \dfrac{f(x)}{x} \ dx + \int_{\varepsilon}^{1} \dfrac{f(x)}{x} \ dx\]

This integral is well defined:

\[\int_{-1}^{1} \dfrac{1}{x} \ dx = 0\]
\[\int_{-1}^{1} \dfrac{x}{x} \ dx = 2\]
\[\int_{-1}^{1} \dfrac{x^2}{x} \ dx = 0\]
\[\int_{-1}^{1} \dfrac{1}{x} \cdot f(x) \ dx = \sum_{j=0}^{\infty} \dfrac{2 \cdot a_{2j+1}}{2j+1}\]

It’s possible to create a quadrature for it:

TODO

Bidimensional integrals

To compute area, momentums and inertias, it’s needed to compute the integral

\[I_{a,b} = \int_{\Omega} x^a \cdot y^b \ dx \ dy\]

Which \(\Omega\) is the defined region with closed boundary \(\Gamma\), \(a\) and \(b\) are natural numbers

By using Green’s thereom, we transform the integral

\[\int_{\Omega} \left(\dfrac{\partial Q}{\partial x} - \dfrac{\partial P}{\partial y}\right) \ dx \ dy = \int_{\Gamma} P \ dx + Q \ dy\]

Without loss of generality, let \(\alpha \in \mathbb{R}\) and take

\[\dfrac{\partial Q}{\partial x} = \alpha \cdot x^a \cdot y^b \Longrightarrow Q = \dfrac{\alpha}{a+1} \cdot x^{a+1} \cdot y^b\]
\[\dfrac{\partial P}{\partial y} = \left(\alpha-1\right) \cdot x^a \cdot y^b \Longrightarrow P = \dfrac{\alpha - 1}{b+1} \cdot x^{a} \cdot y^{b+1}\]

Then

\[I_{a, b} = \dfrac{\alpha - 1}{b+1} \int_{\Gamma} x^{a} \cdot y^{b+1} \ dx + \dfrac{\alpha}{a+1} \int_{\Gamma} x^{a+1} \cdot y^b \ dy\]

This integral is computed in the boundary and the expression depends on \(\alpha\).

In special, by taking \(\alpha = \dfrac{a+1}{a+b+2}\), it’s transformed to

\[(a+b+2) \cdot I_{a, b} = \int_{t_{min}}^{t_{max}} x^a \cdot y^b \cdot \mathbf{p} \times \mathbf{p}' \ dt\]

Computing it can be done by Regular integrals

Polygonal domains

Adding the hypothesis that boundary curves are polygonal, many integrals shown before can be simplified, which are developed bellow.

Let’s say first the polygonal has \(m\) sides, then the parametric space \(t\) can be divided by the knots \(t_0\), \(t_1\), \(\cdots\), \(t_{m-1}\), \(t_m\), which correspond to the vertices.

For an arbitrary interval \(\left[t_k, \ t_{k+1}\right]\), \(\mathbf{p}(t)\) is described as

\[\mathbf{p}(t) = \mathbf{P}_{k} + \tau \cdot \mathbf{V}_k\]
\[\mathbf{p}'(t) = \mathbf{V}_k\]
\[\mathbf{V}_k = \mathbf{P}_{k+1} - \mathbf{P}_{k}\]
\[\tau = \dfrac{t - t_{k}}{t_{k+1} - t_{k}} \in \left[0, \ 1\right]\]

For later use, define

\[\alpha_k = \left\langle \mathbf{P}_k, \ \mathbf{V}_k\right\rangle\]
\[\beta_k = \left\langle \mathbf{V}_k, \ \mathbf{V}_k\right\rangle\]

Bidimensional integrals

It’s our objective computing

\[I_{a, b} = \int_{\Omega} x^a \cdot y^b \ dx \ dy\]

Continuating from Bidimensional integrals, it was transformed to

\[\begin{split}\begin{align*} (a+b+2) \cdot I_{a, b} & = \int_{t_{min}}^{t_{max}} x^a \cdot y^b \cdot \mathbf{p} \times \mathbf{p}' \ dt \\ & = \sum_{k=0}^{m-1} \int_{t_{k}}^{t_{k+1}} x^a \cdot y^b \cdot \mathbf{p} \times \mathbf{p}' \ dt \end{align*}\end{split}\]

For the segment \(\left[t_{k}, \ t_{k+1}\right]\).

\[\mathbf{p}(t) \times \mathbf{p}'(t) = \mathbf{P}_{k} \times \mathbf{P}_{k+1}\]

Hence

\[(a+b+2) \cdot I_{a, b} = \sum_{k=0}^{m-1} \left(x_{k}y_{k+1}-x_{k+1}y_{k}\right) I_{a, b}^{(k)}\]
\[I_{a, b}^{(k)} = \int_{0}^{1} x^a \cdot y^b \ dt\]

The integral can be computed by expanding it and using the beta function:

\[\int_{0}^{1} (1-t)^a \cdot t^b \ dt = \dfrac{1}{a+b+1} \cdot \dfrac{1}{\binom{a+b}{a}}\]

Leading to

\[(a+b+1)\binom{a+b}{a} I_{a, b}^{(k)} = \sum_{i=0}^{a}\sum_{j=0}^{b}\binom{i+j}{j}\binom{a+b-i-j}{b-j}x_{k}^{a-i}x_{k+1}^{i}y_{k}^{b-j}y_{k+1}^{j}\]

Boundary Element Method

For polygonal domains, when the basis functions \(\phi(t)\) are piecewise linear, some computations becomes easier.

Since the source point \(\mathbf{s}_i = \mathbf{p}(\hat{t}_i)\),

  • If \(\hat{t}_i \notin \left[t_{k}, \ t_{k+1}\right]\) then

    \[\mathbf{r}_i(t) = \left(\mathbf{P}_{k}-\mathbf{s}_i\right) + \tau \cdot \mathbf{V}_{k}\]
  • If \(\hat{t}_i \in \left[t_{k}, \ t_{k+1}\right]\)
    \[\mathbf{r}(t) = \left(\tau-\hat{\tau}_i\right) \cdot \mathbf{V}_{k}\]
    \[\hat{\tau}_i = \dfrac{\hat{t}_i - t_{k}}{t_{k+1} - t_{k}}\in \left[0, \ 1\right]\]

Matrix \(\mathbb{A}\)

If the source point \(\mathbf{s}_i\) lies in the middle of the segment

\[\alpha(\mathbf{s}_i) = \pi\]

If the source point \(s_i\) lies in the vertex \(\mathbf{P}_{k}\) then

\[\alpha = \arg\left(\langle\mathbf{V}_{k-1}, \ \mathbf{V}_{k} \rangle + i \cdot \left(\mathbf{V}_{k-1} \times \mathbf{V}_{k}\right)\right)\]

Matrix \(\mathbb{M}\)

\[M_{ij} = \sum_{k=0}^{m-1} \int_{t_{k}}^{t_{k+1}} \varphi_{j} \cdot \dfrac{\mathbf{r} \times \mathbf{p}'}{\left\langle \mathbf{r}, \mathbf{r}\right\rangle} \ dt\]

  • If \(\hat{t}_i \notin \left[t_k, \ t_{k+1}\right]\), then the evaluation is made by Regular integrals

  • If \(\hat{t}_i \in \left[t_k, \ t_{k+1}\right]\)

    \[\mathbf{p}(t) = \mathbf{P}_k + \tau \cdot \mathbf{V}_{k}\]
    \[\mathbf{r}(t) = \left(\tau-\hat{\tau}_i\right) \cdot \mathbf{V}_{k}\]
    \[\mathbf{r} \times \mathbf{p}' = 0\]

    Therefore, the integration over the interval \(\left[t_k, \ t_{k+1}\right]\) can be skipped

Vector \(\mathbf{F}\) for warping

For warping function, the expression \(F_i\) is written as

\[\dfrac{\partial u}{\partial n} = \dfrac{\left\langle \mathbf{p}, \ \mathbf{p}'\right\rangle}{\|\mathbf{p}'\|}\]
\[F_{i} = \sum_{k=0}^{m-1} \int_{0}^{1} \left(\alpha_k + \tau \cdot \beta_k \right) \ln\|\mathbf{r}\| \ d\tau\]

  • If \(\hat{t}_i \notin \left[t_k, \ t_{k+1}\right]\), Regular integrals are used

  • If \(\hat{t}_i \in \left[t_k, \ t_{k+1}\right]\), then
    \[\mathbf{r}(t) = \left(\tau-\hat{\tau}_i\right) \cdot \mathbf{V}_{k}\]
    \[\begin{split}F_{ik} = & \int_{0}^{1} \left(\alpha_k + \tau \beta_k \right) \ln\|\left(\tau-\hat{\tau}_i\right) \cdot \mathbf{V}_k\| \ d\tau \\ = & \left(\alpha_{k} + \dfrac{1}{2}\beta_{k}\right) \cdot \dfrac{1}{2}\ln \beta_k \\ & + \alpha_{k} \int_{0}^{1} \ln |\tau-\hat{\tau}_i| dz \\ & + \beta_k \int_{0}^{1} \tau \cdot \ln |\tau-\hat{\tau}_i| \ dz\end{split}\]

    These two log integrals are computed analiticaly, the expressions are complicated (here and here) and depends on the value of \(\hat{\tau}_i\). Bellow you find a table with some values

    Table 1 Values of logarithm integrals

    \(\hat{\tau}_i\)

    \(\int_0^1 \ln|\tau-\hat{\tau}_i| dz\)

    \(\int_0^1 \tau\ln|\tau-\hat{\tau}_i| dz\)

    \(0\)

    \(-1\)

    \(\frac{-1}{4}\)

    \(\frac{1}{2}\)

    \(-(1+\ln 2)\)

    \(\frac{-1}{2}\left(1+\ln 2\right)\)

    \(1\)

    \(-1\)

    \(\frac{-3}{4}\)

    Therefore, the integral over interval which \(\hat{t}_i\) lies on is made by using analitic values, and singular integrals are not computed.

Vector \(\mathbf{F}\) for shear

TODO

Domain to boundary

In this section we develop the expressions that transform the domain integrals over \(\Omega\) to the boundary integrals over \(\Gamma\)

One example is used in Bidimensional integrals

Torsion constant

The expression used in Torsion constant is

\[\mathbb{J}_{\omega} := \int_{\Omega} y \dfrac{\partial \omega}{\partial x} - x \dfrac{\partial \omega}{\partial y} \ dx \ dy\]

Which is manipulated to

\[\mathbb{J}_{\omega} = \int_{\Omega} \nabla \omega \times \mathbf{p} \ dx \ dy\]

Therefore

\[\mathbb{J}_{\omega} = \int_{t_{min}}^{t_{max}} \omega \cdot \left\langle \mathbf{p}, \ \mathbf{p}'\right\rangle \ dt\]

Torsion center

The used quantities to evaluate the torsion center are

\[Q_{\omega} := \int_{\Omega} \omega \ dx \ dy\]
\[I_{x\omega} := \int_{\Omega} x \omega \ dx \ dy\]
\[I_{y\omega} := \int_{\Omega} y \omega \ dx \ dy\]

To develop the expression, we develop arbitrary expression for \(g\):

\[I_{g\omega} := \int_{\Omega} g \cdot \omega \ dx \ dy\]

Say there’s a function \(q\) such

\[\nabla^2 q = g\]

Since \(\nabla^2 \omega = 0\) and using Green’s identity

\[\begin{split}\begin{align*} I_{g\omega} & = \int_{\Omega} g \cdot \omega \ dx \ dy \\ & = \int_{\Omega} \omega \cdot \nabla^2 q \ dx \ dy \\ & = \int_{\Omega} \omega \cdot \nabla^2 q - q \cdot \nabla^2 \omega \ dx \ dy \\ & = \int_{\Gamma} \omega \cdot \dfrac{\partial q}{\partial n} - q \cdot \dfrac{\partial \omega}{\partial n} \ d\Gamma \\ & = \int_{t_{min}}^{t_{max}} \omega \cdot \left(\nabla q \times \mathbf{p}'\right) dt - \int_{t_{min}}^{t_{max}} q \cdot \langle \mathbf{p}, \ \mathbf{p}'\rangle \ dt \end{align*}\end{split}\]

By selecting \(q\) we can get the expressions:

\[\begin{split}\begin{align*} q = \dfrac{1}{4}\left(x^2+y^2\right) & \Rightarrow g = 1 \\ q = \dfrac{x}{12}\left(x^2+3y^2\right) & \Rightarrow g = x \\ q = \dfrac{y}{12}\left(3x^2+y^2\right) & \Rightarrow g = y \end{align*}\end{split}\]

Multiple curves

So far the developed expressions considered only one curve.

Although it’s enough for most cases, it’s possible to have hollowed domains or with multiple materials

TODO