<jupyter_start><jupyter_text>적용사례 : 오일러 베르누이 외팔보Application Example : Euler-Bernoulli Cantilever<jupyter_code>import ode_solver
import scipy.integrate as si<jupyter_output><empty_output><jupyter_text>오일러 베르누이 외팔보Euler-Bernoulli Cantilever Ref : Wikipedia contributors, 'Euler–Bernoulli beam theory', Wikipedia, The Free Encyclopedia, 19 November 2018, 22:12 UTC, https://en.wikipedia.org/w/index.php?title=Euler%E2%80%93Bernoulli_beam_theory&oldid=869647128 [accessed 24 November 2018] 외팔보는 한쪽 끝에서 처짐과 기울기가 0인 보 이다.A cantilever is a beam with both zero deflection and rotation at one end. 오일러 베르누이 보 이론은 분포하중과 보의 처짐은 다음과 같은 관계를 가진다고 가정한다.Euler-Bernoulli beam theory assumes that the deflection of a beam and the distributed load has following relationship. $$\frac{d^2}{dx^2} \left( EI \frac{d^2w(x)}{dx^2} \right) = q(x)$$ 여기서 $w(x)$와 $q(x)$는 각각 보의 위치 $x$에서의 $z$방향 처짐과 분포하중이다.Here, $w(x)$ and $q$ are, respectively, $z$ directional deflection and distributed load at a location $x$ of the beam. 단순화 하기 위해 $EI$는 일정하다고 가정하자.To simplify, let's assume $EI$ is constant. $$E[Nm^{-2}]I[m^4] \frac{d^4w(x)}{dx^4}[m^{-3}] = q(x)[Nm^{-1}]$$ 상태 변수 $\mathbb{r}(x)$ 를 다음과 같이 정해 보자.Let state variable $\mathbb{r}(x)$ be as follows. $$\mathbb{r}(x) = \begin{pmatrix}r_0\\r_1\\r_2\\r_3\\\end{pmatrix} = \begin{pmatrix}w(x) \\\frac{d}{dx}w(x) \\\frac{d^2}{dx^2}w(x) \\\frac{d^3}{dx^3}w(x) \\\end{pmatrix}= \begin{pmatrix}w(x) \\\theta(x) \\(EI)^{-1}M(x) \\(EI)^{-1}V(x) \\\end{pmatrix}$$ 여기서 $\theta(x), M(x), V(x)$ 는 각각 $x$ 에서의 기울기, 모멘트, 전단력이다.Here, $\theta(x), M(x), V(x)$ are, repectively, slope, moment, and shear force at $x$. 미분해보자.Let's differentiate. $$\frac{d}{dx}\mathbb{r}(x) = \begin{pmatrix}\frac{d}{dx}w(x) \\\frac{d^2}{dx^2}w(x) \\\frac{d^3}{dx^3}w(x) \\\frac{d^4}{dx^4}w(x) \\\end{pmatrix}=\begin{pmatrix}r_1\\r_2\\r_3\\(EI)^{-1}q(x) \\\end{pmatrix}$$ 행렬로 다시 써 보면 다음과 같다.Let's rewrite using matrices. $$\frac{d}{dx}\mathbb{r}(x) = \begin{bmatrix}0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\\0 & 0 & 0 & 0\\\end{bmatrix}\begin{pmatrix}r_0\\r_1\\r_2\\r_3 \\\end{pmatrix}+\begin{pmatrix}0\\0\\0\\(EI)^{-1} \\\end{pmatrix}q(x)$$ $x=0$ 지점에서의 초기조건을 생각해 보자.Let's think about the initial conditions at $x=0$. $$\mathbb{r}(0) = \begin{pmatrix}w(0)[m] \\\frac{d}{dx}w(0)[rad]  \\\frac{d^2}{dx^2}w(0)[m^{-1}] \\\frac{d^3}{dx^3}w(0)[m^{-2}] \\\end{pmatrix}= \begin{pmatrix}w(0)[m] \\\theta(0)[rad] \\(EI)^{-1}[N^{-1}m^{-2}]M(0)[Nm] \\(EI)^{-1}[N^{-1}m^{-2}]V(0)[N] \\\end{pmatrix}= \begin{pmatrix}0 \\0 \\(EI)^{-1}M(0) \\(EI)^{-1}V(0) \\\end{pmatrix}$$ 외팔보의 경우, $V(0)$ 는 $q(x)$ 를 $x=0$에서 $x=L$ 까지 적분한 힘과 평형을 이룰 것이다.In case of a cantilever, $V(0)$ would be in equilibrium with the integration of $q(x)$ from $x=0$ to $x=L$. $$\begin{align}V(0) &+ \int_{x=0}^{x=L}q(x)dx = 0 \\V(0) &= - \int_{x=0}^{x=L}q(x)dx\end{align}$$ $M(0)$는 $q(x)$의 도심 $\bar{x}$에 $V(0)$가 작용할 경우 $x=0$에서의 모멘트이다.$M(0)$ is the moment at $x=0$ when $V(0)$ is located at the centroid $\bar{x}$ of $q(x)$. $$M(0) = \bar{x}V(0)$$ 도심 $\bar{x}$ 는 다음과 같이 구할 수 있다.We can find the cetroid $\bar{x}$ as follows. $$\bar{x} = \frac{1}{L}\int_{x=0}^{x=L}xq(x)dx$$ 모멘트 $M(0)$는 $x$ 에서의 전단력 $V(x)$ 가 $x=0$ 점에 작용하는 모멘트와 균형을 이룬다.The moment $M(0)$ is in equilibrium with the moment by shear force $V(x)$ at $x$ on $x=0$ point. $$M(0) + \int_{x=0}^{x=L}xV(x)dx = 0 \\M(0) + \int_{x=0}^{x=L}x\left(V_0+ \int_{\xi=0}^{\xi=x}q(\xi)d\xi  \right)dx = 0$$ 우선 초기조건을 찾아보자.Let's first find the initial conditions. (Pytel & Kiusalaas, Mechanics of Materials, 2nd Ed., Cengage Learning, 2013, Example 6.1.) 분포 하중은 다음과 같다.Distributed load is as follows.<jupyter_code>w_N_m = 2000

