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 package’s theory can be divided in two parts:

  • Geometric calculations, like integration over areas

  • Solving Poisson’s Equation, used to compute shear and torsion

The geometric calculations, like computing momentum of inertias, are easy developed. But solving the linear PDE 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). So far, no specific book that uses BEM to compute section properties was found, so we develop some of the equations on this section.

The transformations of equations happens mainly from Green’s Theorem and its identities, which are developed in this text. Since the main sources are from these two books, the repeted citations are ommited.

The theory is divided in parts:

  1. Geometric properties

  2. Torsion Properties

  3. Shear properties

  4. Stress and Strain

  5. Numerical Integration

  6. Boundary Element Method

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 that 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 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. For a ring, the radius of gyration matches its radius

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

Principal Axis Properties

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

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 the intersection of the two neutral lines when only bending momentums are applied.

From construction, it’s the same as the Geometric center \(\mathbf{G}\)

\[\mathbf{B} = \left(x_{bc}, \ y_{bc}\right) := \left(x_{gc}, \ y_{gc}\right) = \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.

\[\nabla^2 \omega = 0\]
\[\left\langle \nabla \omega, \ \mathbf{n}\right\rangle = \mathbf{n} \times \mathbf{p}\]

With \(\mathbf{p} = (x, \ y)\) begin a point on the boundary and \(\mathbf{n}\) the normal vector at \(\mathbf{p}\)

This warping function is found by Boundary Element Method apart from a constant \(c_0\), which is later found in Torsion center.

From now on, we suppose it’s already known.

Torsion constant

The torsion constant can be computed

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

With

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

We transform this integral into a boundary one

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

Torsion center

The torsion center \(\mathbf{T}\) is the point such there’s no shear stresses when a torsion moment is applied.

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

The quantities \(x_{tc}\), \(y_{tc}\) and \(c_0\) can be obtained by minimizing the strain energy produced by axial normal warping stresses, which are ignored by Saint-Venant’s theory. Doing so, leads to 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_0 \\ y_0 \\ -x_0\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

\[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\]

These integrals are transformed to the boundary equivalent.

Boundary reformulation of \(Q_{\omega}\), \(I_{x\omega}\) and \(I_{y\omega}\)

Let \(u\) be a function such

\[\nabla^2 u = g(x, y)\]

Select \(u\) respectivelly as

\[g_{1}(x, \ y) = 1 \Longrightarrow u_{1} = \frac{1}{4}(x^2+y^2)\]
\[g_{x}(x, \ y) = x \Longrightarrow u_{x} = \frac{x^3}{6}\]
\[g_{y}(x, \ y) = y \Longrightarrow u_{y} = \frac{y^3}{6}\]

and use Green’s second identity

\[\begin{split}\int_{\Omega} \omega \cdot g \ d\Omega & = \int_{\Omega} \omega \nabla^2 u - u \nabla^2 \omega \ d\Omega \\ & = \oint_{\Gamma} \omega \dfrac{\partial u}{\partial n} \ d\Gamma - u \dfrac{\partial \omega}{\partial n} \ d\Gamma \\ & = \oint_{\Gamma} \omega \dfrac{\partial u}{\partial n} \ d\Gamma - \oint_{\Gamma} u \cdot \langle \mathbf{p}, \ \mathbf{p}'\rangle \ dt\end{split}\]

Transforming to

\[Q_{\omega} = \dfrac{1}{2}\int_{t_{min}}^{t_{max}} \omega \cdot \mathbf{p} \times \mathbf{p}' \ dt - \dfrac{1}{4}\int_{t_{min}}^{t_{max}} \langle \mathbf{p}, \ \mathbf{p} \rangle \cdot \langle \mathbf{p}, \ \mathbf{p}' \rangle \ dt\]
\[I_{x\omega} = \dfrac{1}{2} \oint_{\Gamma} \omega \cdot x^2 \ dy - \dfrac{1}{6}\int_{t_{min}}^{t_{max}} x^3 \cdot \langle \mathbf{p}, \ \mathbf{p}' \rangle \ dt\]
\[I_{y\omega} = \dfrac{-1}{2} \int_{t_{min}}^{t_{max}} \omega \cdot y^2 \ dx - \dfrac{1}{6}\int_{t_{min}}^{t_{max}} y^3 \cdot \langle \mathbf{p}, \ \mathbf{p}' \rangle \ dt\]

Shear properties

Functions

From Saint-venant theory, the functions \(\Psi\) and \(\Phi\) are fundamental to compute shear properties.

\[\begin{split}\begin{bmatrix} \nabla^2 \Psi \\ \nabla^2 \Phi \end{bmatrix} = 2\begin{bmatrix} -I_{\overline{xx}} & I_{\overline{xy}} \\ I_{\overline{xy}} & -I_{\overline{yy}} \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}\end{split}\]

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}\]

Both equations are in fact Poisson equations. We find them by using the Boundary Element Method apart from constants which are computed in Shear center

Shear center

The shear center \(\boldsymbol{S}\) is the point which

\[\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, it may causes torsion.

TODO

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 know how good the obtained value is. It’s mesured with constants \(n\), \(c\), \(k\) and \(m\), depending on the method

\[\left| \int_{0}^{1} f(x) \ dx - \sum_{i=0}^{n-1} w_i \cdot f(x_i) \right| \le \dfrac{c}{n^{k}} \cdot \max_{x \in \left[0, \ 1\right]} f^{(m)}(x)\]

Polynomial 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_{\Gamma} x^a \cdot y^b \cdot \mathbf{p} \times \mathbf{p}' \ dt\]

Computing it can be done by Regular integrals

Polygonal domains

For polygonal domains, \(I_{a, b}\) can be simplified even more. In that case, each segment is a straight line, so

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

which is constant for an arbitrary segment \(i\). Hence

\[(a+b+2) \cdot I_{a, b} = \sum_{i=0}^{n-1} \left(x_{i}y_{i+1}-x_{i+1}y_{i}\right) I_{a, b, i}\]
\[I_{a, b}^{(i)} = \int_{\Gamma_i} 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}^{(i)} = \sum_{j=0}^{a}\sum_{k=0}^{b}\binom{j+k}{k}\binom{a+b-j-k}{b-k}x_{i}^{a-j}x_{i+1}^{j}y_{i}^{b-k}y_{i+1}^{k}\]

For special cases that \(a=0\) or \(b=0\), we get

\[(a+2)(a+1)I_{a,0} = \sum_{i=0}^{n-1} \left(x_{i}y_{i+1}-x_{i+1}y_{i}\right)\dfrac{x_{i+1}^{a+1}-x_{i}^{a+1}}{x_{i+1}-x_{i}}\]
\[(b+2)(b+1)I_{0,b} = \sum_{i=0}^{n-1} \left(x_{i}y_{i+1}-x_{i+1}y_{i}\right)\dfrac{y_{i+1}^{b+1}-y_{i}^{b+1}}{y_{i+1}-y_{i}}\]

Note

It’s possible to have \(x_{i+1} = x_{i}\) or \(y_{i+1} = y_{i}\) in some segment, which leads to divide by zero in \(I_{a,0}\) and \(I_{0,b}\).

In that case, the expression is opened:

\[\dfrac{z_{i+1}^{c+1}-z_{i}^{c+1}}{z_{i+1}-z_{i}} = \sum_{j=0}^{c} z_{i}^{c-j}z_{i+1}^{j}\]

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 Polygonal domain

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

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}(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(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} \in \partial \Omega\), there is a value \(\tau\) such \(\mathbf{p}(\tau) = \mathbf{s}\) and the angle \(\alpha\) is computed by

