Material and structural properties

Start by building the simple planar structure from The structure of Structures:

import numpy as np
nodes = np.array([[0,0,0],[0,1,0],[1,0,0],[1,1,0]]).transpose()
members = np.array([[0,1],[1,2],[2,3],[3,0],[0,2],[1,3]]).transpose()

from tnsgrt import Structure
s = Structure(nodes, members, number_of_strings=4)

which looks like

as generated by the following code:

from matplotlib import pyplot as plt
from tnsgrt.plotter.matplotlib import MatplotlibPlotter
plotter = MatplotlibPlotter()
plotter.plot(s)
fig, ax = plotter.get_handles()
ax.view_init(90,-90)
ax.axis('off')
ax.axis('equal')
plt.show()

Equilibrium

One can perform various calculations on structures. For example one can determine the member forces so that the structure is in equilibrium in various scenarios.

A key method is tnsgrt.structure.Structure.equilibrium(), which calculates the internal forces required to maintain a structure in equilibrium. It is assumed that all nodes act as ball joints.

Unloaded

For example:

reactions = s.equilibrium()

calculates the forces in the members that maintain the structure in equilibrium and return any reactions.

In this unloaded case, no external forces are applied to the structure, and equilibrium is achieved by pretensioning the structure. Because the structure is free of constraints, the reactions

reactions

are all zeros:

array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])
s.member_properties[['lambda_', 'force']]

Running tnsgrt.structure.Structure.equilibrium() updates the structure's properties force and force coefficients, which equals the forces divided by member length. The result of the equilibrium calculation can be found in the updated member properties lambda_ (force coefficient) and force:

s.member_properties[['lambda_', 'force']]

which, in this example, returns:

	lambda_	force
0	1.0	1.000000
1	1.0	1.000000
2	1.0	1.000000
3	1.0	1.000000
4	-1.0	-1.414214
5	-1.0	-1.414214

Pretension is set so that the average force coefficient on all bars is equal to the parameter lambda_bar, which is by default equal to one.

Loaded

In this case an equilibrium is sought in the presence of external forces, given as a 3 x n array as the following one:

f = 0.125*np.array([[0,1,0],[0,-1,0],[0,-2,0],[0,2,0]]).transpose()

Each column is a force vector to be applied on the corresponding node.

The external force array f can then be passed on to the method equilibrium:

reactions = s.equilibrium(f)

resulting in the new set of member forces and force coefficients:

s.member_properties[['lambda_', 'force']]

that returns:

	lambda_	force
0	1.250000e-01	1.250000e-01
1	2.500000e-01	2.500000e-01
2	-3.758633e-11	-3.758633e-11
3	2.500000e-01	2.500000e-01
4	-2.500000e-01	-3.535534e-01
5	-2.500000e-01	-3.535534e-01

The following code produces a visualization of the applied forces superimposed on the structure:

plotter = MatplotlibPlotter()
plotter.plot(s)
fig, ax = plotter.get_handles()
plotter.plot_arrows(s.nodes, f, color='g')
ax.view_init(90,-90)
ax.axis('off')
plt.show()

resulting in a figure like

The forces are represented by the green arrows.

When it is not possible to find a set of internal forces that satisfy the equilibrium conditions an Exception with a message "could not find equilibrium" is produced. For example:

f = 0.125*np.array([[0,1,0],[0,-1,0],[0,-1,0],[0,2,0]]).transpose()
s.equilibrium(f)

can not be in equilibrium only by internal forces.

Stiffness

Once a structure is in equilibrium, its response to forces can be calculated in terms of its stiffness matrix. For that it is necessary to characterize the members' geometry and material properties. The fundamental properties are the member radius, and elasticity modulus:

s.member_properties[['radius', 'inner_radius', 'modulus']]

The current default values for such properties are:

	radius	modulus
0	0.005	2.000000e+11
1	0.005	2.000000e+11
2	0.005	2.000000e+11
3	0.005	2.000000e+11
4	0.010	2.000000e+11
5	0.010	2.000000e+11

For calculating the stiffness matrix of a pretensioned structure, it is also necessary to know the member’s force coefficient and the derived member stiffness property. As seen before, the force coefficient and the force are obtained during the equilibrium calculation:

s.equilibrium()
s.member_properties[['lambda_', 'force', 'stiffness']]

which returns:

	lambda_	force
0	1.0	1.000000
1	1.0	1.000000
2	1.0	1.000000
3	1.0	1.000000
4	-1.0	-1.414214
5	-1.0	-1.414214

Because the stiffness is a “derived” property, it does not get automatically populated, which can be done by calling tnsgrt.structure.Structure.update_member_properties():

s.update_member_properties('stiffness')
s.member_properties[['stiffness']]

to obtain:

	stiffness
0	1.570796e+07
1	1.570796e+07
2	1.570796e+07
3	1.570796e+07
4	4.442883e+07
5	4.442883e+07

After setting the material properties, one can calculate the stiffness model associated with the current equilibrium:

stiffness, _, _ = s.stiffness()

For large models, the stiffness is stored and calculated as sparse arrays. However, for small models, such as this one, the model is stored in dense arrays. The warning message can be suppressed by explicitly setting the parameter storage=dense:

