From 011f7ae73560b690aa1af58d445e1cf0afdf7103 Mon Sep 17 00:00:00 2001
From: eug-vs <eugene@eug-vs.xyz>
Date: Thu, 26 May 2022 13:52:47 +0400
Subject: fix: account for mass source in momentum equations
---
src/simple.py | 15 ++++++++-------
1 file changed, 8 insertions(+), 7 deletions(-)
(limited to 'src/simple.py')
diff --git a/src/simple.py b/src/simple.py
index 3318213..8523db8 100644
--- a/src/simple.py
+++ b/src/simple.py
@@ -2,12 +2,13 @@ import numpy as np
class SIMPLE:
- def __init__(self, shape, bfs_shape, step, Re, alpha=0.8):
+ def __init__(self, shape, bfs_shape, step, Re, alpha_p=0.8, alpha_uv=0.5):
np.set_printoptions(precision=2, floatmode="maxprec", suppress=True)
self.Re = Re
self.nu = 1 / Re
- self.alpha = alpha
+ self.alpha_p = alpha_p
+ self.alpha_uv = alpha_uv
self.step = step
self.bfs_shape = bfs_shape
@@ -58,7 +59,7 @@ class SIMPLE:
a_W * self.u[i][j - 1] +
a_N * self.u[i + 1][j] +
a_S * self.u[i - 1][j] +
- 0#self.b[i][j - 1]
+ self.b[i][j - 1]
) / a_e + self.d_e[i][j] * (self.p_star[i][j - 1] - self.p_star[i][j]) # p - p_e
# Momentum along Y direction
@@ -86,7 +87,7 @@ class SIMPLE:
a_W * self.v[i][j - 1] +
a_N * self.v[i + 1][j] +
a_S * self.v[i - 1][j] +
- 0#self.b[i - 1][j]
+ self.b[i - 1][j]
) / a_n + self.d_n[i][j] * (self.p_star[i - 1][j] - self.p_star[i][j]) # p - p_n
def correct_pressure(self):
@@ -113,17 +114,17 @@ class SIMPLE:
self.b[i][j]
) / a_P
- self.p = self.p_star + self.p_prime * self.alpha
+ self.p = self.p_star + self.p_prime * self.alpha_p
self.p_star = self.p
def correct_velocities(self):
for i in range(self.u.shape[0]):
for j in range(1, self.u.shape[1] - 1):
- self.u[i][j] = self.u_star[i][j] + self.d_e[i][j] * (self.p_prime[i][j - 1] - self.p_prime[i][j])
+ self.u[i][j] = self.u_star[i][j] + self.alpha_uv * self.d_e[i][j] * (self.p_prime[i][j - 1] - self.p_prime[i][j])
for i in range(1, self.v.shape[0] - 1):
for j in range(self.v.shape[1]):
- self.v[i][j] = self.v_star[i][j] + self.d_n[i][j] * (self.p_prime[i - 1][j] - self.p_prime[i][j])
+ self.v[i][j] = self.v_star[i][j] + self.alpha_uv * self.d_n[i][j] * (self.p_prime[i - 1][j] - self.p_prime[i][j])
def iterate(self):
self.solve_momentum_equations()
--
cgit v1.2.3