Imposing reflective boundary conditions in the Linear Euler equations

variational formulation
1

Hello all,

I am trying to impose Reflective boundary conditions to the Linear Euler equations using reflective indices R_in and R_out at the inlet and outlet respectively. I am having an issue with currently defining this in the variational form as I don’t see any effect of changing R_out in the output pressure and velocity plots. I have attached my equations and the variational form that I am currently using. Please let me know your thoughts on how I can better implement the reflective indices in my boundary conditions. Thank you!

image
image822×241 9.41 KB

image
image706×421 17.3 KB

My boundary conditions are:

utflow = DirichletBC(P, 0.0, “x[0] > 1 - DOLFIN_EPS”)
p_in = Expression(“t < 1.0 ? 0.223.142sin(23.142*t) : 0”, t=0.0, degree=2)
inflow = DirichletBC(P, p_in, “x[0] < DOLFIN_EPS”)
bcp = [inflow, outflow]

My variational form is:

F1 = (inner((p_trial - p_old) / (dt * c * c), p_test) * dx + rho * u_old.dx(0) * p_test * dx
+ (c * (1 + R_in) / (1 - R_in) * p_trial) * p_test * ds(1)
+ (-c * (1 + R_out) / (1 - R_out) * p_trial) * p_test * ds(2))

F2 = inner(((u_trial - u_old) * rho) / dt, u_test) * dx + p_old.dx(0) * u_test * dx

2

Please read Read before posting: How do I get my question answered?

Without the actual code, there’s not much one can do. If I were to guess, I’d say that you forgot to assign the facet markers (1 and 2) when creating the ds measure.

3

Hi, thank you for your response. Here is my whole code. I have assigned the ds measure in this way. Please let me know if my approach is right.

EDITED BY FB

4