def q_N_m(x_m):
    return w_N_m


def find_init(x_m, r):
    return np.array([-q_N_m(x_m), x_m*(q_N_m(x_m))])


L_m = 4

E_N_m2 = 200e9

I_mm4 = 113e6
I_m4 = I_mm4 * ((1e-3)**4)

EI_Nm2 = E_N_m2 * I_m4

one_over_EI_1_Nm2 = 1.0 / EI_Nm2


x_m_array = np.linspace(0, L_m, 1000+1)

x, V_M = ode_solver.rk4(find_init, x_m_array, np.array([0, 0]))

V_0_N = V_M[-1][0]
M_0_Nm = V_M[-1][1]
[V_0_N, M_0_Nm]<jupyter_output><empty_output><jupyter_text>구한 초기값을 확인해 보자.Let's check the initial values that we just found.<jupyter_code>assert abs((-w_N_m * L_m)-V_0_N) < 1e-6
assert abs((0.5 * w_N_m * L_m**2)-M_0_Nm) < 1e-6<jupyter_output><empty_output><jupyter_text>이제 처짐 곡선을 구해 보자.Let's find out the deflection curve.<jupyter_code>mat_A = np.array(
    [
        [0, 1, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1],
        [0, 0, 0, 0],
    ]
)

mat_B = np.array([0, 0, 0, one_over_EI_1_Nm2])


def dr_dx(x_m, r):
    return mat_A @ r + mat_B * q_N_m(x_m)


r0 = np.array(
    [
        0, 0, M_0_Nm*one_over_EI_1_Nm2, V_0_N*one_over_EI_1_Nm2
    ]
)


x, r = ode_solver.rk4(dr_dx, x_m_array, r0)

r[-1]


r_array = np.array(r).T
r_array.shape<jupyter_output><empty_output><jupyter_text>엄밀해와 비교Compare with exact solutions  처짐Deflection 처짐 $w(x)$의 이론해는 다음과 같다.Exact solution of the deflection $w(x)$ is as follows. $$w(x) = \frac{w_0 x^2}{24EI}\left(6L^2 -4Lx + x^2 \right)$$<jupyter_code>w_table_6_2 = (w_N_m * (x_m_array ** 2) / (24 * E_N_m2 * I_m4)) * \
    (6 * L_m ** 2 - 4 * L_m * x_m_array + x_m_array**2)
w_table_6_2[-1]


py.plot(x, r_array[0, :], '.', label='numerical')
py.plot(x_m_array, w_table_6_2, label='exact')

py.xlabel('x[m]')
py.ylabel('w(x)[m]')

py.legend(loc=0)

py.grid(True)<jupyter_output><empty_output><jupyter_text>상대오차Relative error<jupyter_code>w_numerical_m = np.interp(x_m_array, x, r_array[0, :])
w_error = nl.norm(w_table_6_2 - w_numerical_m) / nl.norm(w_table_6_2)

