0

I'm trying to solve a 1D wave equation for the pile with periodic BC (periodic load).

I'm pretty sure about my discretization formulas. The only thing I'm not sure about is the periodic BC and time (t) in there ==> sin(omega*t).

When I set it up as it is right now, it's giving me a weird displacement profile. However, if I set it up to be sin(omega*1) or sin(omega*2),... etc, it resembles a sine wave, but it basically means that sin(omega*t), i.e. sin(2*pi*f*t), is equal to 0 when t is an integer value..

I code everything up in Jupyter Notebook together with the visualization part, but the solution is nowhere near the propagating sine wave.

Here is the relevant Python code:

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 document.

trincot
  • 317,000
  • 35
  • 244
  • 286
Tomasz Plewka
  • 1
  • 4
    Welcome to stackoverflow. Best practice is to include the relevant code in your question. See https://stackoverflow.com/help/how-to-ask – Simon.S.A. Feb 14 '19 at 03:54
  • I have extracted the relevant Python code from your document, and embedded it inside the question. As commented above, it is important that the question is self-contained: make sure the problem description is complete without requiring readers to follow hyperlinks. – trincot Feb 14 '19 at 19:21
  • Which values of `f` are you testing with? If that's an integer then `sin(2*pi*f*t)` will return 0 for any integer `t`. – Hoog Feb 14 '19 at 19:27
  • trincot, thank you for doing that. – Tomasz Plewka Feb 15 '19 at 00:40
  • I'm testing it with f = 5, 15, 30, etc. so it's gonna be 0 each time when I do that. However, from what I understand, time (t) should be included as dt (time step) which is rather small number ==> e.g. for dz = 0.25, dt = 0.0001 in this case, so it'll produce non-zero value. The question I keep asking myself is this correct way of thinking and if so, what I am doing wrong with the code (or discretization) that final displacement is not sine-wave-looking one – Tomasz Plewka Feb 15 '19 at 00:45

0 Answers0