Python - memecahkan persamaan balok Bernoulli dengan scipy

Nov 26 2020

Proses menjawab pertanyaan sudah dimulai pada pertanyaan pada link di bawah, tetapi topik tersebut secara khusus tentang mengintegrasikan suatu fungsi, yang telah dijawab. Jadi saya menambahkan pertanyaan baru.

Python - Mengintegrasikan fungsi dan hasil plot

MASALAH: bagaimana menyelesaikan persamaan balok y '' (x) = M (x) / (E * I) menggunakan integrasi scipy.

SOLUSI, milik 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()

SOLUSI TEPAT: solusi eksak dibuat dengan mengintegrasikan persamaan balok dua kali menggunakan integral pasti dan menggunakan kondisi batas untuk menentukan konstanta integral. Semuanya dijelaskan di tautan wiki di atas. Di bawah ini adalah kode untuk memplot y '' (x), y '(x) (kemiringan) dan y (x) (defleksi). Diagram diputar, ujung bebas balok berada pada 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()

SOLUSI APPROXIMATE (JENIS): kode di bawah ini dibuat oleh willcrack. Bentuknya terlihat lebih baik dari pada pertanyaan sebelumnya tetapi nilainya masih kurang baik.

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()

Jawaban

2 gboffi Nov 26 2020 at 09:45

Anda sedang mengintegrasikan persamaan diferensial, pendekatan komputasi Anda dalam satu loop adalah integral tertentu, katakanlah, sub-optimal.

Pendekatan standar di Scipy adalah penggunaan scipy.integrate.solve_ivp, yang menggunakan metode integrasi yang sesuai (secara default, Runge-Kutta 45) untuk memberikan solusi dalam hal objek khusus.

Seperti biasa di bidang integrasi numerik dari persamaan diferensial biasa, metode ini terbatas pada sistem persamaan diferensial orde-1, tetapi persamaan derajat ke-2 Anda dapat diubah menjadi sistem persamaan derajat ke-1 yang memperkenalkan fungsi pembantu