\[\mathbf{v}_0 = -\lim_{\delta \to 0^{+}} \mathbf{p}'\left(\tau - \delta\right)\]
\[\mathbf{v}_1 = \lim_{\delta \to 0^{+}} \mathbf{p}'\left(\tau + \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, the first derivative of \(\mathbf{p}\) is continuous and therefore then \(\alpha = \pi\).

Note

In python code, it’s in fact used alpha = arctan2(cross(v0, v1), inner(v0, v1))

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}\]

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}\]

Which vector \(\mathbf{C}\) is a vector of ones.

The determination exacly of the constant depends on the problem and are better treated in Torsion center and Shear center.

Polygonal domain

For polygonal domains, when the basis functions \(\phi(t)\) are piecewise linear, some computations becomes easier. Let’s say the parametric space \(t\) is 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{V}_k = \mathbf{P}_{k+1} - \mathbf{P}_{k}\]
\[\tau = \dfrac{t - t_{k}}{t_{k+1} - t_{k}} \in \left[0, \ 1\right]\]

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

  • If \(t_i \in \left[t_{k}, \ t_{k+1}\right]\) then

    \[\mathbf{r}(t) = \left(\tau-\tau_i\right) \cdot \left(\mathbf{P}_{k+1} - \mathbf{P}_{k}\right)\]
    \[\tau_i = \dfrac{t_i - t_{k}}{t_{k+1} - t_{k}}\in \left[0, \ 1\right]\]

  • Else

    \[\mathbf{r}(t) = \left(\mathbf{P}_{k}-\mathbf{s}_i\right) + \tau \cdot \left(\mathbf{P}_{k+1} - \mathbf{P}_{k}\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 \(P_{k}\) then

\[\mathbf{v}_0 = \mathbf{P}_{k-1}-\mathbf{P}_{k}\]
\[\mathbf{v}_1 = \mathbf{P}_{k+1}-\mathbf{P}_{k}\]
\[\alpha = \arg\left(\langle\mathbf{v}_0, \ \mathbf{v}_1 \rangle + i \cdot \left(\mathbf{v}_0 \times \mathbf{v}_1\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 \(t_i \notin \left[t_k, \ t_{k+1}\right]\), then the evaluation is made by Regular integrals

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

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

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

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\]

With \(\mathbf{P}_k\) begin the \(k\)-vertex and

\[\mathbf{V}_k = \mathbf{P}_{k+1} - \mathbf{P}_k\]
\[\alpha_k = \left\langle \mathbf{P}_k, \ \mathbf{V}_k\right\rangle\]
\[\beta_k = \left\langle \mathbf{V}_k, \ \mathbf{V}_k\right\rangle\]

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

  • If \(t_i \in \left[t_k, \ t_{k+1}\right]\), then
    \[\tau_i = \dfrac{t_i-t_k}{t_{k+1}-t_{k}} \in \left[0, \ 1\right]\]
    \[\mathbf{V}_k = \mathbf{P}_{k+1} - \mathbf{P}_k\]
    \[\mathbf{p(t)} = \mathbf{P}_k + \tau \cdot \mathbf{V}_{k}\]
    \[\mathbf{r(t)} = \left(\tau-\tau_i\right) \cdot \mathbf{V}_{k}\]
    \[\begin{split}F_{ik} = & \int_{0}^{1} \left(\alpha_k + \tau \beta_k \right) \ln\|\left(\tau-\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-\tau_i| dz \\ & + \beta_k \int_{0}^{1} \tau \cdot \ln |\tau-\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 \(\tau_i\). Bellow you find a table with some values

    Table 1 Values of logarithm integrals

    \(\tau_i\)

    \(\int_0^1 \ln|\tau-\tau_i| dz\)

    \(\int_0^1 \tau\ln|\tau-\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 \(t_i\) lies on is made by using analitic values, and singular integrals are not computed.

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

TODO