0

我正在尝试为具有周期性 BC（周期性负载）的桩求解一维波动方程。

我很确定我的离散化公式。我唯一不确定的是那里的周期性 BC 和时间（t） ==> sin(omega*t)。

当我现在设置它时，它给了我一个奇怪的位移曲线。但是，如果我将其设置为sin(omega*1)or sin(omega*2),... 等，它类似于正弦波，但它基本上意味着sin(omega*t), 即t为整数值sin(2*pi*f*t)时等于 0 ..

我在 Jupyter Notebook 中将所有内容与可视化部分一起编写，但解决方案与传播正弦波相去甚远。

以下是相关的 Python 代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

def oned_wave(L, dz, T, p0, Ep, ro, f):

    """Solve u_tt=c^2*u_xx on (0,L)x(0,T]"""
    """Assume C = 1"""

    p = p0
    E = Ep
    ro = ro
    omega = 2*np.pi*f
    c = np.sqrt(E/ro)
    C = 1                                     # Courant number

    Nz = int(round(L/dz))
    z = np.linspace(0, L, Nz+1)               # Mesh points in space
    dt = dz/c                                 # Time step based on Courant Number
    Nt = int(round(T/dt))
    t = np.linspace(0, Nt*dt, Nt+1)           # Mesh points in time
    C2 = C**2                                 # Help variable in the scheme

    # Make sure dz and dt are compatible with z and t
    dz = z[1] - z[0]
    dt = t[1] - t[0]

    w = np.zeros(Nz+1)                        # Solution array at new time level
    w_n = np.zeros(Nz+1)                      # Solution at 1 time level back
    w_nm1 = np.zeros(Nz+1)                    # Solution at 2 time levels back

    # Set initial condition into w_n
    for i in range(0,Nz+1):
        w_n[i] = 0

    result_matrix = w_n[:]                    # Solution matrix where each row is displacement at given time step

    # Special formula for first time step
    for i in range(1, Nz):
        w[i] = 0.5*C2 * w_n[i-1] + (1 - C2) * w_n[i] + 0.5*C2 * w_n[i+1]

    # Set BC
    w[0] = (1 - C2) * w_n[i] + C2 * w_n[i+1] - C2*dz*((p*np.sin(omega*dt))/E) # this is where, I think, the mistake is: sin(omega*t)
    w[Nz] = 0

    result_matrix = np.vstack((result_matrix, w)) # Append a row to the solution matrix

    w_nm1[:] = w_n; w_n[:] = w                # Switch variables before next step

    for n in range(1, Nt):
        # Update all inner points at time t[n+1]
        for i in range(1, Nz):
            w[i] = - w_nm1[i] + C2 * w_n[i-1] + 2*(1 - C2) * w_n[i] + C2 * w_n[i+1]

        # Set BC
        w[0] = - w_nm1[i] + 2*(1 - C2) * w_n[i] + 2*C2 * w_n[i+1] - 2*dz*((p*np.sin(omega*(dt*n)))/E) # this is where, I think, the mistake is: sin(omega*t)
        w[Nz] = 0

        result_matrix = np.vstack((result_matrix, w)) # Append a row to the solution matrix

        w_nm1[:] = w_n; w_n[:] = w            # Switch variables before next step

    return result_matrix

我的 Jupyter Notebook 文档。

4

0 回答 0