Please use “```” when writing code, otherwise the output is very difficult to read.

1 Like
6 
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from tqdm import tqdm

plt.rcParams['figure.figsize'] = [14, 8]
plt.rcParams.update({'font.size': 24})

set_log_level(30)

n_elements = 100
mesh = UnitIntervalMesh(n_elements)

U = FunctionSpace(mesh, "Lagrange", 1)
P = FunctionSpace(mesh, "Lagrange", 1)

u_trial = TrialFunction(U)
p_trial = TrialFunction(P)

u_test = TestFunction(U)
p_test = TestFunction(P)


u_initial_condition = Expression("0", degree=1)
p_initial_condition = Expression("0", degree=1)

u_old = interpolate(u_initial_condition, U)
p_old = interpolate(p_initial_condition, P)

EndT = 3  # sec
delta_t = 0.0001

dt = Constant(delta_t)

c = Constant(1)
rho = Constant(1.2)

heat_source = Constant(0.0)

R_in = 0.0
R_out = 0.0

outflow = DirichletBC(P, 0.0, "x[0] > 1 - DOLFIN_EPS")
p_in = Expression("t < 1.0 ? 0.2*2*3.14*2*sin(2*3.14*2*t) : 0", t=0.0, degree=2)
inflow = DirichletBC(P, p_in, "x[0] < DOLFIN_EPS")
bcp = [inflow, outflow]

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)

class InletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0, DOLFIN_EPS) and on_boundary

class OutletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0, DOLFIN_EPS) and on_boundary


inlet_boundary = InletBoundary()
inlet_boundary.mark(boundary_markers, 1)

outlet_boundary = OutletBoundary()
outlet_boundary.mark(boundary_markers, 2)

ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

F1 = (inner((p_trial - p_old) / (dt * c * c), p_test) * dx + rho * u_old.dx(0) * p_test * dx
      + (c * (1 + R_in) / (1 - R_in) * p_trial) * p_test * ds(1)
      + (-c * (1 + R_out) / (1 - R_out) * p_trial) * p_test * ds(2))

F2 = inner(((u_trial - u_old) * rho) / dt, u_test) * dx + p_old.dx(0) * u_test * dx

a1 = lhs(F1)
L1 = rhs(F1)

a2 = lhs(F2)
L2 = rhs(F2)

u_solution = Function(U)
p_solution = Function(P)

n_time_steps = int(np.ceil(EndT / delta_t))

u_plot = np.array(u_old.vector().get_local())
p_plot = np.array(p_old.vector().get_local())

x_coords = mesh.coordinates()[:, 0]

boundary_indices = np.where((x_coords < DOLFIN_EPS) | (x_coords > 1 - DOLFIN_EPS))[0]

gf_boundary_all = []

time_current = 0.0
for i in tqdm(range(n_time_steps)):
    time_current += delta_t

    p_in.t = time_current

    solve(a1 == L1, p_solution, bcp)
    solve(a2 == L2, u_solution)

    u_old.assign(u_solution)
    p_old.assign(p_solution)

    u_current = u_solution.vector().get_local()
    p_current = p_solution.vector().get_local()

    f = (p_current / (rho(0) * c(0))) + u_current
    g = (p_current / (rho(0) * c(0))) - u_current

    gf_boundary = np.zeros(len(boundary_indices))
    for idx, boundary_index in enumerate(boundary_indices):
        if abs(f[boundary_index]) > 1E-12:  # Avoid division by zero
            gf_boundary[idx] = abs(g[boundary_index] / f[boundary_index])

    gf_boundary_all.append(gf_boundary)

    u_plot = np.vstack((u_plot, u_current))
    p_plot = np.vstack((p_plot, p_current))

gf_boundary_all = np.array(gf_boundary_all)


plt.figure(1)
plt.contourf(u_plot[::100])
plt.colorbar().set_label('$v$ [-]', rotation=270, labelpad=40)
plt.savefig('./Velocity.png')

plt.figure(2)
plt.contourf(p_plot[::100])
plt.colorbar().set_label('$p$ [-]', rotation=270, labelpad=40)
plt.savefig('./Pressure.png')´´´´
7

Disclaimer: I have no experience with Euler equations.

Still, I think that there is something wrong with your formulation. You are imposing Dirichlet BCs on boundaries 1 and 2, and at the same time adding boundary integrals on those markers (like if they were Neumann BCs in an elliptic problem). The unknown is effectively assigned by the Dirichlet BC, hence it’s expected that you will not see any difference changing the problem definition (e.g., the coefficients).

8

Understood. But does my definition of the ds measure look alright to you?

9

Yes, the definition of the measure is correct.

10

Hello, thanks for your input! Unfortunately, removing the Dirichlet BC is giving rise to graphs with 0 pressure and velocity values for the entire range of R_in and R_out.

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from tqdm import tqdm

plt.rcParams['figure.figsize'] = [14, 8]
plt.rcParams.update({'font.size': 24})

set_log_level(30)

n_elements = 100
mesh = UnitIntervalMesh(n_elements)

# Define a Function Space
U = FunctionSpace(mesh, "Lagrange", 1)
P = FunctionSpace(mesh, "Lagrange", 1)

# Define trial and test functions
u_trial = TrialFunction(U)
p_trial = TrialFunction(P)

u_test = TestFunction(U)
p_test = TestFunction(P)

# The initial condition, u(t=0, x) = sin(pi * x)
u_initial_condition = Expression("0", degree=1)
p_initial_condition = Expression("sin(pi * x[0])", degree=1)
#p_initial_condition = Expression("0", degree=1)

# Discretize the initial condition
u_old = interpolate(u_initial_condition, U)
p_old = interpolate(p_initial_condition, P)

# The time stepping of the implicit Euler discretization (=dt)
EndT = 3  # sec
delta_t = 0.0001

dt = Constant(delta_t)

c = Constant(1)
rho = Constant(1.2)

# The forcing on the rhs of the PDE
heat_source = Constant(0.0)

R_in = 0.5
R_out = 0.5

# Boundary Conditions
#outflow = DirichletBC(P, 0.0, "x[0] > 1 - DOLFIN_EPS")
# inflow = DirichletBC(P, 0.0, "x[0] < DOLFIN_EPS")
p_in = Expression("t < 1.0 ? 0.2*2*3.14*2*sin(2*3.14*2*t) : 0", t=0.0, degree=2)
#inflow = DirichletBC(P, p_in, "x[0] < DOLFIN_EPS")
# bcp = [inflow]
#bcp = [inflow, outflow]

# Define outlet boundary
# Create a MeshFunction for marking boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)

class InletBoundary(SubDomain):
     def inside(self, x, on_boundary):
         return near(x[0], 0.0, DOLFIN_EPS) and on_boundary

 # Define outlet boundary
class OutletBoundary(SubDomain):
     def inside(self, x, on_boundary):
        return near(x[0], 1.0, DOLFIN_EPS) and on_boundary


# Initialize and mark the outlet boundary with a unique ID, e.g., 1
# inlet_boundary = InletBoundary()
# inlet_boundary.mark(boundary_markers, 1)

# outlet_boundary = OutletBoundary()
# outlet_boundary.mark(boundary_markers, 2)

def boundary(x, on_boundary):
    return on_boundary

# Define a Measure for the marked outlet boundary
ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

# Variational formulation
F1 = (inner((p_trial - p_old) / (dt * c * c), p_test) * dx + rho * u_old.dx(0) * p_test * dx
      + (c * (1 + R_in) / (1 - R_in) * p_trial) * p_test * ds(1)
      + (-c * (1 + R_out) / (1 - R_out) * p_trial) * p_test * ds(2))

F2 = inner(((u_trial - u_old) * rho) / dt, u_test) * dx + p_old.dx(0) * u_test * dx

a1 = lhs(F1)
L1 = rhs(F1)

a2 = lhs(F2)
L2 = rhs(F2)

u_solution = Function(U)
p_solution = Function(P)

# Time stepping
n_time_steps = int(np.ceil(EndT / delta_t))

u_plot = np.array(u_old.vector().get_local())
p_plot = np.array(p_old.vector().get_local())

x_coords = mesh.coordinates()[:, 0]

#Boundary points
boundary_indices = np.where((x_coords < DOLFIN_EPS) | (x_coords > 1 - DOLFIN_EPS))[0]

gf_boundary_all = []

time_current = 0.0
for i in tqdm(range(n_time_steps)):
    time_current += delta_t

    # Update pressure boundary condition
    p_in.t = time_current

    solve(a1 == L1, p_solution)#, bcp)
    solve(a2 == L2, u_solution)

    u_old.assign(u_solution)
    p_old.assign(p_solution)

    u_current = u_solution.vector().get_local()
    p_current = p_solution.vector().get_local()

    # Compute f and g for the current time step
    f = (p_current / (rho(0) * c(0))) + u_current
    g = (p_current / (rho(0) * c(0))) - u_current

    gf_boundary = np.zeros(len(boundary_indices))
    for idx, boundary_index in enumerate(boundary_indices):
        if abs(f[boundary_index]) > 1E-12:  # Avoid division by zero
            gf_boundary[idx] = abs(g[boundary_index] / f[boundary_index])

    gf_boundary_all.append(gf_boundary)

    u_plot = np.vstack((u_plot, u_current))
    p_plot = np.vstack((p_plot, p_current))

gf_boundary_all = np.array(gf_boundary_all)

plt.figure(1)
plt.contourf(u_plot[::100])
plt.colorbar().set_label('$v$ [-]', rotation=270, labelpad=40)
plt.savefig('./Velocity01.png')

plt.figure(2)
plt.contourf(p_plot[::100])
plt.colorbar().set_label('$p$ [-]', rotation=270, labelpad=40)
plt.savefig('./Pressure01.png')
11

Which is not surprising, considering that there is no forcing term. I’d suggest that you sort out the weak formulation first: I can’t help with that since I have no specific domain knowledge.