assert 1e-7 > w_error, f"deflection error {w_error*100}% larger than expected"<jupyter_output><empty_output><jupyter_text>전단력Shear force 전단력의 이론해Shear force (in theory)<jupyter_code>sf_N = V_0_N + w_N_m * x_m_array<jupyter_output><empty_output><jupyter_text>전단력의 수치해Numerical solution of the shear force<jupyter_code>sf_numerical_N = np.interp(x_m_array, x, r_array[3, :]) * EI_Nm2


py.plot(x_m_array, sf_numerical_N, '.', label='numerical')
py.plot(x_m_array, sf_N, label='exact')
py.xlabel('x[m]')
py.ylabel('V[N]')

py.legend(loc=0)

py.grid(True)


sf_error = nl.norm(sf_N - sf_numerical_N) / nl.norm(sf_N)

assert 1e-7 > sf_error, f"shear force error ({sf_error}) larger than expected"<jupyter_output><empty_output><jupyter_text>굽힘모멘트Bending moment 굽힘모멘트의 이론해Bending Moment (in theory)<jupyter_code>bm_Nm = M_0_Nm + V_0_N * x_m_array + 0.5 * w_N_m * x_m_array ** 2<jupyter_output><empty_output><jupyter_text>굽힘모멘트의 수치해Numerical solution of the bending moment<jupyter_code>bm_numerical_Nm = np.interp(x_m_array, x, r_array[2, :]) * EI_Nm2


py.plot(x_m_array, bm_numerical_Nm, '.', label='numerical')
py.plot(x_m_array, bm_Nm, label='exact')

py.xlabel('x[m]')
py.ylabel('M[Nm]')

py.legend(loc=0)

py.grid(True)


bm_error = nl.norm(bm_Nm - bm_numerical_Nm) / nl.norm(bm_Nm)

assert 1e-7 > bm_error, f"bending moment error {bm_error} larger than expected"<jupyter_output><empty_output><jupyter_text>기울기Slope 기울기의 엄밀해Exact solution of the slope<jupyter_code>theta_table_6_2 = (w_N_m  / (24 * E_N_m2 * I_m4)) * \
    (12 * L_m ** 2 * x_m_array - 12 * L_m * x_m_array ** 2 + 4 * x_m_array ** 3)
theta_table_6_2[-1]<jupyter_output><empty_output><jupyter_text>기울기의 수치해Numerical solution of the slope<jupyter_code>slope_numerical_radian = np.interp(x_m_array, x, r_array[1, :])


py.plot(x_m_array, py.rad2deg(slope_numerical_radian), '.', label='numerical')
py.plot(x_m_array, py.rad2deg(theta_table_6_2), label='exact')

py.xlabel('x[m]')
py.ylabel('$\\theta$(deg)')

py.grid(True)<jupyter_output><empty_output><jupyter_text>상대 오차 확인Check the relative error of the slope<jupyter_code>theta_error = nl.norm(theta_table_6_2 - slope_numerical_radian) / nl.norm(theta_table_6_2)

assert 1e-7 > bm_error, f"slope moment error {theta_error} larger than expected"<jupyter_output><empty_output><jupyter_text>ref : <NAME>, Numerical Methods for Engineers, 2018, [Online] Available : https://www.ntnu.no/wiki/download/attachments/69242130/main.pdf  연습 문제Exercises 외팔보 상의 분포하중 $q(x)$가 다음과 같을 때 질문에 답하시오:When the distributed load is as follows, answer the question:$$q(x) = 2\left(1 - \frac{1}{L}x \right)$$ 따로 주어지지 않은 값은 위의 예의 값을 사용할 수 있음.You may use values above if not given. 도전 과제 1: $x=0$ 에서의 전단력의 초기 조건을 구하시오.Try This 1: Find the initial condition of the shear force at $x=0$. 도전 과제 2: $x=0$ 에서의 굽힘모멘트의 초기 조건을 구하시오.Try This 2: Find the initial condition of the bending moment at $x=0$. 도전 과제 3: $x$ 에서의 처짐 곡선을 구하시오.Try This 3: Find the deflection curve at $x$. 도전 과제 4: 이론해와 비교해 보시오.Try This 4: Compare with the exact solution.$$w(x)=\frac{w_0 x^2}{120 L \cdot EI}\left( 10L^3 - 10L^2x + 5Lx^2 - x^3 \right)$$  Final Bell마지막 종<jupyter_code># stackoverfow.com/a/24634221
import os
os.system("printf '\a'");<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup