10. Variable Density Flow

We study the solution of the Poisson equation in multi-material flows. The following example shows a 1D problem in a domain \([0,1]\), where the unsimulated air phase occupies the region \([0,0.1]\), incompressible water phase occupies \([0.1,0.6]\times[0.7,1]\), and an incompressible bubble phase occupies \([0.6,0.7]\).

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

cell_count = 100;
pressure = [None]*cell_count;
dx = 1./cell_count;

# create flags
flags = [None]*cell_count;
for i in range(0,cell_count):
    location = (i+.5)*dx
    if(location>=0 and location<=.1):
        flags[i]=0;     # unsimulated air
    elif(location>.1 and location<=.6):
        flags[i]=1;     # water
    elif(location>.6 and location<=.7):
        flags[i]=2;     # bubble
    elif(location>.7 and location<=1):
        flags[i]=1;     # water

# mark variables
index=0
variable_index = [None]*cell_count;
for i in range(0,cell_count):
    if(flags[i]>0):
        variable_index[i]=index;
        index+=1;
    else: variable_index[i]=-1;

number_of_dofs = index;
density = [None]*cell_count;
for i in range(0,cell_count):
    if(flags[i]==1):
        density[i]=10.;         # water
    elif(flags[i]==2):
        density[i]=1000.;       # bubble
    else: density[i]=0;

face_weights = [None]*(cell_count-1);
for i in range(0,cell_count-1):
    left_cell = i;
    right_cell = i+1;
    if(flags[left_cell]==0 and flags[right_cell]==0):
        face_weights[i]=0.;
    elif(flags[left_cell]==0 and flags[right_cell]==1):
        face_weights[i]=1./density[right_cell];
    else: face_weights[i]=2./(density[left_cell]+density[right_cell]);

# create Poisson matrix
A = lil_matrix((number_of_dofs,number_of_dofs));
one_over_dx2 = (1./dx)*(1./dx);
for i in range(0,cell_count):
    if(flags[i]>0):
        A[variable_index[i],variable_index[i]] += face_weights[i-1]*one_over_dx2;
        if(flags[i-1]>0): A[variable_index[i],variable_index[i-1]] -= face_weights[i-1]*one_over_dx2;
        if(i<cell_count-1):
            A[variable_index[i],variable_index[i]] += face_weights[i]*one_over_dx2;
            A[variable_index[i],variable_index[i+1]] -= face_weights[i]*one_over_dx2;
A = A.tocsr()

# create right hand side
b = np.zeros(shape=(number_of_dofs,1))
b[number_of_dofs-2] = 1.;

# solve Ax = b
x = spsolve(A,b);

# output variables
coord = [None]*cell_count;
for i in range(0,cell_count):
    coord[i] = (i+.5)*dx;
    if(flags[i]==0): pressure[i]=0.;
    else: pressure[i]=x[variable_index[i]];

plt.plot(coord,pressure,'ro')
plt.show()