Coverage for docs/source/examples/beam.py: 100%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2A simple beam example with fixed geometry. Solves the discretized
3Euler-Bernoulli beam equations for a constant distributed load
4"""
5import numpy as np
6from gpkit import parse_variables, Model
7from gpkit import parse_variables, Model, ureg
8from gpkit.small_scripts import mag
10eps = 2e-4 # has to be quite large for consistent cvxopt printouts;
11 # normally you'd set this to something more like 1e-20
14class Beam(Model):
15 """Discretization of the Euler beam equations for a distributed load.
17 Variables
18 ---------
19 EI [N*m^2] Bending stiffness
20 dx [m] Length of an element
21 L 5 [m] Overall beam length
23 Boundary Condition Variables
24 ----------------------------
25 V_tip eps [N] Tip loading
26 M_tip eps [N*m] Tip moment
27 th_base eps [-] Base angle
28 w_base eps [m] Base deflection
30 Node Variables of length N
31 --------------------------
32 q 100*np.ones(N) [N/m] Distributed load
33 V [N] Internal shear
34 M [N*m] Internal moment
35 th [-] Slope
36 w [m] Displacement
38 Upper Unbounded
39 ---------------
40 w_tip
42 """
43 @parse_variables(__doc__, globals())
44 def setup(self, N=4):
45 # minimize tip displacement (the last w)
46 self.cost = self.w_tip = w[-1]
47 return {
48 "definition of dx": L == (N-1)*dx,
49 "boundary_conditions": [
50 V[-1] >= V_tip,
51 M[-1] >= M_tip,
52 th[0] >= th_base,
53 w[0] >= w_base
54 ],
55 # below: trapezoidal integration to form a piecewise-linear
56 # approximation of loading, shear, and so on
57 # shear and moment increase from tip to base (left > right)
58 "shear integration":
59 V[:-1] >= V[1:] + 0.5*dx*(q[:-1] + q[1:]),
60 "moment integration":
61 M[:-1] >= M[1:] + 0.5*dx*(V[:-1] + V[1:]),
62 # slope and displacement increase from base to tip (right > left)
63 "theta integration":
64 th[1:] >= th[:-1] + 0.5*dx*(M[1:] + M[:-1])/EI,
65 "displacement integration":
66 w[1:] >= w[:-1] + 0.5*dx*(th[1:] + th[:-1])
67 }
70b = Beam(N=6, substitutions={"L": 6, "EI": 1.1e4, "q": 110*np.ones(6)})
71sol = b.solve(verbosity=0)
72print(sol.table(tables=["cost breakdown"], maxcolumns=6))
73w_gp = sol("w") # deflection along beam
75L, EI, q = sol("L"), sol("EI"), sol("q")
76x = np.linspace(0, mag(L), len(q))*ureg.m # position along beam
77q = q[0] # assume uniform loading for the check below
78w_exact = q/(24*EI) * x**2 * (x**2 - 4*L*x + 6*L**2) # analytic soln
79assert max(abs(w_gp - w_exact)) <= 1.1*ureg.cm
81PLOT = False
82if PLOT: # pragma: no cover
83 import matplotlib.pyplot as plt
84 x_exact = np.linspace(0, L, 1000)
85 w_exact = q/(24*EI) * x_exact**2 * (x_exact**2 - 4*L*x_exact + 6*L**2)
86 plt.plot(x, w_gp, color='red', linestyle='solid', marker='^',
87 markersize=8)
88 plt.plot(x_exact, w_exact, color='blue', linestyle='dashed')
89 plt.xlabel('x [m]')
90 plt.ylabel('Deflection [m]')
91 plt.axis('equal')
92 plt.legend(['GP solution', 'Analytical solution'])
93 plt.show()