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

19 statements  

« prev     ^ index     » next       coverage.py v7.4.0, created at 2024-01-03 16:49 -0500

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 

8 

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

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

11 

12 

13class Beam(Model): 

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

15 

16 Variables 

17 --------- 

18 EI [N*m^2] Bending stiffness 

19 dx [m] Length of an element 

20 L 5 [m] Overall beam length 

21 

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 

28 

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 

36 

37 Upper Unbounded 

38 --------------- 

39 w_tip 

40 

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 } 

67 

68 

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 

73 

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 

79 

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