Python - Résolution de l'équation de faisceau de Bernoulli avec Scipy

Nov 26 2020

Le processus de réponse à la question a déjà commencé dans la question sur le lien ci-dessous, mais ce sujet concernait spécifiquement l'intégration d'une fonction, à laquelle il a été répondu. J'ai donc ajouté une nouvelle question.

Python - Intégrer une fonction et tracer les résultats

LE PROBLÈME: comment résoudre une équation de poutre y '' (x) = M (x) / (E * I) en utilisant scipy integrer.

SOLUTION, gracieuseté de gboffi:

#---------- DESCRIPTION

# cantilever beam with point load P at the free end
# original beam equation: y''(x) = M(x)/(E*I)
# moment equation: M(x) = -P*x
# x goes from the free end to the clamped end

# we have a second order diff eq: y''(x) = x
# we implement a new function:
#      h = y',
#      h' = y'' = M(x) = x

# we get a system of two ODE of first order
#      y' = h
#      h' = x

# we write the equations in vector form
#     Y' = F(x, Y(x)) = F(x,Y)

# we define a function that returns the original values

#----------- CODE

from __future__ import division
from numpy import linspace
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Exact solution, E*Iy = const, y1 = y', y0 = y, 
w = 10  #beam cross sec width (mm)
h = 10  #beam cross sec height (mm)
Iy = (w*h**3)/12   #cross sec moment of inertia (mm^4)
E = 200000   #steel elast modul (N/mm^2)
L = 100  #beam length(mm)
P = 100   #point load (N)

x = linspace(0, L, 51)

y1 = (-P/(2*E*Iy))*x**2+(P*L**2)/(2*E*Iy)
y0 = (-P/(6*E*Iy))*x**3+((P*L**2)/(2*E*Iy))*x-(2*P*L**3)/(6*E*Iy)

# Define the vector function for E=const for integration
def F(x,Y):
    #unpack the vector function
    y = Y[0]
    h = Y[1]
    #compute the derivatives
    dy_dx = h
    dh_dx = (-P/(E*Iy))*x
    #return the vector of derivatives values
    return [dy_dx, dh_dx]

# Numerical solution
s = solve_ivp(
    F, # Y[0]=y0, Y[1]=y1, dy0dx=y1, dy1dx=x
    [L, 0.0], # interval of integration (NB: reversed, because...)
    [0.0, 0.0], # initial conditions (at the 1st point of integ interval)
    t_eval=linspace(L, 0, 101) # where we want the solution to be known
    )

# Plotting
fig, (ax1, ax2) = plt.subplots(2)
ax1.plot(x, y0, label="Exact y")
ax2.plot(x, y1, label="Exact y'")
ax1.plot(s.t[::2], s.y[0][::2], label="Numeric y",  linestyle='', marker='.')
ax2.plot(s.t[::2], s.y[1][::2], label="Numeric y'", linestyle='', marker='.')
plt.show()

SOLUTION EXACTE: la solution exacte est faite en intégrant l'équation de poutre deux fois en utilisant des intégrales définies et en utilisant les conditions aux limites pour définir les constantes intégrales. Tout est expliqué dans le lien wiki ci-dessus. Voici le code pour tracer les y '' (x), y '(x) (pente) et y (x) (déflexion). Le diagramme est retourné, l'extrémité libre de la poutre est à x = 0.

from __future__ import division  #to enable normal floating division
import numpy as np
import matplotlib.pyplot as plt

# Beam parameters
w = 10  #beam cross sec width (mm)
h = 10  #beam cross sec height (mm)
I = (w*h**3)/12   #cross sec moment of inertia (mm^4)
I1 = (w*h**3)/12
E = 200000   #steel elast modul (N/mm^2)
L = 100  #beam length(mm)
F = 100   #force (N)

# Define equations
def d2y_dx2(x):
    return (-F*x)/(E*I)

def dy_dx(x):
    return (1/(E*I))*(-0.5*F*x**2 + 0.5*F*L**2)

def y(x):
    return (1/(E*I))*(-(1/6)*F*(x**3) + (1/2)*F*(L**2)*x - (1/3)*F*(L**3))

# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)

a = 0
b = L
x = np.linspace(a,b,100)

ax1.plot(x, d2y_dx2(x))
ax2.plot(x, dy_dx(x))
ax3.plot(x, y(x))
plt.show()

SOLUTION APPROXIMATIVE (TYPE DE): le code ci-dessous a été créé par willcrack. La forme est meilleure que dans la question précédente mais les valeurs ne sont toujours pas correctes.

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100

# Integration parameters
a = 0.0
b = L

# Define the beam equation
def d2y_dx2(x,y=None):
    return (-F*x)/(E*I)

    
# Define the integration1 - slope
def slope(x):
    slope_res = np.zeros_like(x)
    for i,val in enumerate(x):
        y,err = integrate.quad(f,a,val)
        slope_res[i]=y
    return slope_res

# Define the integration1 - deflection
def defl(x):
    
    defl_res = np.zeros_like(x)
    for i,val in enumerate(x):
        y, err = integrate.dblquad(d2y_dx2,0,val, lambda x: 0, lambda x: val)
        defl_res[i]=y
    return defl_res

# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()

Réponses

2 gboffi Nov 26 2020 at 09:45

Vous intégrez une équation différentielle, votre approche du calcul en boucle des intégrales définies est, disons, sous-optimale.

L'approche standard de Scipy est l'utilisation de scipy.integrate.solve_ivp, qui utilise une méthode d'intégration appropriée (par défaut, Runge-Kutta 45) pour fournir la solution en termes d'objet spécial.

Comme d'habitude dans le domaine de l'intégration numérique d'équations différentielles ordinaires, la méthode se limite à un système d'équation différentielle du 1er ordre, mais votre équation du 2ème degré peut être transformée en un système d'équations du 1er degré introduisant une fonction d'assistance