import panel as pn import scipy.optimize import scipy.linalg import pandas as pd import numpy as np from matplotlib.figure import Figure x_points = [0.] * 13 d_points = [0.] * 13 def elm_k(L, d1, d2, alpha): b = (d2 - d1) / L return 2 * alpha / (35*L**3) * np.array( [ [ 6*(11*L**4*b**4 + 49*L**3*b**3*d1 + 84*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4), L*(19*L**4*b**4 + 84*L**3*b**3*d1 + 147*L**2*b**2*d1**2 + 140*L*b*d1**3 + 105*d1**4), 6*(-11*L**4*b**4 - 49*L**3*b**3*d1 - 84*L**2*b**2*d1**2 - 70*L*b*d1**3 - 35*d1**4), L*(47*L**4*b**4 + 210*L**3*b**3*d1 + 357*L**2*b**2*d1**2 + 280*L*b*d1**3 + 105*d1**4) ], [ L*(19*L**4*b**4 + 84*L**3*b**3*d1 + 147*L**2*b**2*d1**2 + 140*L*b*d1**3 + 105*d1**4), 2*L**2*(3*L**4*b**4 + 14*L**3*b**3*d1 + 28*L**2*b**2*d1**2 + 35*L*b*d1**3 + 35*d1**4), L*(-19*L**4*b**4 - 84*L**3*b**3*d1 - 147*L**2*b**2*d1**2 - 140*L*b*d1**3 - 105*d1**4), L**2*(13*L**4*b**4 + 56*L**3*b**3*d1 + 91*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4) ], [ 6*(-11*L**4*b**4 - 49*L**3*b**3*d1 - 84*L**2*b**2*d1**2 - 70*L*b*d1**3 - 35*d1**4), L*(-19*L**4*b**4 - 84*L**3*b**3*d1 - 147*L**2*b**2*d1**2 - 140*L*b*d1**3 - 105*d1**4), 6*(11*L**4*b**4 + 49*L**3*b**3*d1 + 84*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4), L*(-47*L**4*b**4 - 210*L**3*b**3*d1 - 357*L**2*b**2*d1**2 - 280*L*b*d1**3 - 105*d1**4) ], [ L*(47*L**4*b**4 + 210*L**3*b**3*d1 + 357*L**2*b**2*d1**2 + 280*L*b*d1**3 + 105*d1**4), L**2*(13*L**4*b**4 + 56*L**3*b**3*d1 + 91*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4), L*(-47*L**4*b**4 - 210*L**3*b**3*d1 - 357*L**2*b**2*d1**2 - 280*L*b*d1**3 - 105*d1**4), 2*L**2*(17*L**4*b**4 + 77*L**3*b**3*d1 + 133*L**2*b**2*d1**2 + 105*L*b*d1**3 + 35*d1**4) ] ] ) lengthInput = pn.widgets.FloatInput(name="Bow Length (mm)", value=700, step=1, start=600, end=750) lengthConstInput = pn.widgets.FloatInput(name="Length Constant Diameter/Width (mm)", value=110, step=1, start=60, end=200) dButtInput = pn.widgets.FloatInput(name="Diameter/Width at Butt (mm)", value=8.6, step=0.1, start=6, end=12) dTipInput = pn.widgets.FloatInput(name="Diameter/Width at Tip (mm)", value=5.6, step=0.1, start=3, end=10) shapeInput = pn.widgets.Select(name="Shape", options={ "Round": 0.0490874, "Octagonal": 0.0547379 }) modulusInput = pn.widgets.FloatInput(name="Elastic Modulus (GPa)", value=30, step=0.1, start=10, end=40 ) taperOutput = pn.pane.DataFrame( pd.DataFrame(data = {"x": x_points, "d": d_points}), border=3, width=400, col_space=75, justify="left", formatters=[ lambda x: f"{x:.0f}", lambda x: f"{x:.2f}" ] ) deflectionOutput = pn.pane.Str("") fig0 = Figure(figsize=(6, 4)) plotOutput = pn.pane.Matplotlib(fig0, dpi=144) calculateButton = pn.widgets.Button(name="Calculate", width=300) def calc_taper(): length = lengthInput.value length_constant = lengthConstInput.value d_butt = dButtInput.value d_head = dTipInput.value r = scipy.optimize.root( lambda r: length_constant * (r**12 - 1) / (r - 1) - length, 22 ).x[0] for i in range(13): if i == 0: d_points[i] = d_butt else: x_points[i] = length_constant * (r**i - 1) / (r - 1) d_points[i] = d_butt + (d_head - d_butt) * (i - 1.) / 10. if i == 12: d_points[i] = d_head df = pd.DataFrame(data = { "x": x_points, "d": d_points }) taperOutput.object = df def remesh(): x_nodes = [] d_nodes = [] x_s2 = x_points[11] x_s1 = x_s2 - 575 x_l = x_s2 - 575 / 2 nid_s1 = -1 # storage for node ID of support #1 nid_s2 = -1 # storage for node ID of support #2 nid_l = -1 # storage for node ID of load application tol = lambda xa, xb: abs(xa - xb) < 1e-3 inside = lambda x, xa, xb: (x - xa) * (x - xb) < 0 for x1, x2, d1, d2 in zip( x_points, x_points[1:], d_points, d_points[1:]): x_nodes.append(x1) d_nodes.append(d1) if tol(x_s1, x1): nid_s1 = len(x_nodes) - 1 elif inside(x_s1, x1, x2): x_nodes.append(x_s1) d_nodes.append(d1 + (x_s1 - x1) / (x2 - x1) * (d2 - d1)) nid_s1 = len(x_nodes) - 1 if tol(x_s2, x1): nid_s2 = len(x_nodes) - 1 elif inside(x_s2, x1, x2): x_nodes.append(x_s2) d_nodes.append(d1 + (x_s2 - x1) / (x2 - x1) * (d2 - d1)) nid_s2 = len(x_nodes) - 1 if tol(x_l, x1): nid_l = len(x_nodes) - 1 elif inside(x_l, x1, x2): x_nodes.append(x_l) d_nodes.append(d1 + (x_l - x1) / (x2 - x1) * (d2 - d1)) nid_l = len(x_nodes) - 1 x_nodes.append(x_points[-1]) d_nodes.append(d_points[-1]) return x_nodes, d_nodes, nid_s1, nid_s2, nid_l def calc_stiffness_matrix(x_nodes, d_nodes): k_model= np.zeros((2 * len(x_nodes), 2 * len(x_nodes))) modulus = modulusInput.value * 1000 shape_constant = shapeInput.value for i, (x1, x2, d1, d2) in enumerate( zip(x_nodes, x_nodes[1:], d_nodes, d_nodes[1:])): # Each element connects the two adjacent nodes k_elm = elm_k( L = x2 - x1, d1 = d1, d2 = d2, alpha = shape_constant * modulus ) for ii in range(4): for jj in range(4): k_model[i * 2 + ii, i * 2 + jj] += k_elm[ii,jj] return k_model def solve(k_model, p_model, x_nodes, nid_s1, nid_s2): mask = [i for i, _ in enumerate(p_model) if i != nid_s1 * 2 and i != nid_s2 * 2] p_const = p_model[mask] k_const = k_model[mask, :] k_const = k_const[:, mask] d_const = scipy.linalg.solve(k_const, p_const) d_model = np.zeros(2 * len(x_nodes)) d_model[mask] = d_const return d_model def draw_deflection(x_nodes, d_model): fig0.clear() ax0 = fig0.subplots() ax0.clear() ax0.plot(x_nodes, d_model[0::2]) ax0.grid() ax0.set_title("Deflection, 2 lb load, 575mm span") ax0.set_xlabel("x") ax0.set_ylabel("Vertical Deflection") plotOutput.object = fig0 def calc_deflection(): x_nodes, d_nodes, nid_s1, nid_s2, nid_l = remesh() k_model = calc_stiffness_matrix(x_nodes, d_nodes) # Calculate the load vector p_model = np.zeros(2 * len(x_nodes)) p_model[nid_l * 2] = -8.9075 # 2 lb in N d_model = solve(k_model, p_model, x_nodes, nid_s1, nid_s2) draw_deflection(x_nodes, d_model) deflect_load = -d_model[nid_l * 2] / 25.4 * 1000 deflectionOutput.object = f"Deflection: {deflect_load:.3f}" def calculate(event): calc_taper() calc_deflection() calculateButton.on_click(calculate) pn.Column( pn.FlexBox(*[lengthInput, lengthConstInput, dButtInput, dTipInput, shapeInput, modulusInput]), calculateButton, taperOutput, deflectionOutput, plotOutput ).servable(target='taper_input_panel');