And now it’s all this

I just said what I said and it was wrong
Or was taken wrong

Next post Previous post

Green’s theorem and section properties

January 17, 2018 at 11:07 AM by Dr. Drang

In the last post, I presented a simple Python module with functions for calculating section properties of polygons. Now we’ll go through the derivations of the formulas used in those functions.

The basis for all the formulas is Green’s theorem, which is usually presented something like this:

$$\underset{C}{\oint }Pdx+Qdy=\underset{A}{\iint }(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y})dxdy$$

where $P$ and $Q$ are functions of $x$ and $y$, $A$ is the region over which the right integral is being evaluated, and $C$ is the boundary of that region. The integral on the right is evaluated in accordance with the right-hand rule, i.e., counterclockwise for the usual orientation of the $x$ and $y$ axes.

The section properties of interest are all area integrals. We’ll use Green’s theorem to turn them into boundary integrals and then evaluate those integrals using the coordinates of the polygon’s vertices.

Area

This is the easiest one, but instead of going through the full derivation here, I’ll refer you to this excellent StackExchange page by apnorton and just hit the highlights.

1. The area is defined

   $$A=\underset{A}{\iint }dxdy$$

   and we’ll choose $P=0$ and $Q=x$ as our Greens’ theorem functions. This gives us

   $$A=\underset{A}{\iint }dxdy=\underset{C}{\oint }xdy$$

2. We break the polygonal boundary into a series of straight-line segments, each of which can be parameterized this way:

   $$x=x_i+(x_{i+1}-x_i)t$$ $$y=y_i+(y_{i+1}-y_i)t$$ $$dy=(y_{i+1}-y_1)dt$$

   where the $(x_i,y_i)$ are the coordinates of the vertices.

3. Plugging these equations into the integral, we get

   $$A=\frac{1}{2}{\sum }_{i=0}^{n-1}(x_{i+1}+x_i)(y_{i+1}-y_i)$$

A note on the indexing: The polygon has $n$ vertices, which we’ll number from 0 to $n-1$. The last segment of the boundary goes from $(x_{n-1},y_{n-1})$ to $(x_0,y_0)$. To make this work with the equation, we’ll define $(x_n,y_n)=(x_0,y_0)$.

Let’s compare this with the area function in the module:

python
def area(pts):
  'Area of cross-section.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  s = 0
  for i in range(len(pts) - 1):
    s += x[i]*y[i+1] - x[i+1]*y[i]
  return s/2

We start by checking the pts list to see if the starting and ending items match. If they don’t, we copy the starting item to the end to fit the indexing convention discussed above. We then initialize some variables and execute a loop, summing terms along the way. Rewriting the loop in mathematical terms, we get

$$A=\frac{1}{2}{\sum }_{i=0}^{n-1}{x}_i{y}_{i+1}-x_{i+1}{y}_i$$

This doesn’t look like the equation derived from Green’s theorem, does it? But it’s not too hard to see that they are equivalent. Expanding out the binomial product in the earlier equation gives

$$x_{i+1}{y}_{i+1}-x_{i+1}{y}_i+x_i{y}_{i+1}-x_i{y}_i$$

As we loop through all the values of $i$ from 0 to $n-1$, the leading term of one trip through the loop will cancel the trailing term of the next trip through the loop. Here’s an example for a triangle:

$$(x_1{y}_1-x_1{y}_0+x_0{y}_1-x_0{y}_0)$$ $$+(x_2{y}_2-x_2{y}_1+x_1{y}_2-x_1{y}_1)$$ $$+(x_0{y}_0-x_0{y}_2+x_2{y}_0-x_2{y}_2)$$

After the cancellations, all that’s left are the inner terms, and that’s the formula used in the area function.

The cancellation doesn’t do much for us here, changing from two additions and one multiplication per loop to two multiplications and one addition per loop. But we’ll see this same sort of cancellation in the other section properties, and it will provide greater simplification in those.

Centroid

The centroid is essentially the average position of the area. If a sheet of material of uniform thickness and density were cut into a shape, the centroid would be the center of gravity, the balance point, of that shape. The coordinates of the centroid are defined this way:

