822 lines
38 KiB
Python
822 lines
38 KiB
Python
#!/usr/bin/env python3
|
|
"""Generate a Word document describing the pipe system simulation theory and implementation."""
|
|
|
|
import os
|
|
import sys
|
|
|
|
from docx import Document
|
|
from docx.shared import Pt, Inches, Cm, RGBColor
|
|
from docx.enum.text import WD_ALIGN_PARAGRAPH
|
|
from docx.enum.table import WD_TABLE_ALIGNMENT
|
|
from docx.oxml.ns import qn
|
|
|
|
|
|
def set_cell_shading(cell, color_hex):
|
|
shading_elm = cell._element.get_or_add_tcPr()
|
|
shd = shading_elm.makeelement(qn('w:shd'), {
|
|
qn('w:val'): 'clear',
|
|
qn('w:color'): 'auto',
|
|
qn('w:fill'): color_hex,
|
|
})
|
|
shading_elm.append(shd)
|
|
|
|
|
|
def add_equation(doc, text, bold=False):
|
|
"""Add a centered equation paragraph with Cambria Math font."""
|
|
p = doc.add_paragraph()
|
|
p.alignment = WD_ALIGN_PARAGRAPH.CENTER
|
|
p.paragraph_format.space_before = Pt(6)
|
|
p.paragraph_format.space_after = Pt(6)
|
|
run = p.add_run(text)
|
|
run.font.name = 'Cambria Math'
|
|
run.font.size = Pt(11)
|
|
if bold:
|
|
run.bold = True
|
|
return p
|
|
|
|
|
|
def add_code_block(doc, code):
|
|
"""Add a code block with monospace font and gray background."""
|
|
p = doc.add_paragraph()
|
|
p.paragraph_format.space_before = Pt(4)
|
|
p.paragraph_format.space_after = Pt(4)
|
|
p.paragraph_format.left_indent = Cm(1)
|
|
run = p.add_run(code)
|
|
run.font.name = 'Consolas'
|
|
run.font.size = Pt(9)
|
|
run.font.color.rgb = RGBColor(0x33, 0x33, 0x33)
|
|
return p
|
|
|
|
|
|
def make_table(doc, headers, rows, col_widths=None):
|
|
"""Create a formatted table."""
|
|
table = doc.add_table(rows=1 + len(rows), cols=len(headers))
|
|
table.style = 'Table Grid'
|
|
table.alignment = WD_TABLE_ALIGNMENT.CENTER
|
|
|
|
# Header row
|
|
for j, h in enumerate(headers):
|
|
cell = table.rows[0].cells[j]
|
|
cell.text = h
|
|
set_cell_shading(cell, 'D9E2F3')
|
|
for paragraph in cell.paragraphs:
|
|
for run in paragraph.runs:
|
|
run.bold = True
|
|
run.font.size = Pt(10)
|
|
|
|
# Data rows
|
|
for i, row in enumerate(rows):
|
|
for j, val in enumerate(row):
|
|
cell = table.rows[i + 1].cells[j]
|
|
cell.text = str(val)
|
|
for paragraph in cell.paragraphs:
|
|
for run in paragraph.runs:
|
|
run.font.size = Pt(10)
|
|
|
|
if col_widths:
|
|
for i, w in enumerate(col_widths):
|
|
for row in table.rows:
|
|
row.cells[i].width = Cm(w)
|
|
|
|
return table
|
|
|
|
|
|
def build_document():
|
|
doc = Document()
|
|
|
|
# --- Global style ---
|
|
style = doc.styles['Normal']
|
|
style.font.name = 'Times New Roman'
|
|
style.font.size = Pt(11)
|
|
style.paragraph_format.line_spacing = 1.15
|
|
style.paragraph_format.space_after = Pt(6)
|
|
|
|
for level in range(1, 4):
|
|
hs = doc.styles[f'Heading {level}']
|
|
hs.font.name = 'Times New Roman'
|
|
hs.font.color.rgb = RGBColor(0x1F, 0x3A, 0x5F)
|
|
|
|
# ========================================================================
|
|
# TITLE
|
|
# ========================================================================
|
|
title = doc.add_heading('0D-1D Tank-Pipe Coupled Transient Simulation\nTechnical Documentation', level=0)
|
|
title.alignment = WD_ALIGN_PARAGRAPH.CENTER
|
|
for run in title.runs:
|
|
run.font.size = Pt(22)
|
|
|
|
p = doc.add_paragraph()
|
|
p.alignment = WD_ALIGN_PARAGRAPH.CENTER
|
|
run = p.add_run('0D-1D Tank-Pipe Blowdown Simulation\nTheory, Numerical Methods, and Source Code Description')
|
|
run.font.size = Pt(12)
|
|
run.font.color.rgb = RGBColor(0x66, 0x66, 0x66)
|
|
run.italic = True
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 1. SYSTEM OVERVIEW
|
|
# ========================================================================
|
|
doc.add_heading('1 System Overview', level=1)
|
|
|
|
doc.add_paragraph(
|
|
'This document describes the theoretical foundations, numerical schemes, '
|
|
'and source code implementation of a 0D-1D coupled transient simulation system. '
|
|
'The system models the blowdown process where gas flows from a high-pressure '
|
|
'tank (Tank 1) through a slender pipe into a low-pressure tank (Tank 2).'
|
|
)
|
|
|
|
doc.add_heading('1.1 Physical Scenario', level=2)
|
|
doc.add_paragraph(
|
|
'The system consists of three components:\n'
|
|
'(1) High-pressure tank (upstream, Tank 1): volume V1 = 5 m^3, initial pressure P1 = 5 MPa, temperature T1 = 300 K.\n'
|
|
'(2) Pipe: length L = 1 m, inner diameter D = 5 mm, divided into N = 20 finite-volume cells.\n'
|
|
'(3) Low-pressure tank (downstream, Tank 2): volume V2 = 10 m^3, initial pressure P2 = 2 MPa, temperature T2 = 300 K.\n'
|
|
'\n'
|
|
'The two tanks are connected through the pipe. Due to the pressure difference, '
|
|
'gas flows from the high-pressure side to the low-pressure side. '
|
|
'The flow inside the pipe is one-dimensional compressible flow, while the two tanks, '
|
|
'whose volumes are much larger than the pipe, have spatially uniform internal states '
|
|
'and are modeled with 0D lumped-parameter models.'
|
|
)
|
|
|
|
doc.add_heading('1.2 Modeling Strategy', level=2)
|
|
doc.add_paragraph(
|
|
'The overall approach is "0D-1D coupling":\n'
|
|
'- Tanks: 0D Lumped Parameter Model -- internal state (P, T, rho) is spatially uniform, '
|
|
'varies only with time.\n'
|
|
'- Pipe: 1D Compressible Euler Equations -- captures pressure waves, shock waves, '
|
|
'and expansion waves propagating through the pipe.\n'
|
|
'- Coupling: Ghost Cell Method -- a virtual cell is placed at each end of the pipe, '
|
|
'filled with the current tank state, then the HLL Riemann solver computes boundary fluxes.\n'
|
|
'- Friction: Source Term Method (Operator Splitting) -- Darcy-Weisbach wall friction '
|
|
'added to the momentum equation.'
|
|
)
|
|
|
|
doc.add_heading('1.3 Code Module Structure', level=2)
|
|
|
|
make_table(doc,
|
|
['Module', 'Responsibility', 'Core Class/Function'],
|
|
[
|
|
['config.py', 'Physical constants, geometry, simulation control', '14 constants + assertion checks'],
|
|
['tank.py', '0D tank model', 'Tank class (mass, U, ghost_state, apply_flux)'],
|
|
['riemann.py', 'HLL Riemann numerical flux', 'hll_flux(W_L, W_R, gamma)'],
|
|
['pipe.py', '1D pipe finite-volume model', 'Pipe class (W, primitives, step)'],
|
|
['friction.py', 'Darcy-Weisbach friction factor', 'darcy_friction_factor(Re, eps_D)'],
|
|
['solver.py', 'Time-stepping driver (coupling orchestrator)', 'run(tank1, tank2, pipe, ...)'],
|
|
['output.py', 'Post-processing: storage, plots, animation, report', 'save_history, plot_*, make_pipe_animation'],
|
|
['main.py', 'Entry point: assemble, run, verify, output', 'main()'],
|
|
],
|
|
col_widths=[3.5, 6, 6.5],
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 2. 0D TANK MODEL
|
|
# ========================================================================
|
|
doc.add_heading('2 0D Tank Model', level=1)
|
|
|
|
doc.add_heading('2.1 Basic Assumptions', level=2)
|
|
doc.add_paragraph(
|
|
'The tanks use a lumped-parameter (0D) model with the following assumptions:\n'
|
|
'(1) The gas is an ideal gas with equation of state P = rho * R * T.\n'
|
|
'(2) Internal state is spatially uniform -- density rho, temperature T, and pressure P '
|
|
'are functions of time only.\n'
|
|
'(3) The macroscopic velocity inside the tank is zero (u = 0); kinetic energy is negligible.\n'
|
|
'(4) The tank walls are adiabatic; no heat conduction to the environment.\n'
|
|
'(5) The tank volume is fixed.'
|
|
)
|
|
|
|
doc.add_heading('2.2 Conservation Equations', level=2)
|
|
doc.add_paragraph(
|
|
'For an open-system tank, mass conservation and energy conservation '
|
|
'(open-system first law of thermodynamics) are:'
|
|
)
|
|
add_equation(doc, 'dm/dt = m_dot_in')
|
|
add_equation(doc, 'dU/dt = H_dot_in = m_dot_in * h_t,in')
|
|
doc.add_paragraph(
|
|
'Where:\n'
|
|
'- m = rho * V is the total gas mass in the tank [kg]\n'
|
|
'- U = P*V/(gamma-1) is the total internal energy [J] (since u=0, total energy = internal energy)\n'
|
|
'- m_dot_in is the mass flow rate entering the tank [kg/s]\n'
|
|
'- H_dot_in is the total enthalpy flow rate entering the tank [W]\n'
|
|
'- h_t,in is the specific total enthalpy of the incoming gas [J/kg]\n'
|
|
'\n'
|
|
'Key design: In the code, the Tank primary state is (mass, U), '
|
|
'while rho, T, P are all computed on demand via @property:'
|
|
)
|
|
add_equation(doc, 'rho = mass / V')
|
|
add_equation(doc, 'T = (U / mass) * (gamma - 1) / R')
|
|
add_equation(doc, 'P = rho * R * T')
|
|
doc.add_paragraph(
|
|
'This design avoids state-synchronization bugs after time-step updates -- '
|
|
'only mass and U are updated; all derived quantities are automatically consistent.'
|
|
)
|
|
|
|
doc.add_heading('2.3 Code Implementation (tank.py)', level=2)
|
|
doc.add_paragraph('The Tank constructor computes initial mass and U from (V, P_init, T_init, gamma, R):')
|
|
add_code_block(doc,
|
|
'rho = P_init / (R_gas * T_init)\n'
|
|
'self.mass = rho * V\n'
|
|
'self.U = P_init * V / (gamma - 1)')
|
|
|
|
doc.add_paragraph(
|
|
'The apply_flux(mdot, edot, dt, sign) method performs each time-step update:\n'
|
|
' mass += sign * mdot * dt\n'
|
|
' U += sign * edot * dt\n'
|
|
'Where sign = -1 means gas flows out (upstream Tank 1), '
|
|
'sign = +1 means gas flows in (downstream Tank 2).'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 3. 1D PIPE MODEL
|
|
# ========================================================================
|
|
doc.add_heading('3 1D Pipe Model', level=1)
|
|
|
|
doc.add_heading('3.1 Governing Equations: 1D Compressible Euler Equations', level=2)
|
|
doc.add_paragraph(
|
|
'Gas flow inside the pipe is described by the 1D compressible Euler equations '
|
|
'(without friction):'
|
|
)
|
|
add_equation(doc, 'dW/dt + dF(W)/dx = 0')
|
|
doc.add_paragraph('Where the conservative variable vector W and flux vector F(W) are:')
|
|
add_equation(doc, 'W = [rho, rho*u, rho*E]^T')
|
|
add_equation(doc, 'F = [rho*u, rho*u^2 + P, u*(rho*E + P)]^T')
|
|
doc.add_paragraph(
|
|
'Physical meaning of each component:\n'
|
|
'- W[0] = rho: density [kg/m^3]\n'
|
|
'- W[1] = rho*u: momentum density [kg/(m^2*s)]\n'
|
|
'- W[2] = rho*E: total energy density [J/m^3], where E = e + u^2/2 is specific total energy, '
|
|
'e = P/(rho*(gamma-1)) is specific internal energy\n'
|
|
'\n'
|
|
'The three components of F correspond to mass flux, momentum flux (including pressure), '
|
|
'and energy flux, respectively.'
|
|
)
|
|
doc.add_paragraph(
|
|
'Equation of state (ideal gas closure):'
|
|
)
|
|
add_equation(doc, 'P = (gamma - 1) * (rho*E - 0.5*rho*u^2)')
|
|
add_equation(doc, 'a = sqrt(gamma * P / rho) [speed of sound]')
|
|
|
|
doc.add_heading('3.2 Finite Volume Discretization', level=2)
|
|
doc.add_paragraph(
|
|
'The pipe is uniformly divided into N finite-volume cells, each of width dx = L/N. '
|
|
'Within each cell, the conservative variables W take cell-averaged values '
|
|
'(piecewise-constant reconstruction), giving first-order spatial accuracy.'
|
|
)
|
|
doc.add_paragraph(
|
|
'For the i-th cell [x_{i-1/2}, x_{i+1/2}], integrating the conservation law '
|
|
'yields the semi-discrete form:'
|
|
)
|
|
add_equation(doc, 'dW_i/dt = -(1/dx) * [F_{i+1/2} - F_{i-1/2}]')
|
|
doc.add_paragraph(
|
|
'Where F_{i+1/2} is the numerical flux at the interface between cells i and i+1. '
|
|
'The key challenge is that adjacent cell states are generally discontinuous at the interface '
|
|
'(a Riemann discontinuity), and a Riemann solver is needed to determine the interface flux.'
|
|
)
|
|
|
|
doc.add_heading('3.3 Code Implementation (pipe.py)', level=2)
|
|
doc.add_paragraph(
|
|
'The Pipe class stores the W array with shape (3, N), where each column is one cell\'s '
|
|
'conservative variables. The initial state is computed from (P_init, T_init):'
|
|
)
|
|
add_code_block(doc,
|
|
'rho = P_init / (R_gas * T_init)\n'
|
|
'E_density = P_init / (gamma - 1) # u=0, total energy density = internal\n'
|
|
'W[0, :] = rho # density\n'
|
|
'W[1, :] = 0.0 # momentum density (u=0)\n'
|
|
'W[2, :] = E_density # total energy density')
|
|
doc.add_paragraph(
|
|
'The primitives() method recovers primitive variables (rho, u, P, a) from W. '
|
|
'The step(flux_L, flux_R, dt) method performs one time-step advancement.'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 4. HLL RIEMANN SOLVER
|
|
# ========================================================================
|
|
doc.add_heading('4 HLL Riemann Solver', level=1)
|
|
|
|
doc.add_heading('4.1 The Riemann Problem', level=2)
|
|
doc.add_paragraph(
|
|
'In the finite-volume method, the states on the left and right sides of each cell interface '
|
|
'are generally different, forming a local Riemann problem. For the 1D Euler equations, '
|
|
'the exact Riemann solution contains three waves: a left-going shock/rarefaction, '
|
|
'a contact discontinuity, and a right-going shock/rarefaction.'
|
|
'\n\n'
|
|
'Exact Riemann solvers are computationally expensive, so engineering practice typically uses '
|
|
'approximate Riemann solvers. This system uses the HLL (Harten-Lax-van Leer) scheme, '
|
|
'which is simple and robust.'
|
|
)
|
|
|
|
doc.add_heading('4.2 HLL Scheme Derivation', level=2)
|
|
doc.add_paragraph(
|
|
'The HLL scheme assumes the Riemann fan is separated by two waves (S_L and S_R) '
|
|
'into three regions:\n'
|
|
'- x/t < S_L: left state W_L (undisturbed region ahead of waves)\n'
|
|
'- S_L < x/t < S_R: intermediate state W* ("HLL average state")\n'
|
|
'- x/t > S_R: right state W_R (undisturbed region ahead of waves)\n'
|
|
'\n'
|
|
'The flux at the interface x=0 depends on the signs of S_L and S_R:'
|
|
)
|
|
|
|
doc.add_paragraph('Case 1: If S_L >= 0 (both waves travel rightward), the interface is left of all waves:')
|
|
add_equation(doc, 'F_HLL = F(W_L)')
|
|
|
|
doc.add_paragraph('Case 2: If S_R <= 0 (both waves travel leftward), the interface is right of all waves:')
|
|
add_equation(doc, 'F_HLL = F(W_R)')
|
|
|
|
doc.add_paragraph('Case 3: If S_L < 0 < S_R (interface is between the two waves), weighted average:')
|
|
add_equation(doc, 'F_HLL = [S_R*F_L - S_L*F_R + S_L*S_R*(W_R - W_L)] / (S_R - S_L)')
|
|
|
|
doc.add_heading('4.3 Wave Speed Estimates', level=2)
|
|
doc.add_paragraph(
|
|
'The choice of wave speeds S_L and S_R is critical to the HLL scheme. '
|
|
'This system uses the Davis estimate:'
|
|
)
|
|
add_equation(doc, 'S_L = min(u_L - a_L, u_R - a_R)')
|
|
add_equation(doc, 'S_R = max(u_L + a_L, u_R + a_R)')
|
|
doc.add_paragraph(
|
|
'Where u is the flow velocity and a = sqrt(gamma*P/rho) is the local speed of sound. '
|
|
'This estimate is simple and robust, ensuring all physical signal propagation speeds '
|
|
'are contained within [S_L, S_R].'
|
|
)
|
|
|
|
doc.add_heading('4.4 Code Implementation (riemann.py)', level=2)
|
|
doc.add_paragraph('The hll_flux(W_L, W_R, gamma) function follows this procedure:')
|
|
doc.add_paragraph(
|
|
'(1) Recover primitive variables (rho, u, P, a) from W_L and W_R\n'
|
|
'(2) Compute physical fluxes F_L = F(W_L), F_R = F(W_R)\n'
|
|
'(3) Estimate wave speeds S_L, S_R\n'
|
|
'(4) Return the HLL flux according to the three cases above\n'
|
|
'(5) Input validation: raise ValueError if rho <= 0 or P <= 0'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 5. GHOST CELL BOUNDARY -- THE KEY COUPLING MECHANISM
|
|
# ========================================================================
|
|
doc.add_heading('5 Boundary Conditions: Ghost Cell Method', level=1)
|
|
|
|
doc.add_paragraph(
|
|
'This is one of the most critical design elements of the simulation: how to couple '
|
|
'the 0D tank models with the 1D pipe model. The approach used here is the '
|
|
'Ghost Cell Method. The basic idea is to place a "virtual cell" at each end of the pipe, '
|
|
'fill it with the current tank state, and then call the Riemann solver between the '
|
|
'virtual cell and the first/last real pipe cell to compute the boundary numerical flux.'
|
|
)
|
|
|
|
doc.add_heading('5.1 Ghost State Construction', level=2)
|
|
doc.add_paragraph(
|
|
'For a tank connected to the pipe, the Ghost State is generated by '
|
|
'Tank.ghost_state(). The conservative variables of the ghost cell are:'
|
|
)
|
|
add_equation(doc, 'W_ghost = [rho_tank, 0, P_tank / (gamma-1)]^T')
|
|
doc.add_paragraph(
|
|
'Where:\n'
|
|
'- rho_tank is the current tank density\n'
|
|
'- Momentum rho*u = 0 (stagnation assumption: macroscopic velocity inside the tank is zero)\n'
|
|
'- rho*E = P/(gamma-1) (since u=0, total energy density = internal energy density)\n'
|
|
'\n'
|
|
'The "Stagnation Assumption" is the key here:\n'
|
|
'The tank volume is much larger than the pipe cross-section times pipe length, '
|
|
'so the bulk gas velocity inside the tank is extremely low -- the tank side of the pipe '
|
|
'inlet can be treated as a stagnation state. '
|
|
'The Riemann solver automatically determines the correct flow direction and flux magnitude '
|
|
'based on the state difference between the pipe side and the tank side.'
|
|
)
|
|
|
|
doc.add_heading('5.2 Left Boundary (Tank 1 -> Pipe Inlet)', level=2)
|
|
doc.add_paragraph('The specific computation steps:')
|
|
doc.add_paragraph(
|
|
'(1) Get ghost state from Tank 1: W_ghost_L = tank1.ghost_state()\n'
|
|
'(2) First real pipe cell state: W_pipe_0 = pipe.W[:, 0]\n'
|
|
'(3) Call HLL solver: flux_L = hll_flux(W_ghost_L, W_pipe_0, gamma)\n'
|
|
'(4) Here HLL "left state" = tank ghost cell, "right state" = pipe cell 0\n'
|
|
'\n'
|
|
'The returned flux_L is a 3-component vector [mass flux, momentum flux, energy flux], '
|
|
'with positive direction from left to right (i.e., from tank into pipe).'
|
|
)
|
|
|
|
p = doc.add_paragraph()
|
|
run = p.add_run(
|
|
'Physical interpretation: Since Tank 1 pressure (5 MPa) is much higher than '
|
|
'the pipe initial pressure (2 MPa), the Riemann solver automatically produces '
|
|
'a positive left-to-right flux, driving gas from the high-pressure tank into the pipe.'
|
|
)
|
|
run.italic = True
|
|
|
|
doc.add_heading('5.3 Right Boundary (Pipe Outlet -> Tank 2)', level=2)
|
|
doc.add_paragraph(
|
|
'(1) Last real pipe cell state: W_pipe_{N-1} = pipe.W[:, -1]\n'
|
|
'(2) Get ghost state from Tank 2: W_ghost_R = tank2.ghost_state()\n'
|
|
'(3) Call HLL solver: flux_R = hll_flux(W_pipe_{N-1}, W_ghost_R, gamma)\n'
|
|
'(4) Here HLL "left state" = last pipe cell, "right state" = tank ghost cell\n'
|
|
'\n'
|
|
'flux_R positive direction is also left to right (i.e., from pipe into tank).'
|
|
)
|
|
|
|
doc.add_heading('5.4 Schematic Diagram', level=2)
|
|
doc.add_paragraph(
|
|
'The following diagram shows the spatial layout of the ghost cell method:'
|
|
)
|
|
add_code_block(doc,
|
|
' [Tank 1] | Cell 0 | Cell 1 | ... | Cell N-1 | [Tank 2]\n'
|
|
' (ghost_L) | | | | | (ghost_R)\n'
|
|
' ^ ^\n'
|
|
' flux_L flux_R\n'
|
|
' = hll(ghost_L, = hll(pipe[:,-1],\n'
|
|
' pipe[:,0]) ghost_R)\n'
|
|
)
|
|
doc.add_paragraph(
|
|
'Note: The ghost cells do not occupy physical space. They only provide '
|
|
'"outside-the-pipe" information for the HLL solver to correctly compute boundary fluxes. '
|
|
'The ghost cell W values are re-frozen (snapshot) at the start of each time step '
|
|
'to ensure boundary conditions remain consistent within a single time step.'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 6. FLUX DOUBLING
|
|
# ========================================================================
|
|
doc.add_heading('6 Flux Doubling: Guaranteeing Conservation', level=1)
|
|
|
|
doc.add_paragraph(
|
|
'This is the core mechanism ensuring physical conservation laws '
|
|
'(mass conservation, energy conservation) in the entire simulation system.'
|
|
)
|
|
|
|
doc.add_heading('6.1 The Problem', level=2)
|
|
doc.add_paragraph(
|
|
'In 0D-1D coupling, the pipe boundary flux must be used not only to update '
|
|
'the pipe internal cell states, but also to simultaneously update the adjacent '
|
|
'tank\'s mass and energy. If the pipe and the tank use different flux calculations '
|
|
'to update their own states, then the total system mass and energy cannot be exactly '
|
|
'conserved -- an unphysical phenomenon of "mass appearing or disappearing" '
|
|
'would occur at the interface.'
|
|
)
|
|
|
|
doc.add_heading('6.2 Solution: Same Flux Updates Both Sides', level=2)
|
|
doc.add_paragraph(
|
|
'The approach in this system is: compute the boundary HLL flux only once, '
|
|
'then use the same flux vector to simultaneously update:\n'
|
|
'(1) The pipe boundary cell via the finite-volume formula (in pipe.step using flux_L, flux_R)\n'
|
|
'(2) The tank mass and energy (in tank.apply_flux)\n'
|
|
'\n'
|
|
'This is called "Flux Doubling".'
|
|
)
|
|
|
|
doc.add_paragraph('Using the left boundary as an example, the key logic in solver.py:')
|
|
add_code_block(doc,
|
|
'# Phase 3: compute left boundary flux (only once)\n'
|
|
'flux_L = hll_flux(W_ghost_L, pipe.W[:, 0], gamma)\n'
|
|
'\n'
|
|
'# Phase 4: pipe uses this flux for update (as left interface flux)\n'
|
|
'pipe.step(flux_L, flux_R, dt)\n'
|
|
' -> inside: W[:,0] -= (dt/dx) * (flux_int[:,0] - flux_L)\n'
|
|
'\n'
|
|
'# Phase 5: tank also uses the SAME flux for update\n'
|
|
'fL_A = flux_L * pipe.area # flux x cross-section area = physical rate\n'
|
|
'tank1.apply_flux(mdot=fL_A[0], edot=fL_A[2], dt=dt, sign=-1)\n'
|
|
' -> inside: mass -= fL_A[0] * dt\n'
|
|
' U -= fL_A[2] * dt')
|
|
|
|
doc.add_heading('6.3 Why Does This Guarantee Conservation?', level=2)
|
|
doc.add_paragraph(
|
|
'Consider mass conservation at the left boundary. In one time step dt:\n'
|
|
'- Pipe cell 0 gains mass from left interface flux = +flux_L[0] * A * dt\n'
|
|
'- Tank 1 loses mass from the same flux = flux_L[0] * A * dt\n'
|
|
'\n'
|
|
'The two cancel exactly. The same applies to the right boundary and the energy component. '
|
|
'Therefore, the total system quantity (Tank 1 + all Pipe cells + Tank 2) is strictly '
|
|
'conserved at every time step, with errors only from floating-point roundoff, '
|
|
'typically at the 10^(-15) level.'
|
|
)
|
|
|
|
doc.add_paragraph(
|
|
'Verification result: After simulation, the relative errors of total mass and total '
|
|
'energy are ~6.5 x 10^(-16) (mass) and ~1.1 x 10^(-15) (energy), i.e., '
|
|
'machine-epsilon level.'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 7. TIME STEPPING
|
|
# ========================================================================
|
|
doc.add_heading('7 Time-Stepping Algorithm', level=1)
|
|
|
|
doc.add_heading('7.1 Explicit Euler Time Integration', level=2)
|
|
doc.add_paragraph(
|
|
'This system uses first-order explicit Euler time integration. For each pipe cell:'
|
|
)
|
|
add_equation(doc, 'W_i^(n+1) = W_i^(n) - (dt/dx) * [F_{i+1/2}^(n) - F_{i-1/2}^(n)]')
|
|
doc.add_paragraph(
|
|
'For the tanks:'
|
|
)
|
|
add_equation(doc, 'mass^(n+1) = mass^(n) + sign * m_dot * dt')
|
|
add_equation(doc, 'U^(n+1) = U^(n) + sign * E_dot * dt')
|
|
|
|
doc.add_heading('7.2 CFL Condition', level=2)
|
|
doc.add_paragraph(
|
|
'Explicit time integration has a stability restriction: the time step cannot exceed '
|
|
'the time for a signal to traverse one grid cell. '
|
|
'The CFL (Courant-Friedrichs-Lewy) condition is:'
|
|
)
|
|
add_equation(doc, 'dt = CFL * dx / max_i(|u_i| + a_i)')
|
|
doc.add_paragraph(
|
|
'Where CFL is a safety factor in (0, 1] (this system uses CFL = 0.5), '
|
|
'and max_i(|u_i| + a_i) is the maximum characteristic speed across all pipe cells. '
|
|
'Additionally, dt is clamped to ensure it does not exceed the remaining simulation time '
|
|
't_end - t.'
|
|
)
|
|
|
|
doc.add_heading('7.3 Complete 7-Phase Time Step', level=2)
|
|
doc.add_paragraph(
|
|
'Each time step follows this execution flow (corresponding to the run() function in solver.py):'
|
|
)
|
|
|
|
make_table(doc,
|
|
['Phase', 'Operation', 'Description'],
|
|
[
|
|
['Phase 1', 'CFL dt calculation', 'dt = CFL*dx/max(|u|+a), clamped to t_end-t'],
|
|
['Phase 2', 'Freeze ghost cells', 'W_ghost_L = tank1.ghost_state()\nW_ghost_R = tank2.ghost_state()'],
|
|
['Phase 3', 'Compute boundary HLL fluxes', 'flux_L = hll(ghost_L, pipe[:,0])\nflux_R = hll(pipe[:,-1], ghost_R)'],
|
|
['Phase 4', 'Advance pipe', 'pipe.step(flux_L, flux_R, dt)\nincludes internal HLL fluxes + friction source'],
|
|
['Phase 5', 'Advance tanks (flux doubling)', 'tank1.apply_flux(flux_L*A, sign=-1)\ntank2.apply_flux(flux_R*A, sign=+1)'],
|
|
['Phase 6', 'Advance time', 't += dt, step += 1'],
|
|
['Phase 7', 'Record history', 'Store P1, T1, P2, T2, pipe.W'],
|
|
],
|
|
col_widths=[2, 4, 9],
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 8. PIPE STEP DETAIL
|
|
# ========================================================================
|
|
doc.add_heading('8 Pipe Single-Step Update Detail (pipe.step)', level=1)
|
|
|
|
doc.add_heading('8.1 Flux Computation and State Update', level=2)
|
|
doc.add_paragraph(
|
|
'The pipe.step(flux_L, flux_R, dt) method is the core of the 1D pipe model. '
|
|
'It receives the two boundary fluxes from the solver, internally computes N-1 '
|
|
'interior interface fluxes, and then updates all cells using the finite-volume formula. '
|
|
'The detailed steps:'
|
|
)
|
|
|
|
doc.add_paragraph(
|
|
'(1) Snapshot current state: W_snap = W.copy() (ensures same time level)\n'
|
|
'\n'
|
|
'(2) Compute N-1 interior interface fluxes:\n'
|
|
' For k = 0, 1, ..., N-2:\n'
|
|
' flux_int[:, k] = hll_flux(W_snap[:, k], W_snap[:, k+1], gamma)\n'
|
|
'\n'
|
|
'(3) Update first cell (Cell 0):\n'
|
|
' Left face = flux_L (from Tank 1 ghost cell)\n'
|
|
' Right face = flux_int[:, 0]\n'
|
|
' W[:, 0] = W_snap[:, 0] - (dt/dx) * (flux_int[:, 0] - flux_L)\n'
|
|
'\n'
|
|
'(4) Update interior cells (Cell 1 to Cell N-2):\n'
|
|
' W[:, i] = W_snap[:, i] - (dt/dx) * (flux_int[:, i] - flux_int[:, i-1])\n'
|
|
'\n'
|
|
'(5) Update last cell (Cell N-1):\n'
|
|
' Left face = flux_int[:, N-2]\n'
|
|
' Right face = flux_R (from Tank 2 ghost cell)\n'
|
|
' W[:, N-1] = W_snap[:, N-1] - (dt/dx) * (flux_R - flux_int[:, N-2])'
|
|
)
|
|
|
|
doc.add_heading('8.2 Spatial Relationship of Fluxes and Cells', level=2)
|
|
add_code_block(doc,
|
|
' flux_L flux_int[0] flux_int[1] flux_int[N-2] flux_R\n'
|
|
' | | | ... | |\n'
|
|
' v v v v v\n'
|
|
' | Cell 0 | Cell 1 | Cell 2 | ... | Cell N-1 |\n'
|
|
' | W[:,0] | W[:,1] | W[:,2] | | W[:,N-1] |\n'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 9. FRICTION
|
|
# ========================================================================
|
|
doc.add_heading('9 Wall Friction: Source Term Method', level=1)
|
|
|
|
doc.add_heading('9.1 Modified Governing Equations', level=2)
|
|
doc.add_paragraph(
|
|
'With wall friction, the 1D Euler equations become a system with source terms:'
|
|
)
|
|
add_equation(doc, 'dW/dt + dF(W)/dx = S')
|
|
doc.add_paragraph('Where the source term vector S is:')
|
|
add_equation(doc, 'S = [0, -f/D * rho*u*|u|/2, 0]^T')
|
|
doc.add_paragraph(
|
|
'Meaning of each component:\n'
|
|
'- S[0] = 0: Friction does not create or destroy mass\n'
|
|
'- S[1] = -f/D * rho*u*|u|/2: Darcy-Weisbach wall friction force (per unit volume)\n'
|
|
' -- f is the Darcy friction factor (dimensionless)\n'
|
|
' -- D is the pipe inner diameter\n'
|
|
' -- |u| ensures the drag direction always opposes the flow direction\n'
|
|
'- S[2] = 0: Under the adiabatic wall assumption, kinetic energy dissipated by friction '
|
|
'is entirely converted to internal energy; total energy (internal + kinetic) remains unchanged\n'
|
|
'\n'
|
|
'Note: S[2] = 0 means wall friction does not change the system total energy. '
|
|
'Friction decelerates the flow (momentum decreases), but the lost kinetic energy '
|
|
'is converted to internal energy via frictional heating; their sum remains constant. '
|
|
'Therefore, even with friction enabled, the system total energy remains strictly conserved.'
|
|
)
|
|
|
|
doc.add_heading('9.2 Operator Splitting', level=2)
|
|
doc.add_paragraph(
|
|
'To maintain code clarity and modularity, friction source terms are handled via '
|
|
'operator splitting: within each time step, first complete the source-free flux update '
|
|
'(Euler equation part), then separately apply the source term. '
|
|
'This is equivalent to Lie Splitting:'
|
|
)
|
|
add_equation(doc, 'W* = W^(n) - (dt/dx)*[F_{i+1/2} - F_{i-1/2}] (flux step)')
|
|
add_equation(doc, 'W^(n+1)[1] = W*[1] + dt * S[1] (source step, momentum only)')
|
|
doc.add_paragraph(
|
|
'Since S[0] = S[2] = 0, the source step only updates the momentum component W[1] = rho*u. '
|
|
'Implementation:'
|
|
)
|
|
add_code_block(doc,
|
|
'if self.mu > 0:\n'
|
|
' rho_s = self.W[0, :]\n'
|
|
' u_s = self.W[1, :] / rho_s\n'
|
|
' abs_u = np.abs(u_s)\n'
|
|
' Re = rho_s * abs_u * self.D / self.mu\n'
|
|
' f = darcy_friction_factor(Re, self.eps_D)\n'
|
|
' S_mom = -f / self.D * rho_s * u_s * abs_u / 2.0\n'
|
|
' self.W[1, :] += dt * S_mom')
|
|
|
|
doc.add_heading('9.3 Darcy Friction Factor Calculation', level=2)
|
|
doc.add_paragraph(
|
|
'The friction factor f depends on the Reynolds number Re and relative wall roughness eps/D:'
|
|
)
|
|
add_equation(doc, 'Re = rho * |u| * D / mu')
|
|
doc.add_paragraph(
|
|
'Where mu is the dynamic viscosity. Depending on Re, three regimes are used:'
|
|
)
|
|
|
|
make_table(doc,
|
|
['Flow Regime', 'Re Range', 'Friction Factor Formula'],
|
|
[
|
|
['Laminar', 'Re < 2300', 'f = 64 / Re'],
|
|
['Transition', '2300 <= Re <= 4000', 'Linear blend of laminar and turbulent:\n'
|
|
'f = (1-alpha)*f_lam + alpha*f_turb\nalpha = (Re - 2300) / 1700'],
|
|
['Turbulent', 'Re > 4000', 'Colebrook-White implicit equation:\n'
|
|
'1/sqrt(f) = -2*log10(eps/(3.7*D) + 2.51/(Re*sqrt(f)))\n'
|
|
'Solved via Swamee-Jain initial guess + 10 fixed-point iterations'],
|
|
],
|
|
col_widths=[2.5, 4, 9],
|
|
)
|
|
|
|
doc.add_paragraph(
|
|
'The current default configuration uses smooth pipe walls (ROUGHNESS = 0), '
|
|
'in which case the Colebrook-White equation simplifies to the smooth-pipe implicit friction law.'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 10. CONSERVATION PROOF
|
|
# ========================================================================
|
|
doc.add_heading('10 Conservation Analysis and Verification', level=1)
|
|
|
|
doc.add_heading('10.1 Total System Mass', level=2)
|
|
add_equation(doc, 'M_total = m_tank1 + m_tank2 + SUM_i(rho_i * A * dx)')
|
|
doc.add_paragraph(
|
|
'Where the summation runs over all N pipe cells. Due to the flux doubling mechanism, '
|
|
'within each time step:\n'
|
|
'- Tank 1 mass loss = flux_L[0] * A * dt\n'
|
|
'- Cell 0 mass gain = (flux_L[0] - flux_int[0]) * (A*dt/dx) * dx = difference\n'
|
|
'- ... internal cell flux differences exactly "pass" mass through ...\n'
|
|
'- Cell N-1 mass gain = (flux_int[N-2] - flux_R[0]) * ...\n'
|
|
'- Tank 2 mass gain = flux_R[0] * A * dt\n'
|
|
'\n'
|
|
'This is a Telescoping Sum: all intermediate terms cancel, and the boundary terms '
|
|
'are exactly offset by the tank updates. Therefore M_total is strictly invariant.'
|
|
)
|
|
|
|
doc.add_heading('10.2 Total System Energy', level=2)
|
|
add_equation(doc, 'E_total = U_tank1 + U_tank2 + SUM_i((rho*E)_i * A * dx)')
|
|
doc.add_paragraph(
|
|
'The analysis is entirely analogous to mass conservation. '
|
|
'Note: friction source term S[2] = 0, so friction does not affect energy conservation.'
|
|
)
|
|
|
|
doc.add_heading('10.3 Numerical Verification Results', level=2)
|
|
doc.add_paragraph(
|
|
'Conservation checks are performed in both main.py and tests/test_integration.py:'
|
|
)
|
|
make_table(doc,
|
|
['Conserved Quantity', 'Initial Value', 'Final Value', 'Relative Error'],
|
|
[
|
|
['Total mass', '5.226485 x 10^2 kg', '5.226485 x 10^2 kg', '~6.5 x 10^(-16)'],
|
|
['Total energy', '1.125001 x 10^8 J', '1.125001 x 10^8 J', '~1.1 x 10^(-15)'],
|
|
],
|
|
)
|
|
doc.add_paragraph(
|
|
'Relative errors are at the 10^(-15) level, i.e., IEEE 754 double-precision '
|
|
'machine epsilon, confirming the correctness of the flux doubling mechanism.'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 11. SUMMARY OF DATA FLOW
|
|
# ========================================================================
|
|
doc.add_heading('11 Data Flow Summary', level=1)
|
|
|
|
doc.add_paragraph(
|
|
'The following summarizes the data flow between modules in one complete time step:'
|
|
)
|
|
add_code_block(doc,
|
|
'+-----------------------------------------------------------+\n'
|
|
'| solver.run() |\n'
|
|
'| |\n'
|
|
'| (1) CFL: a_max = pipe.max_wave_speed() |\n'
|
|
'| dt = CFL * dx / a_max |\n'
|
|
'| |\n'
|
|
'| (2) Ghost: W_gL = tank1.ghost_state() |\n'
|
|
'| W_gR = tank2.ghost_state() |\n'
|
|
'| |\n'
|
|
'| (3) Boundary flux: |\n'
|
|
'| fL = hll_flux(W_gL, pipe.W[:,0]) <-- riemann.py |\n'
|
|
'| fR = hll_flux(pipe.W[:,-1], W_gR) <-- riemann.py |\n'
|
|
'| |\n'
|
|
'| (4) Pipe update: |\n'
|
|
'| pipe.step(fL, fR, dt) |\n'
|
|
'| +-- internal HLL fluxes (riemann.py) |\n'
|
|
'| +-- finite-volume update W |\n'
|
|
'| +-- friction source term (friction.py) |\n'
|
|
'| |\n'
|
|
'| (5) Tank update (same fL, fR): |\n'
|
|
'| tank1.apply_flux(fL*A, sign=-1) |\n'
|
|
'| tank2.apply_flux(fR*A, sign=+1) |\n'
|
|
'| |\n'
|
|
'| (6)(7) t += dt, record history |\n'
|
|
'+-----------------------------------------------------------+\n'
|
|
)
|
|
|
|
doc.add_page_break()
|
|
|
|
# ========================================================================
|
|
# 12. PARAMETERS
|
|
# ========================================================================
|
|
doc.add_heading('12 Simulation Parameter Summary', level=1)
|
|
|
|
doc.add_heading('12.1 Gas Properties', level=2)
|
|
make_table(doc,
|
|
['Parameter', 'Symbol', 'Value', 'Unit'],
|
|
[
|
|
['Heat capacity ratio', 'gamma', '1.4', '--'],
|
|
['Gas constant', 'R', '287.0', 'J/(kg*K)'],
|
|
['Dynamic viscosity', 'mu', '1.8 x 10^(-5)', 'Pa*s'],
|
|
],
|
|
)
|
|
|
|
doc.add_heading('12.2 Geometry and Initial Conditions', level=2)
|
|
make_table(doc,
|
|
['Parameter', 'Symbol', 'Value', 'Unit'],
|
|
[
|
|
['High-pressure tank volume', 'V1', '5.0', 'm^3'],
|
|
['High-pressure tank initial pressure', 'P1', '5.0', 'MPa'],
|
|
['High-pressure tank initial temperature', 'T1', '300.0', 'K'],
|
|
['Low-pressure tank volume', 'V2', '10.0', 'm^3'],
|
|
['Low-pressure tank initial pressure', 'P2', '2.0', 'MPa'],
|
|
['Low-pressure tank initial temperature', 'T2', '300.0', 'K'],
|
|
['Pipe length', 'L', '1.0', 'm'],
|
|
['Pipe inner diameter', 'D', '5.0', 'mm'],
|
|
['Wall roughness', 'eps', '0.0', 'm (smooth pipe)'],
|
|
],
|
|
)
|
|
|
|
doc.add_heading('12.3 Numerical Control', level=2)
|
|
make_table(doc,
|
|
['Parameter', 'Symbol', 'Value', 'Description'],
|
|
[
|
|
['Number of cells', 'N', '20', 'Pipe uniformly divided into 20 finite-volume cells'],
|
|
['Cell size', 'dx', '50 mm', '= L / N'],
|
|
['CFL number', 'CFL', '0.5', 'Time-step safety factor'],
|
|
['Simulation end time', 't_end', '0.1', 'seconds'],
|
|
],
|
|
)
|
|
|
|
# ========================================================================
|
|
# SAVE
|
|
# ========================================================================
|
|
out_path = os.path.join('docs', 'pipe_system_simulation_technical_doc.docx')
|
|
os.makedirs('docs', exist_ok=True)
|
|
doc.save(out_path)
|
|
print(f'Document saved to: {out_path}')
|
|
return out_path
|
|
|
|
|
|
if __name__ == '__main__':
|
|
build_document()
|