Coverage for docs\source\examples\beam.py: 0%

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 

7# from gpkit import parse_variables, Model, ureg 

8# from gpkit.small_scripts import mag 

9 

10eps = 2e-4 # has to be quite large for consistent cvxopt printouts; 

11 # normally you'd set this to something more like 1e-20 

12 

13 

14class Beam(Model): 

15 """Discretization of the Euler beam equations for a distributed load. 

16 

17 Variables 

18 --------- 

19 EI [N*m^2] Bending stiffness 

20 dx [m] Length of an element 

21 L 5 [m] Overall beam length 

22 

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 

29 

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 

37 

38 Upper Unbounded 

39 --------------- 

40 w_tip 

41 

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 } 

68 

69 

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

73# w_gp = sol("w") # deflection along beam 

74# 

75# L, EI, q = sol("L"), sol("EI"), sol("q") 

76# x = np.linspace(0, mag(L), len(q))*ureg.m # position along beam 

77# q = q[0] # assume uniform loading for the check below 

78# w_exact = q/(24*EI) * x**2 * (x**2 - 4*L*x + 6*L**2) # analytic soln 

79# assert max(abs(w_gp - w_exact)) <= 1.1*ureg.cm 

80 

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