stiffness, _, _ = s.stiffness(storage='dense')

WARNING: setting storage='dense' for large models is not advised.

The stiffness model can be used to calculate various quantities of interest. One example is the calculation of the approximate displacements generated in response to a set of nodal forces. For the nodal forces:

f = 0.125*np.array([[0,1,0],[0,-1,0],[0,-2,0],[0,2,0]]).transpose()

these approximate displacements can be, in principle, calculated using:

stiffness.displacements(f)

which in this case results in an error.

The failure of the above procedure is due to the singularity of the current stiffness model. This can be visualized by calculating the model's eigenvalues and eigenvectors:

d, v = stiffness.eigs()

In this case, because there are not enough constraints in the possible nodal displacements of the structure, we encounter various eigenvalues which are numerically close to zero:

returns:

-6.237207e-09
-4.329203e-10
415459e-11
183017e-10
478545e-09
290895e-09
000000e+00
141592e+07
141593e+07
141593e+07

Six of these are the so-called “rigid body modes”, associated to the three rigid translations and three rigid rotations of the structure. They can be “removed” by applying certain constraints to the set of allowed displacements. Enforcement of these constraints can be done by passing the parameter apply_rigid_body_constraint=True when calculating the stiffness model:

stiffness, _, _ = s.stiffness(storage='dense', apply_rigid_body_constraint=True)

To see that the six near zero eigenvalues of the stiffness matrix have been removed by the rigid body constraints recalculate:

d, v = stiffness.eigs()
d

to obtain:

Interestingly, in this case, there still remains one eigenvalue that is much smaller than the remaining ones. We will deal with this eigenvalue later.

For now, even though the smallest eigenvalue is small, the resulting stiffness matrix is not singular, and therefore suitable for computing displacements. This time:

x = stiffness.displacements(f)
x

successfully calculates the resulting approximate displacements:

-2.20468248e-09, -2.20468248e-09,  2.20468248e-09,  2.20468248e-09
 1.77419161e-09, -1.77419161e-09, -5.75306493e-09,  5.75306493e-09
 4.02657501e-18, -4.02657481e-18,  4.02657460e-18, -4.02657419e-18

successfully calculates the resulting approximate displacements, which can be visualized, after much enlargement, along with the applied forces in the figure

in which the forces are in green and the vectors indicating the resulting displacement are in yellow. This figure is generated by the code:

X = f
Y = 5e7*x

plotter = MatplotlibPlotter()
plotter.plot(s)
plotter.plot_arrows(s.nodes, f, color='g')
plotter.plot_arrows(s.nodes, 5e7*x, color='y')
fig, ax = plotter.get_handles()
ax.view_init(90,-90)
ax.axis('off')
plt.show()

Back to the small eigenvalue, which is sometimes associated with what is called a soft mode, in this case it appeared because the structure is planar, and its ball joints offer little resistance to out-of-plane forces. Indeed, the eigenvector associated with the eigenvalue is visualed by the following code:

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
plotter.plot_arrows(s.nodes, 0.2*v[:,0].reshape((3, 4), order='F'), color='y')
ax.view_init(30,-60)
ax.axis('off')
ax.axis('equal')
plt.show()

as the following figure:

in which it can be seen to be a pair of “couples” in the out-of-plane z-direction.

Constraining the node displacements to be planar "eliminates" such mode, as in:

stiffness, _, _ = s.stiffness(storage='dense', apply_rigid_body_constraint=True, apply_planar_constraint=True)

resulting in a structure in which:

d, v = stiffness.eigs()
d

equals:

indicating that there are no soft modes.

Of course one should expect no impact in the displacements if the forces do not have out-of-plane components and:

x = stiffness.displacements(f)
x

indeed returns displacements that are very similar to the ones calculated before, because the external forces do not have any non-zero z-components.

Node constraints

The above structures were free in space. A structure can be attached to a point in a reference frame by adding constraints to nodes.

Start by rebuilding the example structure used above:

from tnsgrt.structure import Structure
s = Structure(nodes, members, number_of_strings=4)

Fixed node constraints

The simplest type of node constraint is a fixed node constraint. Fixed node constraints mean that there can be no displacement of the node. It can be added to a node by setting the property 'constraint' to a tnsgrt.stiffness.NodeConstraint(). For example:

from tnsgrt.stiffness import NodeConstraint
s.set_node_properties([0, 3], 'constraint', NodeConstraint())

states that nodes 0 and 3 are to be considered fixed. Fixed node constraints are represented by spheres when you plot the structure:

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
ax.view_init(90,-90,0)
ax.axis('off')
ax.axis('equal')
plt.show()

Equilibrium in a structure with constraints is often different than equilibrium in a free structure. For example:

reactions = s.equilibrium()

produces the force coefficients

s.member_properties['lambda_']

  1.000000e+00
  1.000000e+00
  1.000000e+00
  1.062277e-10
 -1.000000e+00
 -1.000000e+00
Name: lambda\_, dtype: float64

and non-zero reaction forces at nodes 0 and 3

reactions[:, [0, 3]]

