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