$$x_c=\frac{1}{A}\underset{A}{\iint }xdxdy$$ $$y_c=\frac{1}{A}\underset{A}{\iint }ydxdy$$

Let’s derive the formula for $x_c$ for a polygon; the derivation of the formula for $y_c$ will be similar.

In applying Green’s theorem, we’ll take $P=0$ and $Q=\frac{1}{2}{x}^2$. Therefore,

$$x_c=\frac{1}{A}\underset{A}{\iint }xdxdy=\underset{C}{\oint }\frac{1}{2}{x}^{2}dy$$

Breaking the polygonal boundary into straight-line segments and using the same parametric equations as before, we get an integral that looks like this

$${\int }_0^1\frac{1}{2}[(x_i^2+2x_i(x_{i+1}-x_i)t+(x_{i+1}-x_i{)}^2{t}^2](y_{i+1}-y_i)dt$$

for each segment. This integral evaluates to

$$\frac{1}{6}(x_{i+1}^2+x_{i+1}{x}_i+x_i^2)(y_{i+1}-y_i)$$

so our formula for the centroid is

$$x_c=\frac{1}{6A}{\sum }_{i=0}^{n-1}(x_{i+1}^2+x_{i+1}{x}_i+x_i^2)(y_{i+1}-y_i)$$

As we found in the formula for area, the leading and trailing terms in the expansion of this product cancel out as we loop through the sum, leaving us with

$$x_c=\frac{1}{6A}{\sum }_{i=0}^{n-1}-y_i{x}_{i+1}^2+y_{i+1}{x}_{i+1}{x}_i-y_i{x}_{i+1}{x}_i+y_{i+1}{x}_i^2$$

This looks like a mess, but it can be factored into a more compact form:

$$x_c=\frac{1}{6A}{\sum }_{i=0}^{n-1}(x_{i+1}+x_i)(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

The expression for the other centroidal coordinate is as you’d expect:

$$y_c=\frac{1}{6A}{\sum }_{i=0}^{n-1}(y_{i+1}+y_i)(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

These are the formulas used in the centroid function.

python:
def centroid(pts):
  'Location of centroid.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  sx = sy = 0
  a = area(pts)
  for i in range(len(pts) - 1):
    sx += (x[i] + x[i+1])*(x[i]*y[i+1] - x[i+1]*y[i])
    sy += (y[i] + y[i+1])*(x[i]*y[i+1] - x[i+1]*y[i])
  return sx/(6*a), sy/(6*a)

Moments and product of inertia

You may be familiar with moments and products of inertia from dynamics, where the terms are related to the distribution of mass in a body. The moments and product of inertia we’ll be talking about here—more properly called the second moments of area—are mathematically similar and refer to the distribution of area across a planar shape.

The moments and product of inertia that matter in beam bending are taken about the centroidal axis (i.e., a set of $x$ and $y$ axes with the origin at the centroid of the shape). Since we don’t know where the centroid is when we set up our coordinate system, our list of vertex points don’t work off that basis. But we can still calculate the centroidal moments and product of inertia by using these formulas:

$$I_{xx}=\underset{A}{\iint }(y-y_c{)}^{2}dxdy$$ $$I_{yy}=\underset{A}{\iint }(x-x_c{)}^{2}dxdy$$ $$I_{xy}=\underset{A}{\iint }(y-y_c)(x-x_c)dxdy$$

We’ll concentrate on $I_{yy}$; the other two will be similarly derived.

First, let’s expand the square inside the integral and see what we get:

$$I_{yy}=\underset{A}{\iint }{x}^{2}dxdy-2x_c\underset{A}{\iint }xdxdy+x_c^2\underset{A}{\iint }dxdy$$

The integral in the second term is $A{x}_c$ and the integral in the third term is just $A$. Putting this together, we get¹

$$I_{yy}=\underset{A}{\iint }{x}^{2}dxdy-A{x}_c^2$$

Since we already have formulas for $x$ and $x_c$, we can concentrate on the integral in the first term on the right.

Returning to Green’s theorem, we’ll use $P=0$ and $Q=\frac{1}{3}{x}^3$, giving us

$$\underset{A}{\iint }{x}^{2}dxdy=\underset{C}{\oint }\frac{1}{3}{x}^{3}dy$$

Once again, we break the polygonal boundary into straight-line segments and use parametric equations to define the segments. For each segment, we’ll get the following integral:

$${\int }_0^1\frac{1}{3}[x_i^3+3x_i^2(x_{i+1}-x_1)t+3x_i(x_{i+1}-x_1{)}^2{t}^2+(x_{i+1}-x_1{)}^3{t}^3](y_{i+1}-y_i)dt$$

This integral evaluates to

$$\frac{1}{12}[x_{i+1}^3+x_{i+1}^2{x}_i+x_{i+1}{x}_i^2+x_i^3](y_{i+1}-y_i)$$

giving us

$$\underset{A}{\iint }{x}^{2}dxdy=\frac{1}{12}{\sum }_{i=0}^{n-1}[x_{i+1}^3+x_{i+1}^2{x}_i+x_{i+1}{x}_i^2+x_i^3](y_{i+1}-y_i)$$

Once again, if we expand out the product inside the sum, we’ll find that the leading and trailing terms cancel as we work through the loop. That gives us

$$\frac{1}{12}{\sum }_{i=0}^{n-1}-y_i{x}_{i+1}^3+y_{i+1}{x}_{i+1}^2{x}_i-y_i{x}_{i+1}^2{x}_i+y_{i+1}{x}_{i+1}{x}_i^2-y_i{x}_{i+1}{x}_i^2+y_{i+1}{x}_i^3$$

And that long expression can be factored, leaving

$$\underset{A}{\iint }{x}^{2}dxdy=\frac{1}{12}{\sum }_{i=0}^{n-1}(x_{i+1}^2+x_{i+1}{x}_i+x_i^2)(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

Similar² derivations give us

$$\underset{A}{\iint }{y}^{2}dxdy=\frac{1}{12}{\sum }_{i=0}^{n-1}(y_{i+1}^2+y_{i+1}{y}_i+y_i^2)(x_i{y}_{i+1}-x_{i+1}{y}_i)$$ $$\underset{A}{\iint }xydxdy=\frac{1}{24}{\sum }_{i=0}^{n-1}(2x_{i+1}{y}_{i+1}+x_{i+1}{y}_i+x_i{y}_{i+1}+2x_i{y}_i)(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

These formulas, and the terms accounting for the location of the centroid, are in the function inertia.

python:
def inertia(pts):
  'Moments and product of inertia about centroid.'

  if pts[0] != pts[-1]:
    pts = pts + pts[:1]
  x = [ c[0] for c in pts ]
  y = [ c[1] for c in pts ]
  sxx = syy = sxy = 0
  a = area(pts)
  cx, cy = centroid(pts)
  for i in range(len(pts) - 1):
    sxx += (y[i]**2 + y[i]*y[i+1] + y[i+1]**2)*(x[i]*y[i+1] - x[i+1]*y[i])
    syy += (x[i]**2 + x[i]*x[i+1] + x[i+1]**2)*(x[i]*y[i+1] - x[i+1]*y[i])
    sxy += (x[i]*y[i+1] + 2*x[i]*y[i] + 2*x[i+1]*y[i+1] + x[i+1]*y[i])*(x[i]*y[i+1] - x[i+1]*y[i])
  return sxx/12 - a*cy**2, syy/12 - a*cx**2, sxy/24 - a*cx*cy

This older post explains the use of the moment of inertia in beam bending, but I avoided the trickier bits associated with the product of inertia and principal axes. We’ll cover them in the next post.

Next post Previous post

Site search

Meta

• drdrang at leancrew
• Blog archive
• RSS feed
• JSON feed
• Mastodon
• GitHub repositories

Recent posts

Credits

This work is licensed under a Creative Commons Attribution-Share Alike 3.0 Unported License.

© 2005–2023, Dr. Drang