# Time loop
try:
while T > tsp and t < T:
timer_dt.start()
update_counter(counters)
# Update current time
t += tsp
# Update boundary conditions : only if time-dependent
# parabolic_profile.t = t; tim.t = t; num_cycle.cycle = int(t / t_period)
# for ui, value in inflow.items():
# inflow[ui][0].v = evaluate_boundary_val(param_LSPV); inflow[ui][1].v = evaluate_boundary_val(param_LIPV)
# inflow[ui][2].v = evaluate_boundary_val(param_RSPV); inflow[ui][3].v = evaluate_boundary_val(param_RIPV)
if problem_physics['solve_FSI'] == True:
timer_si.start()
Lm_f.assign(interpolate_nonmatching_mesh_delta(fsi_interpolation, Lm_[1], FS['fluid'][2], interpolation_fx, "F"))
si += timer_si.stop()
timer_s1.start()
# print(BLUE % "1: Predict tentative velocity step", flush = True)
A1, b1 = flow.assemble_tentative_velocity(u_, p_, Lm_f, dt)
flow.solve_tentative_velocity(A1, u_[0], b1, bcs['velocity'])
s1 += timer_s1.stop()
timer_s2.start()
# print(BLUE % "2: Pressure correction step", flush = True)
b2 = flow.assemble_pressure_correction(u_, p_, Lm_f, dt)
flow.solve_pressure_correction(p_[0], b2, bcs['pressure'])
s2 += timer_s2.stop()
timer_s3.start()
# print(BLUE % "3: Velocity correction step", flush = True)
b3 = flow.assemble_velocity_correction(u_, p_, dt)
flow.solve_velocity_correction(u_[0], b3, bcs['velocity'])
s3 += timer_s3.stop()
assigner_uv.assign(uv, [u_[0][ui] for ui in range(u_components)])
# ---------------------------------------------------------------------------------
if problem_physics['solve_FSI'] and problem_physics['solve_temperature'] == True:
timer_si.start()
LmTf_.assign(interpolate_nonmatching_mesh_delta(fsi_interpolation, LmTs_[1], FS['fluid_temp'][0], interpolation_fx, "F"))
si += timer_si.stop()
timer_s4.start()
# print(BLUE % "4: Energy conservation step", flush = True)
if problem_physics['solve_temperature'] == True:
A4, b4 = flow_temp.assemble_temperature(T_, uv, LmTf_, dt)
flow_temp.solve_temperature(A4, T_[0], b4, bcs['temperature'])
s4 += timer_s4.stop()
# ---------------------------------------------------------------------------------
if problem_physics['solve_FSI'] == True:
timer_si.start()
uf_.assign(interpolate_nonmatching_mesh_delta(fsi_interpolation, uv, FS['lagrange'][0], interpolation_fx, "S"))
si += timer_si.stop()
timer_s5.start()
# print(BLUE % "5: Solid momentum eq. step", flush = True)
if problem_physics['solve_FSI'] == True:
a5 = solid.assemble_solid_problem(problem_physics['compressible_solid'], Dp_, mix, uf_, Lm_[1], dt)
try:
solid.solve_solid_displacement(solid_mesh_R.mesh, problem_physics['compressible_solid'], a5, Dp_[1], mix, ps_, p_[0], bcs['solid'])
except:
solid.change_initial_guess(Dp_[1], mix)
solid.solve_solid_displacement(solid_mesh_R.mesh, problem_physics['compressible_solid'], a5, Dp_[1], mix, ps_, p_[0], bcs['solid'])
Dp_[0].vector().axpy(1.0, Dp_[1].vector())
# solid.compute_jacobian(J_, Dp_[0])
us_.vector().zero()
us_.vector().axpy(1/float(dt), Dp_[1].vector())
s5 += timer_s5.stop()
# ---------------------------------------------------------------------------------
timer_s6.start()
# print(BLUE % "6: Lagrange multiplier (fictitious force) step", flush = True)
if problem_physics['solve_FSI'] == True:
a6, b6 = lagrange.assemble_lagrange_multiplier(Lm_, us_, uf_, dt)
lagrange.solve_lagrange_multiplier(a6, Lm_[0], b6)
s6 += timer_s6.stop()
# ---------------------------------------------------------------------------------
if problem_physics['solve_FSI'] and problem_physics['solve_temperature'] == True:
timer_si.start()
Ts_[0].assign(interpolate_nonmatching_mesh_delta(fsi_interpolation, T_[0], FS['solid_temp'][1], interpolation_fx, "S"))
si += timer_si.stop()
timer_s7.start()
# print(BLUE % "7: Solid temperature based lagrange multiplier step", flush = True)
if problem_physics['solve_FSI'] and problem_physics['solve_temperature'] == True:
a7, b7 = solid_temp.assemble_solid_temperature_lagrange_multiplier(Ts_, uf_, dt)
solid_temp.solve_solid_temperature_lagrange_multiplier(a7, LmTs_[0], b7)
s7 += timer_s7.stop()
# Solver parameters
krylov_solvers=dict(
monitor_convergence=False,
report=False,
error_on_nonconvergence=True,
nonzero_initial_guess=True,
maximum_iterations=300,
absolute_tolerance=1e-8)
# Solver dictionaries
tentative_velocity_solver=dict(
solver_type='bicgstab',
preconditioner_type='jacobi')
solid_displacement_parameters = {"newton_solver":{"linear_solver":solid_momentum_solver['solver_type'], "preconditioner":'hypre_amg', "report":True, \
"error_on_nonconvergence":True, "absolute_tolerance":1e-15, "relative_tolerance":1e-6, "maximum_iterations":50}}
restart = False # Restart parameter
# Physics of the problem
# ---------------------------------------------------------------------
problem_physics = dict(
solve_temperature = True, # enter "True" if you want to solve for temperature
solve_FSI = True, # enter "True" if you want to solve for fluid-structure interaction
compressible_solid = True, # enter "True" if compressible: Also remember to specify compressibility (Ld)
# enter "False" if incompressible
viscous_dissipation = False, # Heat release due to viscous gradients
body_force = False, # Gravitational force (uniform volumetric force)
)
def f_dir(dim): # Body force direction : -ve y direction (by default)
n = -1*tensors.unit_vector(1, dim)
return n
interpolation_fx = 'phi4' # Delta-function interpolation for FSI problems
# FEM stabilization and constants
# ---------------------------------------------------------------------
stabilization_parameters = dict(
# Navier-stokes
SUPG_NS = False, # explicit
PSPG_NS = False, # explicit
crosswind_NS = False, # implicit
backflow_NS = False,
# Energy-equation
SUPG_HT = False, # explicit
crosswind_HT = False # implicit
)
alpha = Constant(0.85) # SUPG/PSPG stabilization constant
C_cw = Constant(0.7) # Crosswind stabilization constant (As per R Codina : quadratic elements: 0.35, for linear elements: 0.7)
# Physical parameters
# ---------------------------------------------------------------------
physical_parameters = dict(
g = 9.81, # Gravity (m/s2)
# Fluid
rho_f = 1, # Density (kg/m3)
nu = 1, # Dynamic viscosity (kg/m.s)
Spht_f = 1, # Specific heat (J/kg.C)
K_f = 1, # Thermal conductivity (W/m.C)
# Solid
rho_s = 10, # Density (kg/m3)
Sm = 0, # Shear modulus (N/m2)
Ld = 0, # Compressibility (N/m2)
Spht_s = 0.11, # Specific heat (J/kg.C)
K_s = 1.2 # Thermal conductivity (W/m.C)
)
def calc_non_dimensional_solid_properties(g, rho_f, nu, Spht_f, K_f, rho_s, Sm, Ld, Spht_s, K_s, Lsc, Vsc, T0, Tm, Tsc):
rho = rho_s/rho_f
Spht = Spht_s/Spht_f
K = K_s/K_f
Ld = Ld/(rho_f*Vsc*Vsc)
Sm = Sm/(rho_f*Vsc*Vsc)
return rho, Spht, K, Ld, Sm
# Characteristic scales
# ---------------------------------------------------------------------
characteristic_scales = dict(
Lsc = 1, # m
Vsc = 1, # m/s
T0 = -1*52, # lower_temp (C)
Tm = 37 # higher_temp (c)
)
# Temporal control
# ---------------------------------------------------------------------
time_control = dict(
C_no = 0.35, # Maximum possible Courant number
dt = 0.0025, # Time-step: constant throughout runtime if adjustable-timestep is "False"
T = 100, # Total runtime
adjustable_timestep = True # Calculate time-step using max Courant no. during runtime: used to accelerate temporal solution
)
# FEM degree of variables
# ---------------------------------------------------------------------
fem_degree = dict(
velocity_degree = 2,
pressure_degree = 1,
temperature_degree = 2,
displacement_degree = 2,
lagrange_degree = 1
)
# Non-dimensional numbers
# ---------------------------------------------------------------------
def calc_non_dimensional_numbers(g, rho_f, nu, Spht_f, K_f, rho_s, Sm, Ld, Spht_s, K_s, Lsc, Vsc, T0, Tm, Tsc):
Re = rho_f*(Vsc*Lsc)/nu
Pr = (Spht_f*nu)/K_f
Ec = (Vsc*Vsc)/(Spht_f*(Tm-T0))
Fr = Vsc/sqrt(g*Lsc)
return Re, Pr, Ec, Fr
# Enter "True" if you want to post-process data
# ---------------------------------------------------------------------
post_process = True
# File printing / solid-remeshing control
# ---------------------------------------------------------------------
print_control = dict(
a = 40, # for printing variables and restart files
b = 50, # for post processing data
c = 20, # for simulation_wall_time text file
d = 5, # for remeshing solid current-configuration mesh
e = 20 # for runtime_tsp_courant_no_stats text file
)
# If 2D problem?: Do u want to calculate stream function and vorticity! # Note to self: streamfunction is not defined for 3D.
# ---------------------------------------------------------------------
calc_stream_function = True
class PeriodicDomain(SubDomain):
def inside(self, x, on_boundary):
return bool(x[2] < DOLFIN_EPS and x[2] > -DOLFIN_EPS and on_boundary)
def map(self, x, y):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - 1.0
constrained_domain = PeriodicDomain() # None
parabolic_profile = Expression('6.0*x[1]*(4.1 - x[1])/(4.1*4.1)', degree=2)
class Point_pressure(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 4.2) and near(x[1], 5.) and near(x[2], 2.)
# Boundary conditions
def fluid_create_boundary_conditions(fluid_mesh, inflow, **V):
boundaries = fluid_mesh.get_mesh_boundaries()
# velocity
bcu_left_x = DirichletBC(V['fluid'][0], parabolic_profile, boundaries, 1)
bcu_bottom_x = DirichletBC(V['fluid'][0], Constant(0), boundaries, 2)
bcu_top_x = DirichletBC(V['fluid'][0], Constant(0), boundaries, 4)
bcu_x = [bcu_left_x, bcu_bottom_x, bcu_top_x]
bcu_left_y = DirichletBC(V['fluid'][0], Constant(0), boundaries, 1)
bcu_bottom_y = DirichletBC(V['fluid'][0], Constant(0), boundaries, 2)
bcu_top_y = DirichletBC(V['fluid'][0], Constant(0), boundaries, 4)
bcu_y = [bcu_left_y, bcu_bottom_y, bcu_top_y]
bcu_left_z = DirichletBC(V['fluid'][0], Constant(0), boundaries, 1)
bcu_bottom_z = DirichletBC(V['fluid'][0], Constant(0), boundaries, 2)
bcu_z = [bcu_left_z, bcu_bottom_z]
bcu = [bcu_x, bcu_y, bcu_z]
# pressure
bcp_right = DirichletBC(V['fluid'][1], Constant(0), boundaries, 3)
bcp = [bcp_right]
# Streamfunction
wall = 'on_boundary'
bcPSI = DirichletBC(V['fluid'][1], 0, wall)
bcs = dict(velocity = bcu, pressure = bcp, streamfunction = bcPSI)
if problem_physics['solve_temperature'] == True:
# temperature
bcT_left = DirichletBC(V['fluid_temp'][0], Constant(1), boundaries, 1)
bcT_top = DirichletBC(V['fluid_temp'][0], Constant(0), boundaries, 4)
bcT = [bcT_left, bcT_top]
bcs.update(temperature = bcT)
return bcs
def solid_create_boundary_conditions(solid_mesh_R, compressible_solid, dt, **V):
boundaries = solid_mesh_R.get_mesh_boundaries()
# Note to self: Boundary conditions are for incremental displacement (delta D)
# Solid
if compressible_solid == False:
bcx_cylinder = DirichletBC(V['solid'][1].sub(0), Constant((0, 0, 0)), boundaries, 1)
elif compressible_solid == True:
bcx_cylinder = DirichletBC(V['solid'][0], Constant((0, 0, 0)), boundaries, 1)
bcx = [bcx_cylinder]
return bcx
# Initial conditions
def fluid_create_initial_conditions(u_, p_, T_):
# Velocity / pressure
for i in range(3):
u_[i][0].vector()[:] = 0.0
u_[i][1].vector()[:] = 0.0
u_[i][2].vector()[:] = 0.0
p_[i].vector()[:] = 0.0
# Temperature
for i in range(3):
T_[i].vector()[:] = 0.0
def solid_create_initial_conditions(Dp_, mix, dt):
# Solid pressure (only defined for incompressible solid)
assign(mix.sub(1), interpolate(Constant(0), mix.sub(1).function_space().collapse()))
# Cumulative displacement
Dp_[0].vector()[:] = 0.0
# Incremental displacement (delta D)
Dp_[1].vector()[:] = 0.0 # V_init*dt
Dp_[2].vector()[:] = 0.0 # V_init*dt
assign(mix.sub(0), interpolate(Expression(('0.0', '0.0', '0.0'), degree = 2), mix.sub(0).function_space().collapse()))