У меня серьезные проблемы с решением перевода трех связанных уравнений в python.
Эти 3 уравнения происходят из DE 4-го порядка, используемого для вычисления изгибающего момента подводного трубопровода, который был заменен / преобразован вУравнения 1-го порядка и связаны друг с другом.
Латексные уравнения
W1, w2 и w3 - это три секции конвейера, которые все получены из интегрирования первогоуравнение.w1 начинается в x = 0 и заканчивается в оценочном x = I1, w2 начинается в x = I1 и заканчивается в x = I1 + I2, а w3 заканчивается в I1 + I2 + I3.
Система имеет 15 неизвестныхограничения и 15 граничных условий условия, для которых я хочу решить.
Граничные условия
- w_1 (0) = 0 (прогиб при x_1 = 0 равен 0), D = 0
- w_1 '(0) = 0 (наклон при x_1 = 0 равен 0), C = 0
- w_1'' (0) = 0 (изгибающий момент при x_1 = 0 равен 0), B = 0
- w_1 (I_1) = -p
- w_2 (0)= p, H = p
- w_1 '(I_1) = w_2' (0)
- w_1 '' (I_1) = w_2 '' (0)
- w_2 (I_2) = 0
- w_3 (0) = 0
- w_2 '(I_2) = w_3' (0)
- w_2 '' (I_2) = w_3 '' (0)
- w_2 '' '(I_2) = w_3' '' (0)
- w_3 (I_3) = d
- w_3 (I_3) '= 0
- w_3 (I_3)' '= 0
Я посмотрел онлайн и попробовал некоторые коды самостоятельно, но до сих пор, к сожалению, я не получил никаких полезных результатов,Пример моего кода приведен ниже, где для неизвестных переменных даны 15 начальных / предположительных значений.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import scipy.optimize as opt
#-------------- Pipe and hook plough parameters ------------#
#
#-------------- PIPE ------------#
#
p = -0.4 # HOOK ELEVATION [m]
d = 1 # BURIAL DEPTH PIPELINE [m]
D_o = 0.324 # OUTSIDE DIAMETER PIPELINE [m]
WT = 0.0254 # WALL THICKNESS PIPELINE [m]
D_i = D_o - 2*WT # INSIDE DIAMETER PIPELINE [m]
A = 0.25*np.pi*(D_o**2 - D_i**2) # CROSS SECTIONAL AREA [m]
E_e = 2.1*10**11 # YOUNGS MODULUS
sig_max = 413*(10**6) # SYMS
gamma = 1.5 # SAFETY FACTOR Note: 0.67 SYMS
#
#-------------- ENIVIRONMENT ------------#
#
g = 9.81 # GRAVITATIONAL ACCELERATION [m/s2]
rho_s = 7840 # STEEL DENSITY [kg/m3]
rho_w = 1025 # WATER DENSITY [kg/m3]
gamma_wet = 18500 # SUBMERGED UNIT WEIGHT SAND [N/m3] Note: North sea soil, fine to medium sand.
gamma_water = 10250 # SPECIFIC WEIGHT WATER [N/m3]
phi_g = 37 # HOOK ANGLE TO NEUTRAL [deg] Note: tan(37) = 0.754
K_0 = 0.8 # LATERAL EARTH PRSSURE COEFFICIENT Note: can be between 0.5 an 0.8
#
#-----------------------------------------------------------#
#-------------- Initialisation calculations ------------#
phi_deg= np.tan((phi_g/360)*np.pi*2) #
q_subm = (rho_s*A*g) - (rho_w*np.pi*(D_o/2)**2)*g # SUBMERGED WEIGHT [N/m]
q_grnd = D_o*(gamma_wet-gamma_water)*d # SOIL WEIGHT [N/m]
q_fric = K_0*phi_deg*(gamma_wet-gamma_water)*(d**2) # SOIL FRICTION [N/m]
Q_grnd = q_grnd/d # [N/m]
Q_fric = q_fric/(d**2) # [N/m]
W = (np.pi*(D_o**4 - D_i**4))/(32*D_o) # SECTION MODULUS [m3]
I = (np.pi*(D_o**4 - D_i**4))/(64) # MOMENT OF INERTIA [m3]
EI = E_e*I # [Nm2]
M_max = (sig_max*W)/gamma # MAXIMUM ALLAWABLE ELASTIC BENDING MOMENT [Nm]
R_min = EI/M_max # MINIMUM ALLAWABLE ELASTIC RADIUS OF CURVATURE [m]
beta = np.sqrt(np.sqrt((Q_fric+Q_grnd)/(4*EI))) # SUBSTITUTION FOR IN DIFF EQ. [1/m]
#
#-----------------------------------------------------------#
#Test if initialisation calculations are correct
print('q_subm',q_subm,'\n','q_grnd',q_grnd,'\n','q_fric',q_fric,'\n','Q_grnd',Q_grnd,'\n','Q_fric',Q_fric,'\n','W',W,'\n','I',I,'\n','EI',EI,'\n','M_max',M_max,'\n','R_min',R_min,'\n','beta',beta)
def main():
nobservations = 20
# Guess values for the unkown variables
I1 = 46.483
I2 = 5.916
I3 = 21.90
A = -3.828*10**(-5)
B = 0
C = 0
D = 0
E = -1.928*10**(-4)
F = 4.270*10**(-3)
G = 0.049
H = p
K = 0.052
L = 0.1
M = -0.021
N = 0.787
f, x1, x2, x3 = generate_data(nobservations, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N)
print ('Non-linear results (should be {}, {}, {}):'.format(I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N))
print (nonlinear_invert(f, x1, x2, x3))
def generate_data(nobservations, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N, noise_level=0.01):
x, y, z = np.random.random((3, nobservations))
noise = noise_level * np.random.normal(0, noise_level, nobservations)
f = func(x, y, z, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N) + noise
return f, x, y, z
def func(x1, x2, x3, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N):
f = (A*x1**3 + B*x1**2 + C*x1 + D + (q_subm*x1**4)/(24*EI))
#I would like to include/couple these two equations
f2 = (E*x2**3 + F*x2**2 + G*x2 + H + (q_subm*x2**4)/(24*EI))
f3 = K*np.exp(np.sqrt(2)*beta*x3)+L*np.exp(-np.sqrt(2)*beta*x3)+M*np.cos(np.sqrt(2)*beta*x3)+N*np.sin(np.sqrt(2)*beta*x3)-(q_subm/(Q_grnd+Q_fric))
return f
def nonlinear_invert(f, x1, x2, x3):
# "curve_fit" expects the function to take a slightly different form...
def wrapped_func(observation_points, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N):
x1, x2, x3 = observation_points
return func(x1, x2, x3, I1, I2, I3, A, B, C, D, E, F, G, H, K, L, M, N)
xdata = np.vstack([x1, x2, x3])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
main()
Кто-нибудь может быть так любезен, чтобы помочь мне?