import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

Exercise potential energy statically indeterminate structure#

Find the displacement in C using potential energy and compare it to the result from solving the differential equation.

../_images/structure7.svg

Fig. 1 Structure#

Solve using potential energy#

Let’s start by defining the trial function. Let’s assume a parabolic displacement field with the shape \(u\left(x\right) = a + b x + c x^2\).

Exercise (Boundary conditions kinematical admissible trial function)

The displacement field must be kinematically admissible. Which boundary conditions must be enforced for that?

Exercise (Kinematical admissibility trial function)

Now, let’s use SymPy to rewrite the trial function so that it is kinematically admissible.

import sympy as sym
a, b, c, x, L, EA, F = sym.symbols('a, b, c, x, L, EA, F')
u_trial_general = a + b * x + c * x **2

eq1 = sym.Eq(u_trial_general.subs(x, #value to substitute for x), #right hand side of equation)
eq2 = #...
sol = sym.solve((eq1,eq2),(a,b))
display(sol)
u_trial = u_trial_general.subs(sol)
display(u_trial)

Let’s check your solution graphically for a value of \(1\) for \(c\) and \(5\) for \(L\):

import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_formats = ['svg']
%matplotlib inline
u_trial_func = sym.lambdify(x, u_trial.subs(sol).subs(c,1).subs(L,5))
x_plot = np.linspace(0, 10,100)
fig, ax = plt.subplots()
ax.plot(u_trial_func(x_plot),x_plot)
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('u(x)')
ax.set_ylabel('x')
fig.patch.set_alpha(0.0)
ax.patch.set_alpha(0.0)

Exercise (epsilon)

For the strain energy we needed the strain distribution instead of the displacement. The strain is the derivative of the displacement, which can be found with SymPy:

For this exercise you just have to run the cell below.

epsilon_trial = sym.diff(u_trial,x)
display(epsilon_trial)

Exercise (strain energy)

Now let’s evaluate the strain energy:

E_v_trial = sym.integrate(#term to integrate,(x,#integrate from, #integrate to)) + ...
display(E_v_trial)

Exercise (Work by force)

What’s the work done by the force \(F\)?

Exercise (Potential energy)

Find the relation for the potential energy.

V = #expression for V
display(V)

Exercise (Stationary strain energy)

Now, let’s impose the stationary condition for the strain energy in an equation which we can use to solve for the unknown \(c\).

eq3 = sym.Eq(#left hand side of equation, #right hand side of equation)
display(eq3)
c_sol = sym.solve(eq3,c)[0]
display(c_sol)
u_trial_sol = u_trial.subs(c,c_sol)
display(u_trial_sol)

Let’s check your solution graphically for a value of \(5\) m for \(L\), \(2000\) kN for \(EA\) and \(30\) kN for \(F\):

%matplotlib inline
L_test = 5
EA_test = 2000
F_test = 30
u_trial_func = sym.lambdify(x, u_trial_sol.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
x_plot = np.linspace(0, L_test*2,100)
fig, ax = plt.subplots()
ax.plot(u_trial_func(x_plot),x_plot)
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('u(x)')
ax.set_ylabel('x')
fig.patch.set_alpha(0.0)
ax.patch.set_alpha(0.0)

Exercise (Normal force)

What is the normal force distribution?

epsilon = #...

N_trial_sol_AC = #
N_trial_sol_BC = #

display(N_trial_sol_AC.simplify())
display(N_trial_sol_BC.simplify())

Let’s check your solution graphically for a value of \(5\) m for \(L\), \(2000\) kN for \(EA\) and \(30\) kN for \(F\):

%matplotlib inline
L_test = 5
EA_test = 2000
F_test = 30
fig, ax = plt.subplots()

N_trial_func_AC = sym.lambdify(x, N_trial_sol_AC.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
N_trial_func_BC = sym.lambdify(x, N_trial_sol_BC.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
x_plot_AC = np.linspace(0, L_test,2)
x_plot_BC = np.linspace(L_test, L_test*2,2)
ax.plot(N_trial_func_AC(x_plot_AC),x_plot_AC)
ax.plot(N_trial_func_BC(x_plot_BC),x_plot_BC)
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('N(x)')
ax.set_ylabel('x')
fig.patch.set_alpha(0.0)
ax.patch.set_alpha(0.0)

Solve with differential equation#

Now, let’s solve the same structure using the differential equation.

Exercise (Generalised relations)

Given a load \(q=0\), integrate twice to find the normal force, strain and displacement relations with integration constants:

C_1, C_2, C_3, C_4 = sym.symbols('C_1, C_2, C_3, C_4')

q_AC = 0
q_BC = 0

N_AC = C_1 + sym.integrate(#..., x)
N_BC = C_3 + #...

display(N_AC)
display(N_BC)

epsilon_AC = #...
epsilon_BC = #...

display(epsilon_AC)
display(epsilon_BC)

u_AC = C_2 + #...
u_BC = C_3 + #...

display(u_AC)
display(u_BC)

Exercise (Boundary conditions)

What are the boundary conditions?

Exercise (Solve boundary conditions)

Now, let’s use SymPy to solve for the boundary conditions using sym.solve similarly as you did in the exercise for the kinematical admissibility trial function and find the solution for \(u_\text{AC}\), \(u_\text{BC}, \)N_\text{AC}\( and \)N_\text{BC}

eq4 = #...
#...

C_sol = sym.solve((eq4,#...),[#])
display(sol)

u_AC_sol = u_AC.subs(C_sol)
u_BC_sol = u_BC.subs(C_sol)
display(u_AC_sol)
display(u_BC_sol)
N_AC_sol = #
N_BC_sol = #
display(N_AC_sol)
display(N_BC_sol)

Let’s check your solution graphically for a value of \(5\) m for \(L\), \(2000\) kN for \(EA\) and \(30\) kN for \(F\):

%matplotlib inline
u_AC_func = sym.lambdify(x, u_AC_sol.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
u_BC_func = sym.lambdify(x, u_BC_sol.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
x_plot_AC = np.linspace(0, L_test,100)
x_plot_BC = np.linspace(L_test, L_test*2,100)
fig, ax = plt.subplots()
ax.plot(u_AC_func(x_plot_AC),x_plot_AC)
ax.plot(u_BC_func(x_plot_BC),x_plot_BC)
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('u(x)')
ax.set_ylabel('x')
fig.patch.set_alpha(0.0)
ax.patch.set_alpha(0.0)
N_AC_func = sym.lambdify(x, N_AC_sol.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
N_BC_func = sym.lambdify(x, N_BC_sol.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
x_plot_AC = np.linspace(0, L_test,2)
x_plot_BC = np.linspace(L_test, L_test*2,2)
fig, ax = plt.subplots()
ax.plot(N_AC_func(x_plot_AC)*np.ones(2),x_plot_AC)
ax.plot(N_BC_func(x_plot_BC)*np.ones(2),x_plot_BC)
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('u(x)')
ax.set_ylabel('x')
fig.patch.set_alpha(0.0)
ax.patch.set_alpha(0.0)

Comparison#

Now, let’s compare the plots of both solution strategies, you can try different values:

from ipywidgets import widgets, interact
%matplotlib widget

fig, ax = plt.subplots()
fig.canvas.toolbar_visible = False
def func(L_test, EA_test, F_test):
    ax.clear()
    u_trial_func = sym.lambdify(x, u_trial_sol.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    u_AC_func = sym.lambdify(x, u_AC_sol.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    u_BC_func = sym.lambdify(x, u_BC_sol.subs(sol).subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    x_plot = np.linspace(0, L_test*2,100)
    x_plot_AC = np.linspace(0, L_test,2)
    x_plot_BC = np.linspace(L_test, L_test*2,2)
    ax.plot(u_AC_func(x_plot_AC),x_plot_AC,label='u_AC ODE')
    ax.plot(u_BC_func(x_plot_BC),x_plot_BC, label='u_BC ODE')
    ax.plot(u_trial_func(x_plot),x_plot, label='u trial')
    ax.spines['right'].set_color('none')
    ax.spines['top'].set_color('none')
    ax.spines['bottom'].set_position('zero')
    ax.spines['left'].set_position('zero')
    ax.set_xlabel('u(x)')
    ax.set_ylabel('x')
    ax.legend(loc='best')
    fig.patch.set_alpha(0.0)
    ax.patch.set_alpha(0.0)

interact(func, L_test = widgets.FloatSlider(min=0.1, max=5, value=2, step=0.1, description="Length (m)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}),
         EA_test = widgets.FloatSlider(min=1000, max=10000, value=2000, step=100, description="EA (kN)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}),
         F_test = widgets.FloatSlider(min=10, max=100, value=30, step=1, description="Force (kN)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}) );
%matplotlib widget
fig, ax = plt.subplots()
fig.canvas.toolbar_visible = False

def func(L_test, EA_test, F_test):
    ax.clear()
    N_trial_func_AC = sym.lambdify(x, N_trial_sol_AC.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    N_trial_func_BC = sym.lambdify(x, N_trial_sol_BC.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    N_AC_func = sym.lambdify(x, N_AC_sol.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    N_BC_func = sym.lambdify(x, N_BC_sol.subs(L,L_test).subs(EA,EA_test).subs(F,F_test))
    x_plot = np.linspace(0, L_test*2,2)
    x_plot_AC = np.linspace(0, L_test,2)
    x_plot_BC = np.linspace(L_test, L_test*2,2)
    ax.plot(N_AC_func(x_plot_AC)*np.ones(2),x_plot_AC,label='u_AC ODE')
    ax.plot(N_BC_func(x_plot_BC)*np.ones(2),x_plot_BC, label='u_BC ODE')
    ax.plot(N_trial_func_AC(x_plot_AC),x_plot_AC, label='u_AC trial')
    ax.plot(N_trial_func_BC(x_plot_BC),x_plot_BC, label='u_BC trial')
    ax.spines['right'].set_color('none')
    ax.spines['top'].set_color('none')
    ax.spines['bottom'].set_position('zero')
    ax.spines['left'].set_position('zero')
    ax.set_xlabel('N(x)')
    ax.set_ylabel('x')
    ax.legend(loc='best')
    fig.patch.set_alpha(0.0)
    ax.patch.set_alpha(0.0)

interact(func, L_test = widgets.FloatSlider(min=0.1, max=5, value=2, step=0.1, description="Length (m)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}),
         EA_test = widgets.FloatSlider(min=1000, max=10000, value=2000, step=100, description="EA (kN)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}),
         F_test = widgets.FloatSlider(min=10, max=100, value=30, step=1, description="Force (kN)",readout_format='.1f',style= {'description_width': '200px'},layout = {'width': '450px'}) );

Exercise (Comparison)

Would you consider the trial function solution for proper estimates for displacements? How about normal forces?