0

I need to calculate the vertical natural frequency and the horizontal natural frequency of a beam. In order to do this I have been given two values of the vertical and horizontal natural frequency which I need to approximate as closey as possible. Now in order to do this I have to use the moment of inertia in the y and z direction and the surface area. So the values that I need to calculate in this instance are the width(B), the height(H) and the thickness(t_h and t_v). For this I have used two solvers. The first one computes the values for B and H and the second one computes the values for t_h and t_v. Now the first solver works just fine and gives me numerical values. However, the second one somehow does not and only gives me []. Below I have put the code in that I have used. I cannot seem to understand what I am doing wrong in this code and I'm hoping that somebody can help me with his, because I have been at it for a couple of days and cannot seem to figure it out.

import numpy as np
from sympy import symbols, Eq, solve

w_2_vb_g=4.341
w_2_hb_g=2.601
E=210E9
rho=7700
L=100
kL_2=7.85
kL_3=11.00

B,H,t_h,t_v=symbols('B H t_h t_v', real=True, positive=True)
w_2_vb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*H*B**3))/(rho*(B*H)))**(1/2))), w_2_vb_g)
w_2_hb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*B*H**3))/(rho*(H*B)))**(1/2))), w_2_hb_g)
Solution=solve((w_2_hb_rec,w_2_vb_rec), (B,H))
print(Solution)

print("Solutions:")
for var, val in Solution.items():
    print(f"{var}: {val}")

list_from_dict = list(Solution.items())
print("List from dict:", list_from_dict)

B=round(float(list_from_dict[1][1]*1.03),3)
H=round(float(list_from_dict[0][1]*1.03),3)

B_list=np.array([t_h, (H-2*t_h), B, t_v])
H_list=np.array([B, t_v, t_h, (H-2*t_h)])
D0_list=np.array([0, (0.5*B-0.5*t_v), (0.5*H-0.5*t_h),0])

w_2_vb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[0]*H_list[0]**3)+B_list[0]*H_list[0]*D0_list[0]**2)+2*(((1/12)*B_list[1]*H_list[1]**3)+B_list[1]*H_list[1]*D0_list[1]**2))/(rho*(((2*(B_list[0]*H_list[0])))+(2*(B_list[1]*H_list[1]))))**(1/2))))), w_2_vb_g)
w_2_hb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[2]*H_list[2]**3)+B_list[2]*H_list[2]*D0_list[2]**2)+2*(((1/12)*B_list[3]*H_list[3]**3)+B_list[3]*H_list[3]*D0_list[3]**2))/(rho*(((2*(H_list[2]*B_list[2])))+(2*(H_list[3]*H_list[3]))))**(1/2))))), w_2_hb_g)
Solution_beam=solve((w_2_hb,w_2_vb), (t_h,t_v))
print(Solution_beam)

Below You can see the output that my console gives:

{H: 0.279980237079456, B: 0.467279588297547}
Solutions:
H: 0.279980237079456
B: 0.467279588297547
List from dict: [(H, 0.279980237079456), (B, 0.467279588297547)]
[]

I have already tried different solver such as linsolve() or solve_undetermined_coeffs but neither has worked. This could be attributed to my code and that it doesn't work in this case but as far as I can see it should. I have also considered by splitting the equations in multiple parts so as to make it more comprehensive for the solver. However, this is not possible because I have not been given a starting value for the moment of inertia or the surface area which are the only two options in which it could've worked. Hopefully somebody can help, many thanks in advance!

1 Answer 1

0

As other questions on SO (for example this one) have also demonstrated, the symbolic solvers such as solve from sympy are not designed for solving such numerical problems as your issue here. This is due to multiple reasons in this case, e.g., the use of big numbers and that symbolic solvers can not handle non-linearity in general.

Defining the second function numerically results in calculating t_h and t_v:

def equations(vars):
    t_h, t_v = vars

    B_list = np.array([t_h, (H_val - 2 * t_h), B_val, t_v])
    H_list = np.array([B_val, t_v, t_h, (H_val - 2 * t_h)])
    D0_list = np.array([0, (0.5 * B_val - 0.5 * t_v), (0.5 * H_val - 0.5 * t_h), 0])

    term1 = (E * (2 * (((1 / 12) * B_list[0] * H_list[0] ** 3) + B_list[0] * H_list[0] * D0_list[0] ** 2) +
                    2 * (((1 / 12) * B_list[1] * H_list[1] ** 3) + B_list[1] * H_list[1] * D0_list[1] ** 2))) / \
            (rho * ((2 * (B_list[0] * H_list[0])) + (2 * (B_list[1] * H_list[1]))))

    term2 = (E * (2 * (((1 / 12) * B_list[2] * H_list[2] ** 3) + B_list[2] * H_list[2] * D0_list[2] ** 2) +
                    2 * (((1 / 12) * B_list[3] * H_list[3] ** 3) + B_list[3] * H_list[3] * D0_list[3] ** 2))) / \
            (rho * ((2 * (H_list[2] * B_list[2])) + (2 * (H_list[3] * H_list[3]))))

    eq_vb = ((((kL_2) ** 2) / (L ** 2)) * (term1 ** (1 / 2))) - w_2_vb_g
    eq_hb = ((((kL_2) ** 2) / (L ** 2)) * (term2 ** (1 / 2))) - w_2_hb_g

    return [eq_vb, eq_hb]

#initial guesses for solver
t_h_init, t_v_init = 0.01, 0.01

solution_beam = fsolve(equations, [t_h_init, t_v_init])

t_h, t_v = solution_beam

print(f"Computed values: t_h = {t_h:.4f}, t_v = {t_v:.4f}")
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.