array([[ 1.00000000e+00, -1.00000000e+00],
       [ 2.04629879e-15,  2.04105833e-15],
       [ 0.00000000e+00,  0.00000000e+00]])

Note how the reaction forces replace the bottom string force (member 3), which now has a zero force coefficient. This behavior is the result of the equilibrium calculation that is set to minimize the sum of the force coefficients in all members.

The internal member force vectors can be calculated by multiplying the force coefficients by the member vectors as in

f_member = s.member_properties['lambda_'].values * s.get_member_vectors()

Those internal forces are plotted (in green) along with the structure and the reactions (in yellow) in the next figure:

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
plotter.plot_arrows(s.nodes[:,s.members[0,:]], .2*f_member, color='g')
plotter.plot_arrows(s.nodes[:,s.members[1,:]], -.2*f_member, color='g')
plotter.plot_arrows(s.nodes, .2*reactions, color='y')
ax.view_init(90,-90,0)
ax.axis('off')
ax.axis('equal')
plt.show()

Constraints are key in the presence of external forces. For example, consider the set of forces:

f = np.zeros((3, 4))
fz = np.array([[0, 1, 0]]).transpose()
f[:, [1,2]] = -fz
f

array([[ 0.,  0.,  0.,  0.],
       [ 0., -1., -1.,  0.],
       [ 0.,  0.,  0.,  0.]])

that represent a compressive load applied at the top nodes 1 and 2. It would be possible for the unconstrained structure to be in equilibrium with such forces. This is not a problem for the structure with constraints, for which

reactions = s.equilibrium(f)

produces the force coefficients

s.member_properties['lambda_']

 -1.780101e-09
  1.000000e+00
 -1.780101e-09
  1.777714e-10
 -1.000000e+00
 -1.000000e+00
Name: lambda\_, dtype: float64

and the reactions at the nodes 0 and 3:

reactions[:, [0, 3]]

array([[ 1., -1.],
       [ 1.,  1.],
       [ 0.,  0.]])

which are plotted along with the structure and the reactions in the next figure:

f_member = s.member_properties['lambda_'].values * s.get_member_vectors()

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
plotter.plot_arrows(s.nodes[:,s.members[0,:]], .2*f_member, color='g')
plotter.plot_arrows(s.nodes[:,s.members[1,:]], -.2*f_member, color='g')
plotter.plot_arrows(s.nodes, .2*reactions, color='y')
plotter.plot_arrows(s.nodes, .2*f, color='b')
ax.view_init(90,-90,0)
ax.axis('off')
ax.axis('equal')
plt.show()

In the presence of external forces, three strings now have zero force coefficients, and the reactions directly oposes the bar forces.

Node constraints with 2 degrees of freedom

It is also possible to impose constraints while preserving certain degrees of freedom. For example:

s.set_node_properties(3, 'constraint', NodeConstraint(constraint=np.array([[0, 1, 0]])))

constrain the displacements of node 3 to lie in the plane orthogonal to the y-axis. Node 0 continues to have a fixed node constraint.

When the structure is plotted, this constraint is represented by a small rectangle in the normal plane. The next figure is tilted a bit to show that.

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
ax.view_init(120,-90)
ax.axis('off')
ax.axis('equal')
plt.show()

In this type of constraint, node displacements are allowed only on the constraint plane, with the node sliding along it.

In the presence of the same compressive forces as before, equilibrium

reactions = s.equilibrium(f)

results in the member force coefficients

s.member_properties['lambda_']

 -4.454949e-11
  1.000000e+00
 -4.454694e-11
  1.000000e+00
 -1.000000e+00
 -1.000000e+00
Name: lambda\_, dtype: float64

and reactions at nodes 0 and 3

reactions[:, [0,3]]

array([[6.14180213e-15, 0.00000000e+00],
       [1.00000000e+00, 1.00000000e+00],
       [0.00000000e+00, 0.00000000e+00]])

The corresponding internal forces and reactions are plotted next:

f_member = s.member_properties['lambda_'].values * s.get_member_vectors()

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
plotter.plot_arrows(s.nodes[:,s.members[0,:]], .2*f_member, color='g')
plotter.plot_arrows(s.nodes[:,s.members[1,:]], -.2*f_member, color='g')
plotter.plot_arrows(s.nodes, .2*reactions, color='y')
plotter.plot_arrows(s.nodes, .2*f, color='b')
ax.view_init(120,-90,0)
ax.axis('off')
ax.axis('equal')
plt.show()

Note that because the constraint at node 3 can only apply forces in the y-direction, the bottom strings carry forces, which the two vertical strings do not.

Node constraints with one degree of freedom

Finally, it is also possible to constrain motion along a line, as in

s.set_node_properties(3, 'constraint', NodeConstraint(displacement=np.array([[1, 0, 0]]).transpose()))

which is illustrated in the next figure by the line at node 3:

plotter = MatplotlibPlotter()
plotter.plot(s)
_, ax = plotter.get_handles()
ax.view_init(120,-90)
ax.axis('off')
ax.axis('equal')
plt.show()

Equilibrium is the same as in the case of the planar constraint.