Coverage for src / sdynpy / fem / sdynpy_beam.py: 85%
104 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Functions for creating mass and stiffness matrices for beam structures.
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import numpy as np
25import scipy.linalg as la
28def beamkm(node_coords, element_connectivity, bend_direction_1,
29 ae, jg, ei1, ei2, mass_per_length, tmmi_per_length):
30 '''Compute mass and stiffness matrices for the given beam model
32 This function computes the mass and stiffness matrices for the beam model
33 defined by the given nodal coordinates, element connectivity, orientation
34 parameters, and elastic properties.
36 Parameters
37 ----------
38 node_coords : ndarray
39 An nx3 array where node_coords[i,:] is the x,y,z coordinate of the i-th
40 node. The number of rows is equal to the number of nodes in the model.
41 element_connectivity : ndarray
42 An mx2 array where element_connectivity[i,:] are the row indices of
43 node_coords corresponding to the nodes in the i-th element. The number
44 of rows is equal to the number of elements in the model
45 bend_direction_1 : ndarray
46 An mx3 array where bend_direction_1[i,:] is the x,y,z coordinates of a
47 vector that defines the first bending axis of the beam. The second
48 bending axis is computed from the cross product of the beam axis and
49 the first bending axis. The number of rows is equal to the number of
50 elements in the model
51 ae : ndarray
52 A length m array where ae[i] is the axial stiffness of element i. This
53 can be computed from the area of the element (A) times the elastic
54 modulus (E). The length of the array should be equal to the number of
55 elements in the model.
56 jg : ndarray
57 A length m array where jg[i] is the torsional stiffness of element i.
58 This can be computed from the torsional constant (J) times the shear
59 modulus (G). The length of the array should be equal to the number of
60 elements in the model.
61 ei1 : ndarray
62 A length m array where ei1[i] is the bending stiffness about axis 1 of
63 element i. This can be computed from the second moment of area of the
64 beam's cross section about bending axis 1 (I1) times the elastic
65 modulus (E). The length of the array should be equal to the number of
66 elements in the model.
67 ei2 : ndarray
68 A length m array where ei2[i] is the bending stiffness about axis 2 of
69 element i. This can be computed from the second moment of area of the
70 beam's cross section about bending axis 2 (I2) times the elastic
71 modulus (E). The length of the array should be equal to the number of
72 elements in the model.
73 mass_per_length : ndarray
74 A length m array where mass_per_length[i] is the linear density of
75 element i. This can be computed from the density of the material times
76 the cross sectional area of the beam. The length of the array should
77 be equal to the number of elements in the model.
78 tmmi_per_length : ndarray
79 A length m array where tmmi_per_length[i] is the torsional mass moment
80 of inertia per unit length of element i. This can be computed from the
81 density of the material times the polar moment of inertia of the beam's
82 cross-section. The length of the array should be equal to the number
83 of elements in the model.
85 Returns
86 -------
87 K : ndarray
88 The stiffness matrix of the beam model, which will have size 6nx6n
89 where n is the number of nodes in the model.
90 M : ndarray
91 The mass matrix of the beam model, which will have size 6nx6n where n
92 is the number of nodes in the model.
94 Notes
95 -----
96 The degrees of freedom in the mass and stiffness matrices will be ordered
97 as follows:
98 [disp_x_0, disp_y_0, disp_z_0, rot_x_0, rot_y_0, rot_z_0,
99 disp_x_1, disp_y_1, disp_z_1, rot_x_1, rot_y_1, rot_z_1,
100 ...
101 disp_x_n, disp_y_n, disp_z_n, rot_x_n, rot_y_n, rot_z_n]
102 '''
103 # Validate inputs
104 try:
105 number_of_nodes = node_coords.shape[0]
106 if node_coords.ndim != 2:
107 raise ValueError(
108 'node_coords should be a 2D array with shape nx3 where n is the number of nodes in the model')
109 if node_coords.shape[1] != 3:
110 raise ValueError(
111 'node_coords should have shape nx3 where n is the number of nodes in the model')
112 except AttributeError:
113 raise ValueError('node_coords should be a numpy.ndarray')
114 try:
115 number_of_elements = element_connectivity.shape[0]
116 if element_connectivity.ndim != 2:
117 raise ValueError(
118 'element_connectivity should be a 2D array with shape nx2 where n is the number of elements in the model')
119 if element_connectivity.shape[1] != 2:
120 raise ValueError(
121 'element_connectivity should have shape nx2 where n is the number of elements in the model')
122 except AttributeError:
123 raise ValueError('element_connectivity should be a numpy.ndarray')
124 # Check that the element properties that have been passed are the correct
125 # shape
126 try:
127 if bend_direction_1.shape != (number_of_elements, 3):
128 raise ValueError(
129 'bend_direction_1 should be a 2D array with shape (element_connectivity.shape[0],3)')
130 except AttributeError:
131 raise ValueError('bend_direction_1 should be a numpy.ndarray')
132 for val in [ae, jg, ei1, ei2, mass_per_length, tmmi_per_length]:
133 try:
134 if val.shape != (number_of_elements,):
135 raise ValueError(
136 'Element Properties (ae,jg,ei1,ei2,mass_per_length,tmmi_per_length) should be 1D arrays with length element_connectivity.shape[0]')
137 except AttributeError:
138 raise ValueError(
139 'Element Properties (ae,jg,ei1,ei2,mass_per_length,tmmi_per_length) should be numpy.ndarray')
140 # Initialize
141 K = np.zeros((6 * number_of_nodes, 6 * number_of_nodes))
142 M = np.zeros((6 * number_of_nodes, 6 * number_of_nodes))
144 # Initialize Element Matrices
145 PA = np.array(
146 [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
147 [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])
148 PT = np.array(
149 [[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
150 [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]])
151 PB1 = np.array(
152 [[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
153 [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
154 [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
155 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])
156 PB2 = np.array(
157 [[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
158 [0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0],
159 [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
160 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0]])
161 # Loop through all elements and assemble mass and stiffness matrices
162 for i, (node_indices, bd1, AE, JG, EI1, EI2, rho, rhoT) in enumerate(zip(element_connectivity,
163 bend_direction_1, ae, jg, ei1, ei2,
164 mass_per_length, tmmi_per_length)):
165 # Get element length
166 dx = node_coords[node_indices[1], :] - node_coords[node_indices[0], :]
167 L = np.linalg.norm(dx)
168 # Compute Element Matrices
169 KelemA = AE / L * np.array([[1, -1],
170 [-1, 1]])
171 KelemT = JG / L * np.array([[1, -1],
172 [-1, 1]])
173 KelemB1 = EI1 / L**3 * np.array([[12, 6 * L, -12, 6 * L],
174 [6 * L, 4 * L**2, -6 * L, 2 * L**2],
175 [-12, -6 * L, 12, -6 * L],
176 [6 * L, 2 * L**2, -6 * L, 4 * L**2]])
177 KelemB2 = EI2 / L**3 * np.array([[12, 6 * L, -12, 6 * L],
178 [6 * L, 4 * L**2, -6 * L, 2 * L**2],
179 [-12, -6 * L, 12, -6 * L],
180 [6 * L, 2 * L**2, -6 * L, 4 * L**2]])
181 Kelem = PA.T @ KelemA @ PA + PT.T @ KelemT @ PT + PB1.T @ KelemB1 @ PB1 + PB2.T @ KelemB2 @ PB2
183 MelemA = rho * L * np.array([[1 / 3, 1 / 6],
184 [1 / 6, 1 / 3]])
185 MelemT = rhoT * L * np.array([[1 / 3, 1 / 6],
186 [1 / 6, 1 / 3]])
187 MelemB1 = rho * L * np.array([[13 / 35, 11 / 210 * L, 9 / 70, -13 / 420 * L],
188 [11 / 210 * L, 1 / 105 * L**2, 13 / 420 * L, -1 / 140 * L**2],
189 [9 / 70, 13 / 420 * L, 13 / 35, -11 / 210 * L],
190 [-13 / 420 * L, -1 / 140 * L**2, -11 / 210 * L, 1 / 105 * L**2]])
191 MelemB2 = rho * L * np.array([[13 / 35, 11 / 210 * L, 9 / 70, -13 / 420 * L],
192 [11 / 210 * L, 1 / 105 * L**2, 13 / 420 * L, -1 / 140 * L**2],
193 [9 / 70, 13 / 420 * L, 13 / 35, -11 / 210 * L],
194 [-13 / 420 * L, -1 / 140 * L**2, -11 / 210 * L, 1 / 105 * L**2]])
195 Melem = PA.T @ MelemA @ PA + PT.T @ MelemT @ PT + PB1.T @ MelemB1 @ PB1 + PB2.T @ MelemB2 @ PB2
196 # Compute Directions
197 d0 = dx / np.linalg.norm(dx)
198 d1 = bd1 / np.linalg.norm(bd1)
199 d2 = np.cross(d0, d1)
200 d2 = d2 / np.linalg.norm(d2)
201 d1 = np.cross(d2, d0)
202 d1 = d1 / np.linalg.norm(d1)
203 C = np.array([d0, d1, d2])
205 A = np.zeros((12, 12))
206 A[0:3, 0:3] = C
207 A[3:6, 3:6] = C
208 A[6:9, 6:9] = C
209 A[9:12, 9:12] = C
211 i1 = np.arange(6 * node_indices[0], 6 * (node_indices[0] + 1))
212 i2 = np.arange(6 * node_indices[1], 6 * (node_indices[1] + 1))
214 K[np.concatenate((i1, i2))[:, np.newaxis], np.concatenate((i1, i2))] += A.T @ Kelem @ A
215 M[np.concatenate((i1, i2))[:, np.newaxis], np.concatenate((i1, i2))] += A.T @ Melem @ A
217 return K, M
220'''
221node_coords = np.array([np.linspace(0,1,10),
222 np.zeros(10),
223 np.zeros(10)]).T
225element_connectivity = np.array([np.arange(0,9),np.arange(1,10)]).T
227bend_direction_1 = np.array([np.zeros(element_connectivity.shape[0]),np.zeros(element_connectivity.shape[0]),np.ones(element_connectivity.shape[0])]).T
229b = 0.01
230h = 0.01
231E = 69e9
232nu = 0.33
233rho = 2830
234I1 = b*h**3/12
235I2 = b**3*h/12
236G = E/(2*(1-nu))
237J = I1+I2
238Ixx_per_L = (1/12)*rho*b*h*(b**2+h**2)
240ae = b*h*E * np.ones(element_connectivity.shape[0])
241jg = J*G * np.ones(element_connectivity.shape[0])
242ei1 = E*I1 * np.ones(element_connectivity.shape[0])
243ei2 = E*I2 * np.ones(element_connectivity.shape[0])
244mass_per_length = rho*b*h * np.ones(element_connectivity.shape[0])
245tmmi_per_length = Ixx_per_L * np.ones(element_connectivity.shape[0])
247K,M = beamkm(node_coords,element_connectivity,bend_direction_1,ae,jg,ei1,ei2,mass_per_length,tmmi_per_length)
249# Perform an eigenanalysis
250lam,phi = la.eigh(K,M)
251lam[0:6] = 0
252fn = np.sqrt(lam)/(2*np.pi)
253'''
256def rect_beam_props(E, rho, nu, b, h=None, nelem=None):
257 '''Return beam keyword dictionary for a rectangular beam.
259 Returns parameters that can be used in the beamkm function for a beam of
260 uniform rectangular cross-section.
262 Parameters
263 ----------
264 E : float
265 Young's Modulus
266 rho : float
267 Density
268 nu : float
269 Poisson's Ratio
270 b : float
271 Thickness of the beam in the direction of bend_direction_1
272 h : float
273 Thickness of the beam in the direction of bend_direction_2. If it is
274 not specified, a square cross section is assumed.
275 nelem : int
276 Number of elements. This will repeat the values in the dict nelem
277 times, as is required by beamkm. If not specified, only a single value
278 will be output for each entry in the dict.
280 Returns
281 -------
282 kwargs : dict
283 A dictionary that can be used in beamkm(**kwargs).
285 Notes
286 -----
287 Steel : SI
288 E = 200e9 # [N/m^2],
289 nu = 0.25 # [-],
290 rho = 7850 # [kg/m^3]
291 Aluminum : SI
292 E = 69e9 # [N/m^2],
293 nu = 0.33 # [-],
294 rho = 2830 # [kg/m^3]
295 '''
296 if h is None:
297 h = b
298 A = b * h
299 I1 = b * h**3 / 12
300 I2 = b**3 * h / 12
301 G = E / (2 * (1 - nu))
302 J = I1 + I2
303 Ixx_per_L = (1 / 12) * rho * b * h * (b**2 + h**2)
304 return_dict = {}
305 return_dict['ae'] = A * E
306 return_dict['jg'] = J * G
307 return_dict['ei1'] = E * I1
308 return_dict['ei2'] = E * I2
309 return_dict['mass_per_length'] = rho * A
310 return_dict['tmmi_per_length'] = Ixx_per_L
311 if nelem is not None:
312 for key in return_dict:
313 return_dict[key] = return_dict[key] * np.ones((nelem,), float)
314 return return_dict
317def beamkm_2d(length, width, height, nnodes, E, rho, nu, axial=True):
318 '''A simplified 2D beam with uniform cross section and linear materials
320 Parameters
321 ----------
322 length : float
323 Length of the beam
324 width : float
325 Width of the beam, perpendicular to the bending plane
326 height : float
327 Height of the beam, in the bending direction.
328 nnodes : int
329 Number of nodes in the beam (number of elements + 1)
330 E : float
331 Elastic modulus
332 rho : float
333 Density
334 nu : float
335 Poisson's Ratio
336 axial : bool
337 Keep axial motions (default = True)
339 Returns
340 -------
341 K : np.ndarray
342 Stiffness matrix of the beam
343 M : np.ndarray
344 Mass matrix of the beam
346 Notes
347 -----
348 Dof ordering is axial translation, bending translation, and bending
349 rotation repeated for each node. Axial is X, bending is Z, and rotation is
350 about Y. If axial == False, Dof ordering is bending translation, bending
351 rotation.
352 '''
353 node_coords = np.linspace(0, length, nnodes)
354 coords_3d = np.zeros((nnodes, 3))
355 coords_3d[:, 0] = node_coords
356 node_indices = np.arange(nnodes)
357 elems = np.array((node_indices[:-1], node_indices[1:])).T
358 bend_direction_1 = np.zeros((len(elems), 3))
359 bend_direction_1[:, 2] = 1
361 props = rect_beam_props(E, rho, nu, width, height, nelem=len(elems))
363 K, M = beamkm(coords_3d, elems, bend_direction_1, **props)
365 # Eliminate everything except in bend direction 1. Keep translation in X, Z,
366 # and rotation about y
367 if axial:
368 keep_dofs = ([0, 2, 4] + np.arange(nnodes)[:, np.newaxis] * 6).flatten()
369 else:
370 keep_dofs = ([2, 4] + np.arange(nnodes)[:, np.newaxis] * 6).flatten()
372 K = K[keep_dofs[:, np.newaxis], keep_dofs[np.newaxis, :]]
373 M = M[keep_dofs[:, np.newaxis], keep_dofs[np.newaxis, :]]
374 return (K, M)