diff options
| -rw-r--r-- | src/main.py | 45 |
1 files changed, 32 insertions, 13 deletions
diff --git a/src/main.py b/src/main.py index 5e0d682..fffb570 100644 --- a/src/main.py +++ b/src/main.py @@ -8,7 +8,7 @@ figure, axes = plt.subplots() cb = None -Re = 100 +Re = 400 nu = 1 / Re domain_size = (1, 2) step = 0.05 @@ -32,9 +32,7 @@ v = np.zeros(shape=shape, dtype=float) v_star = np.zeros(shape=shape, dtype=float) v_new = np.zeros(shape=shape, dtype=float) -p = np.zeros(shape=shape, dtype=float) -p_star = np.zeros(shape=shape, dtype=float) -p_new = np.zeros(shape=shape, dtype=float) +p = np.random.rand(shape[0], shape[1]) d_e = np.zeros(shape=shape, dtype=float) d_n = np.zeros(shape=shape, dtype=float) @@ -46,19 +44,21 @@ def u_boundary(t, y): return 1 if t_m < t else 0.5 * (np.sin(0.5 * PI * (2 * t / t_m - 1)) + 1) return 6 * f(t) * (y - h_c) * (1 - y) / (1 - h_c)**2 +def assert_rule2(value): + assert(value > -0.01) # > 0 # Loop error = 1 precision = 10 ** -7 -t = 0.35 +t = 0.4 iteration = 0 while error > precision: iteration += 1 # Inflow boundary condition for i in range(N): - y = domain_size[0] - i * step + y = i * step u_star[i][0] = u_boundary(t, y) v_star[i][0] = 0 @@ -69,6 +69,13 @@ while error > precision: u_star[0][j] = 0 u_star[N - 1][j] = 0 + # Backwards-facing step boundary conditions (same as sides) + for i in range(int(h_c / step)): + for j in range(int(l_c / step)): + u_star[i][j] = 0 + v_star[i][j] = 0 + p[i][j] = 0 + # x-momentum for i in range(1, N - 1): for j in range(1, M - 1): @@ -81,6 +88,10 @@ while error > precision: a_W = +0.5 * u_W * step + nu a_N = -0.5 * v_N * step + nu a_S = +0.5 * v_S * step + nu + assert_rule2(a_E) + assert_rule2(a_W) + assert_rule2(a_N) + assert_rule2(a_S) a_e = 0.5 * step * (u_E - u_W + v_N - v_S) + 4 * nu A_e = -step @@ -102,6 +113,10 @@ while error > precision: a_W = +0.5 * u_W * step + nu a_N = -0.5 * v_N * step + nu a_S = +0.5 * v_S * step + nu + assert_rule2(a_E) + assert_rule2(a_W) + assert_rule2(a_N) + assert_rule2(a_S) a_n = 0.5 * step * (u_E - u_W + v_N - v_S) + 4 * nu A_n = -step @@ -111,13 +126,6 @@ while error > precision: v_star[i][j] = (a_E * v[i][j + 1] + a_W * v[i][j - 1] + a_N * v[i - 1][j] + a_S * v[i + 1][j]) / a_n + d_n[i][j] * (p[i][j] - p[i + 1][j]) - # Backwards-facing step boundary conditions (same as sides) - for i in range(int(h_c / step)): - for j in range(int(l_c / step)): - u_star[i][j] = 0 - v_star[i][j] = 0 - p_new[i][j] = 100 - # Pressure correction p_c = np.zeros(shape=shape, dtype=float) for i in range(1, N - 1): @@ -126,6 +134,10 @@ while error > precision: a_W = -d_e[i][j-1] * step a_N = -d_n[i-1][j] * step a_S = -d_n[i][j] * step + assert_rule2(a_E) + assert_rule2(a_W) + assert_rule2(a_N) + assert_rule2(a_S) a_P = a_E + a_W + a_N + a_S b[i][j] = step * (-(u_star[i][j] - u_star[i][j-1]) + (v_star[i][j] - v_star[i-1][j])) @@ -146,6 +158,13 @@ while error > precision: u_new[i][j] = u_star[i][j] + alpha * d_e[i][j] * (p_c[i + 1][j] - p_c[i][j]) v_new[i][j] = v_star[i][j] + alpha * d_n[i][j] * (p_c[i][j] - p_c[i + 1][j]) + # Backwards-facing step boundary conditions enforce + for i in range(int(h_c / step)): + for j in range(int(l_c / step)): + u_new[i][j] = 0 + v_new[i][j] = 0 + p_new[i][j] = 0 + # Continuity residual as error measure error = 0 for i in range(N): |