# Coverage for docs/source/examples/beam.py: 100%

## 19 statements

, created at 2022-07-28 12:35 -0400

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 ----------------------------

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 >= th_base,

52 w >= w_base

53 ],

54 # below: trapezoidal integration to form a piecewise-linear

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

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()