Как правильно решить одномерное волновое уравнение, чтобы получить профиль смещения (периодическая проблема граничных условий)? - PullRequest
Новая
       22

Как правильно решить одномерное волновое уравнение, чтобы получить профиль смещения (периодическая проблема граничных условий)?

0 голосов
/ 14 февраля 2019

Я пытаюсь решить 1D волновое уравнение для кучи с периодической БК (периодической нагрузкой).

Я почти уверен в своих формулах дискретизации.Единственное, в чем я не уверен, так это периодические BC и время ( t ) там ==> sin(omega*t).

Когда я настраиваю его так, как сейчас, это дает мне странный профиль смещения.Однако, если я установлю его на sin(omega*1) или sin(omega*2), ... и т. Д., Это будет похоже на синусоидальную волну, но в основном это означает, что sin(omega*t), то есть sin(2*pi*f*t), равно 0, когда t является целочисленным значением.

Я кодирую все в 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

Документ My Jupyter Notebook .

...