%path = "physics/mechanics/forces" %kind = kinda["examples"] %level = 10 #in school years
This pin jointed truss is the initial example take from akabaila
It can be analysed using force and moment vectors, because resulting linear equations are neither underdetermined nor overdetermined, but determined. The truss is said to be statically determinate.
We want to find the forces along the members.
To find all the forces along the members we use:
no revolution: all moments must be matched by reacting moments
no translation: all forces in all nodes add to zero.
We will solve this here using Python
with numpy
and scipy
and more precisely we will do with these functions:
1 2 3 4 5 6 | from numpy import dot, cross, array, allclose, transpose, finfo, double;
from scipy.linalg import norm, solve;
unit = lambda v: v/norm(v)
R90 = array([[0, 1], [-1, 0]]) # rotate 90 degrees
V = array
eps = finfo(double).eps
|
First we need to describe the problem world, i.e. we need to find the variables. By variable I mean the real thing here. A variable consists of values. When using a value of a variable, then we do this via reference, more precisely via index into the list of values representing the variable.
The following describes the system. I use captial letters for variables and small letters for references to values of the variables. The names are short, one letter if possible. I take the first letter of the english word. N is Nodes, E is Edges, F is Forces, S is Support nodes. n, e, f and s reference values of these variables. A value of Nodes consists of x and y, which reference the external variables X and Y (the values of). By external I mean that they are not specified in the code.
1 2 3 4 5 6 7 8 9 | #nodes = x, y
N = [V((a*8,b*6)) for a,b in [(0,0),(0,1),(1,0),(1,1),(2,0),(2,1),(3,0),(3,1)]];
#edge = (n1,n2), n_i = indix into N
E = [(0, 1), (1, 3), (1, 2), (0, 2), (2, 3), (3, 5), (3, 4),
(2, 4), (4, 5), (5,7), (4, 6), (4, 7), (6, 7)]
#external forces = index into N for node of application, vector
F = [(2,V((0,-100))), (4,V((0,-150)))]
#support points = indices into N
S = [0,6];
|
Now let’s find the forces along the edges.
No revolution.
We need to make the moment created around one support point zero by constructing a force at the other support point. If there were more than one other support point for an axis, the system would be overdetermined, which we don’t handle here.
1 2 3 4 5 6 7 8 9 10 11 | def react_to_mp_at_q(mp,q):
"""p != q are any nodes.
m stands for moment. mp is the moment around node p.
"""
dp=N[q]-N[mp[0]]
norm_dp = norm(dp)
if norm_dp == 0:
return V((0,0))
ndp = dp/norm_dp
fq = mp[1]*dot(R90,ndp)/norm_dp
return -fq
|
No translation
We distribute the forces to a node to those edges not having a force associated yet. In our 2D case we need two such edges. One is OK, if the force is exactly in that direction.
The force placed on an edge via this distribution will be forwarded to the other node, but there the direction must be changed: An edge under tension will pull from both nodes and a contracted edge will push into both nodes.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | def distribute(f,es,q):#f = sum of forces on edges es to node q
ies = [i for i in range(len(E)) if q in E[i]]
mat = []
eo = []#edge, other node
for e in ies:
if e not in es:
#E[e]
t = [tt for tt in E[e] if tt==q][0]#this
o = [tt for tt in E[e] if tt!=q][0]#other
d0 = unit(N[o]-N[t])
mat.append(d0)
eo.append((e,o))
A = transpose(array(mat))
dim = len(f)
if len(eo)==dim:
r = solve(A,f)
for i in range(len(r)):
ff = r[i]*mat[i]
yield ff, eo[i]#even if ff==0
elif len(eo) > dim:
for i,v in enumerate(mat):
angle = dot(v,f)/norm(v)/norm(f)
if abs(angle) < eps or abs(angle+1) < 4 * eps: #same direction
yield f,eo[i]
return
raise ValueError('node %i overdetermined'%q)
else:
if allclose(unit(f),mat[0]):
yield f, eo[0]
else:
raise ValueError('node %i underdetermined'%q)
|
The above distribute
needs the edges along which forces come into the node.
We keep track of the edges with forces in a {node, [(edge,force)..]}
dictionary.
Initially this is empty. We add the external forces and the forces due to the moments.
Then we distribute forces in unbalanced nodes.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | def no_revolution():
EF = dict([(p,[]) for p in range(len(N))])
for p,ff in F:
EF[p].append(([],ff))
for i in range(len(S)):
for j in range(len(S)):
if j != i:
p = S[i]
q = S[j]
mp = (p,sum([cross(ff[1],(N[ff[0]]-N[p])) for ff in F]))
fq = react_to_mp_at_q(mp,q)
EF[q].append(([],fq))
return EF
def no_translation(EF):
_sum = lambda tt: [reduce(lambda x,y:x+y,t) for t in zip(*tt)]
unbalanced = lambda:[(i,v) for i,v in [(i,_sum(EF[i])) for i in EF]
if v and not allclose(norm(v[1]),0)]
u = unbalanced()
while len(u)>0:
q,(es,f) = u[0]
dist=list(distribute(f,es,q))
for ff,eo in dist:
EF[q].append(([eo[0]],-ff)) #q is this node of edge eo[0]
EF[eo[1]].append(([eo[0]],ff)) #eo[1] is the other node of edge eo[0]
u = unbalanced()
return EF
def format_ef(EF):
from itertools import chain
from pprint import pformat
e_f = list(chain.from_iterable([[(tuple(e),norm(f)) for e,f in EF[i]] for i in EF]))
e_f = list([(e[0],round(f,2)) for e,f in set(e_f) if e])
e_f.sort()
return pformat(e_f)
EF = no_revolution()
EF = no_translation(EF)
format_ef(EF)
|
Here are the resulting forces along the edges:
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 |
---|---|---|---|---|---|---|---|---|---|---|---|---|
116.67 |
155.56 |
194.44 |
0.0 |
16.67 |
177.78 |
27.78 |
155.56 |
0.0 |
177.78 |
0.0 |
222.22 |
133.33 |