"""Virtual flush bench: tolerances, inlet mixing, reversible wall retention and sensor lag.
The wall model is a conservative linear two-capacity exchange model, not fitted chemistry.
"""
from pathlib import Path
import hashlib, itertools, json, math
import numpy as np
from scipy.sparse import block_diag, diags
from scipy.sparse.linalg import factorized
from transport import assemble, crossing, run_transport
from flow import solve_flow
from model import export_chambers
from support import dump, csvfile, table
ROOT=Path(__file__).resolve().parent
def pressure(length,a,b,q,mu):
if not all(math.isfinite(x) for x in [length,a,b,q,mu]) or min(length,a,q,mu)<=0 or b<=a:
raise ValueError('Finite positive dimensions, flow, viscosity and clearance required')
return 8*mu*length*q/(math.pi*(b**4-a**4-(b*b-a*a)**2/math.log(b/a)))
def run_case(length,a,b,q,diffusion,region,*,nz=80,nr=24,dt=.25,capacity=0,exchange_time=30,inlet_volume=0,sensor_tau=0,end_time=900):
if min(capacity,inlet_volume,sensor_tau)<0 or exchange_time<=0 or dt<=0:raise ValueError('Nonnegative capacities/lags; positive exchange time and timestep required')
m=assemble(length,a,b,q,diffusion,nz,nr)
n=nz*nr;dz=length/nz;dr=(b-a)/nr
# Every axial strip has separate inner- and outer-wall capacity.
areas=np.r_[np.full(nz,2*math.pi*a*dz),np.full(nz,2*math.pi*b*dz)]
attached=np.r_[np.arange(nz)*nr,np.arange(nz)*nr+nr-1]
if capacity>0:
wc=capacity*areas;weights=np.r_[m.volumes,wc]
op=block_diag((m.operator,np.zeros((2*nz,2*nz))),format='lil')
for i,(k,cap) in enumerate(zip(attached,wc)):
w=n+i;g=cap/exchange_time
op[k,k]-=g;op[w,w]-=g;op[k,w]+=g;op[w,k]+=g
else:weights=m.volumes;op=m.operator
solve=factorized(diags(weights,format='csc')-dt*op.tocsc())
z_overlap=np.maximum(0,np.minimum(m.z+dz/2,region['z_max_mm'])-np.maximum(m.z-dz/2,region['z_min_mm']))
r_hi=np.minimum(m.r+dr/2,a+region['radial_depth_mm']);r_lo=np.maximum(m.r-dr/2,a)
probe=(z_overlap[:,None]*(math.pi*np.maximum(0,r_hi**2-r_lo**2))[None,:]).ravel();probe/=probe.sum()
c=np.ones(len(weights));initial=float(weights.sum());incoming=outgoing=0.;sensor=1.;rows=[];maximum_error=0.;minimum=1.;maximum=1.
def append(t,cin):
fluid=float(np.dot(m.volumes,c[:n])/m.initial_mass);near=float(np.dot(probe,c[:n]));outlet=float(np.dot(m.ring_flow,c[n-nr:n])/q)
wall=float(np.dot(weights[n:],c[n:])) if capacity else 0.
mass=float(np.dot(weights,c));err=(mass+outgoing-incoming-initial)/initial
rows.append({'time_s':t,'inlet_fraction':cin,'whole_fluid_fraction':fluid,'probe_fraction':near,'sensor_fraction':sensor,'outlet_fraction':outlet,'total_inventory_fraction':mass/initial,'wall_inventory_mm3':wall,'integrated_inlet_mm3':incoming,'integrated_outlet_mm3':outgoing,'mass_balance_fraction':err})
return err
append(0,1 if inlet_volume else 0)
tau_in=inlet_volume/q if inlet_volume else 0
for step in range(1,int(end_time/dt)+1):
t=step*dt;old=c;rhs=weights*c
# Exact step-average exponential inlet preserves its integrated added mass.
cin=tau_in/dt*(math.exp(-(t-dt)/tau_in)-math.exp(-t/tau_in)) if tau_in else 0.
rhs[:nr]+=dt*m.ring_flow*cin
c=solve(rhs)
near=float(np.dot(probe,c[:n]));previous_near=float(np.dot(probe,old[:n]))
# Sensor uses trapezoidal input with exact zero-order response on each small step.
sensor=near if sensor_tau==0 else sensor*math.exp(-dt/sensor_tau)+.5*(near+previous_near)*(-math.expm1(-dt/sensor_tau))
incoming+=dt*q*cin;outgoing+=dt*float(np.dot(m.ring_flow,c[n-nr:n]))
err=append(t,cin);maximum_error=max(maximum_error,abs(err));minimum=min(minimum,float(c.min()));maximum=max(maximum,float(c.max()))
if max(rows[-1]['total_inventory_fraction'],near,sensor)<.002:break
times=[x['time_s'] for x in rows]
metrics={f't99_{name}_s':crossing(times,[x[key] for x in rows]) for name,key in [('fluid','whole_fluid_fraction'),('probe','probe_fraction'),('sensor','sensor_fraction'),('inventory','total_inventory_fraction'),('outlet','outlet_fraction')]}
metrics.update(maximum_mass_balance_fraction=maximum_error,minimum_fraction=minimum,maximum_fraction=maximum,initial_inventory_mm3=initial,initial_wall_inventory_mm3=float(weights[n:].sum()) if capacity else 0,elapsed_s=rows[-1]['time_s'],nz=nz,nr=nr,dt_s=dt)
return rows,metrics
def main():
import os;os.chdir(ROOT)
cfg=json.loads(Path('inputs.json').read_text());v=json.loads(Path('virtual-bench-inputs.json').read_text());out=Path('outputs/virtual-bench');out.mkdir(parents=True,exist_ok=True)
checks=[]
def rec(id,description,expected,observed,passed,method):checks.append(dict(id=id,description=description,expected=expected,observed=observed,passed=bool(passed),method=method))
length=cfg['length_mm'];a=cfg['probe_radius_mm'];b=cfg['chamber_radii_mm'][0];q=cfg['flow_mm3_s'];mu=cfg['viscosity_Pa_s'];D=cfg['diffusivity_mm2_s'];region=cfg['probe_region'];robust=v['robust_chamber_radius_mm']
envelope=[]
for label,nom_b in [('original_nominal',b),('tolerance_allowance',robust)]:
for ai,bi,qi,mui in itertools.product([-1,0,1],[-1,0,1],v['flow_range_mm3_s'],v['viscosity_range_Pa_s']):
aa=a+ai*v['probe_radius_tolerance_mm'];bb=nom_b+bi*v['chamber_radius_tolerance_mm'];p=pressure(length,aa,bb,qi,mui)
re=cfg['density_kg_m3']*qi/(math.pi*(bb**2-aa**2))*2*(bb-aa)*1e-6/mui
envelope.append({'design':label,'probe_radius_mm':aa,'chamber_radius_mm':bb,'flow_mm3_s':qi,'viscosity_Pa_s':mui,'radial_clearance_mm':bb-aa,'pressure_Pa':p,'Re':re,'clearance_pass':bb-aa>=cfg['minimum_radial_clearance_mm']-1e-12,'pressure_pass':p<=cfg['pressure_budget_Pa']})
csvfile(out/'dimensional-envelope.csv',envelope)
old=[x for x in envelope if x['design']=='original_nominal'];new=[x for x in envelope if x['design']=='tolerance_allowance']
rec('original_clearance_failure','Apply declared manufacturing tolerances to original 0.5 mm nominal gap','minimum gap below 0.5 mm is detected',min(x['radial_clearance_mm'] for x in old),any(not x['clearance_pass'] for x in old),'Full factorial exact annular hydraulic model')
rec('allowance_clearance','2.525 mm nominal radius covers +/-0.015 body and +/-0.010 probe radii','all 81 cases have >= 0.5 mm clearance',min(x['radial_clearance_mm'] for x in new),all(x['clearance_pass'] for x in new),'Worst independent radial tolerance sum, no statistical cancellation')
rec('allowance_pressure','Unchanged 10 Pa pressure budget for revised nominal radius','all 81 cases below 10 Pa',max(x['pressure_Pa'] for x in new),all(x['pressure_pass'] for x in new),'Exact developed annular Stokes solution')
rec('regime','Tolerance/flow/viscosity range remains creeping flow','max Re < 0.1',max(x['Re'] for x in envelope),max(x['Re'] for x in envelope)<.1,'Hydraulic diameter 2(b-a)')
worst=max(new,key=lambda x:x['pressure_Pa']);f=solve_flow([0,length],[worst['chamber_radius_mm']]*2,worst['flow_mm3_s'],worst['viscosity_Pa_s'],inner=worst['probe_radius_mm'],nz=80,nr=24)
pd=f.diagnostics()['pressure_drop_Pa'];err=abs(pd/worst['pressure_Pa']-1)
rec('fem_worst_corner','Check revised worst hydraulic case against independent FEM','relative error < 1e-5',err,err<1e-5,'P2/P1 Stokes solve at revised tolerance corner')
cad=export_chambers(length,a,[robust]);dump(out/'cad-checks.json',cad)
rec('cad_roundtrip','Revised nominal assembly exports and reimports','valid body/probe solids; exact fluid volume',cad,all(x['body_valid'] and x['probe_valid'] for x in cad),'CadQuery solid validity and STEP roundtrip')
# Response cases deliberately retain original selected fluid geometry, so existing results can be compared directly.
cases=[('ideal_passive',0,30,0,0),('sensor_5s',0,30,0,5),('sensor_30s',0,30,0,30),('inlet_mix_20mm3',0,30,20,5),('wall_fast',.02,30,0,5),('wall_slow',.1,120,20,30)]
summaries=[];all_traces={}
for name,capacity,tau,vin,ts in cases:
rows,metrics=run_case(length,a,b,q,D,region,capacity=capacity,exchange_time=tau,inlet_volume=vin,sensor_tau=ts)
all_traces[name]=rows;csvfile(out/(name+'-trace.csv'),rows)
summaries.append({'case':name,'capacity_length_mm':capacity,'wall_exchange_s':tau,'inlet_mixed_volume_mm3':vin,'sensor_tau_s':ts,**metrics,'sensor_within_90s':metrics['t99_sensor_s'] is not None and metrics['t99_sensor_s']<=v['response_limit_s']})
rec('mass_'+name,'Fluid + reversible wall + inlet/outlet conservation in '+name,'max normalized residual < 1e-10',metrics['maximum_mass_balance_fraction'],metrics['maximum_mass_balance_fraction']<1e-10,'Sparse mass-form backward Euler; exchange fluxes cancel in coupled capacities')
rec('bounded_'+name,'Concentration bounds in '+name,'all fluid/wall concentrations in [0,1] +/- 1e-10',[metrics['minimum_fraction'],metrics['maximum_fraction']],metrics['minimum_fraction']>=-1e-10 and metrics['maximum_fraction']<=1+1e-10,'Every solved state checked')
# Recover the already implemented conservative transport solver when wall, inlet lag and sensor lag vanish.
m=assemble(length,a,b,q,D,80,24);base=run_transport(m,.25,{**region,'inner_radius_mm':a});base_t=summaries[0]['t99_probe_s'];reldiff=abs(base_t/base['t99_probe_s']-1)
rec('zero_retention_limit','New coupled model reduces to existing passive model','relative t99 probe error < 1e-10',reldiff,reldiff<1e-10,'Independent public run_transport execution at same mesh and time step')
# Re-run the difficult wall/inlet/sensor case with both grid and timestep refinement.
refinements=[]
for nz,nr,dt in [(80,24,.25),(160,48,.25),(160,48,.125)]:
if nz==80:met={k:summaries[-1][k] for k in ['t99_fluid_s','t99_probe_s','t99_sensor_s','t99_inventory_s','maximum_mass_balance_fraction']}
else:
rows,met=run_case(length,a,b,q,D,region,nz=nz,nr=nr,dt=dt,capacity=.1,exchange_time=120,inlet_volume=20,sensor_tau=30)
if dt==.125:csvfile(out/'wall-slow-refined-trace.csv',rows)
refinements.append({'nz':nz,'nr':nr,'dt_s':dt,**{k:met[k] for k in ['t99_fluid_s','t99_probe_s','t99_sensor_s','t99_inventory_s','maximum_mass_balance_fraction']}})
csvfile(out/'response-refinement.csv',refinements)
for label,first,last in [('grid',refinements[0],refinements[1]),('timestep',refinements[1],refinements[2])]:
change=max(abs(last[k]/first[k]-1) for k in ['t99_fluid_s','t99_probe_s','t99_sensor_s','t99_inventory_s'])
rec('refine_'+label,'Slow-wall case '+label+' sensitivity','max t99 change < 0.03',change,change<.03,'80x24 to 160x48 and 0.25 s to 0.125 s, original threshold retained')
# A separately solved isolated two-capacity exchange benchmark.
vf,vw,g=10.,2.,.04;dt=.01;duration=20.;A=np.array([[-g/vf,g/vf],[g/vw,-g/vw]]);c=np.array([1.,0.]);M=np.linalg.inv(np.eye(2)-dt*A)
for _ in range(round(duration/dt)):c=M@c
equilibrium=vf/(vf+vw);exact=np.array([equilibrium+vw/(vf+vw)*math.exp(-g*(1/vf+1/vw)*duration),equilibrium-equilibrium*math.exp(-g*(1/vf+1/vw)*duration)])
exchange_error=float(np.max(np.abs(c-exact)))
rec('exchange_benchmark','Isolated fluid/wall exchange approaches conserved equilibrium','max absolute error < 5e-5',exchange_error,exchange_error<5e-5,'Backward Euler compared with independent exact two-capacity exponential')
mass=float(vf*c[0]+vw*c[1]);rec('exchange_mass','Closed two-capacity benchmark conserves initial tracer','absolute difference < 1e-10',abs(mass-vf),abs(mass-vf)<1e-10,'Independent dense update and analytical conserved mass')
# Pressure instrument demonstration: differential elevation can exceed the chamber pressure margin.
rng=np.random.default_rng(v['seed']);inst=v['pressure_instrument'];true_p=pressure(length,a,b,q,mu);hydro=cfg['density_kg_m3']*9.80665*inst['tap_elevation_error_mm']*1e-3
readings=[]
for i in range(20):
raw=true_p+hydro+inst['zero_offset_Pa']+rng.normal(0,inst['noise_sigma_Pa'])
readings.append({'repeat':i+1,'raw_Pa':float(raw),'corrected_Pa':float(raw-hydro-inst['zero_offset_Pa']),'true_synthetic_Pa':true_p,'synthetic':True})
csvfile(out/'synthetic-pressure-observations.csv',readings)
mean=float(np.mean([x['corrected_Pa'] for x in readings]));rec('pressure_correction','Recover known synthetic hydraulic pressure after elevation/zero correction','absolute mean error < 0.08 Pa',mean-true_p,abs(mean-true_p)<.08,'20 seeded readings, known imposed elevation and zero bias')
late=summaries[2]['t99_sensor_s'];rec('sensor_lag_challenge','30 s sensor time constant changes the reference conclusion','sensor t99 exceeds 90 s while passive probe t99 is below 90 s',late,late>90 and summaries[0]['t99_probe_s']<90,'Causal first-order sensor driven by the spatial probe average')
wall=summaries[-1];rec('retention_challenge','Slow wall inventory is retained in response accounting','inventory t99 exceeds passive-fluid t99',wall['t99_inventory_s'],wall['t99_inventory_s']>summaries[0]['t99_fluid_s'],'Reversible exchange with explicitly assumed capacity and time constant')
invalid=0
for args in [(20,2,2,5,.01),(20,2,2.5,-1,.01),(20,2,2.5,5,float('nan'))]:
try:pressure(*args)
except ValueError:invalid+=1
rec('invalid_inputs','Reject invalid annular hydraulic conditions','3 invalid conditions rejected',invalid,invalid==3,'Zero gap, negative flow, nonfinite viscosity')
# Adopt the tolerance-allowance geometry and explicitly qualify a longer flush decision.
# This is a finite synthetic qualification set, not an inference about an unspecified real sensor.
selected_runs=[]
for ai,bi in itertools.product([-1,1],[-1,1]):
aa=a+ai*v['probe_radius_tolerance_mm'];bb=robust+bi*v['chamber_radius_tolerance_mm']
rows,met=run_case(length,aa,bb,4.5,D,region,capacity=.1,exchange_time=120,inlet_volume=20,sensor_tau=30,end_time=900)
at600=next(x for x in rows if abs(x['time_s']-600)<1e-8)
selected_runs.append({'probe_radius_mm':aa,'chamber_radius_mm':bb,'flow_mm3_s':4.5,'capacity_length_mm':.1,'wall_exchange_s':120,'inlet_volume_mm3':20,'sensor_tau_s':30,'check_time_s':600,'sensor_fraction':at600['sensor_fraction'],'fluid_fraction':at600['whole_fluid_fraction'],'inventory_fraction':at600['total_inventory_fraction'],'outlet_fraction':at600['outlet_fraction'],'t99_inventory_s':met['t99_inventory_s'],'t99_sensor_s':met['t99_sensor_s'],'maximum_mass_balance_fraction':met['maximum_mass_balance_fraction']})
selected_worst=max(selected_runs,key=lambda x:x['inventory_fraction'])
fine_rows,fine_met=run_case(length,selected_worst['probe_radius_mm'],selected_worst['chamber_radius_mm'],4.5,D,region,nz=160,nr=48,dt=.125,capacity=.1,exchange_time=120,inlet_volume=20,sensor_tau=30,end_time=900)
fine_at600=next(x for x in fine_rows if abs(x['time_s']-600)<1e-8)
csvfile(out/'selected-flush-corners.csv',selected_runs);csvfile(out/'selected-flush-fine-trace.csv',fine_rows)
max_fraction=max(max(x['sensor_fraction'],x['fluid_fraction'],x['inventory_fraction'],x['outlet_fraction']) for x in selected_runs)
fine_fraction=max(fine_at600[k] for k in ['sensor_fraction','whole_fluid_fraction','total_inventory_fraction','outlet_fraction'])
rec('selected_flush_decision','Adopted 600 s flush on tolerance-allowance geometry','all fluid, sensor, inventory and outlet fractions <= 0.01 at four geometry corners',max_fraction,max_fraction<=.01,'Four new 80x24 runs at 4.5 mm3/s, 0.1 mm wall capacity, 120 s exchange, 20 mm3 inlet and 30 s sensor')
rec('selected_flush_refined','Verify the slowest selected corner at finer resolution','all four fractions <= 0.01 at 600 s',fine_fraction,fine_fraction<=.01,'Additional 160x48 run at 0.125 s; full trace retained')
decision={'nominal_chamber_radius_mm':robust,'nominal_probe_radius_mm':a,'minimum_gap_mm':min(x['radial_clearance_mm'] for x in new),'maximum_hydraulic_pressure_Pa':worst['pressure_Pa'],'flush_time_s':600,'qualified_flush_flow_mm3_s':4.5,'wall_capacity_length_mm':.1,'wall_exchange_s':120,'inlet_mixed_volume_mm3':20,'sensor_time_constant_s':30,'diffusivity_mm2_s':D,'maximum_old_fraction_at_flush_end':max(max_fraction,fine_fraction),'qualification':'Four independent geometry-tolerance corners plus finer slowest-corner solve. Parameters are explicit synthetic qualification values, not an assertion about unmeasured sensors. The original 90 s desired response is rejected for the slow-retention condition; 600 s is the selected operating wait for this virtual fixture.'}
dump(out/'selected-operating-envelope.json',decision)
csvfile(out/'response-summary.csv',summaries)
acceptance=[{'requirement':'Minimum 0.5 mm physical radial clearance with declared tolerances','outcome':'ORIGINAL FAIL; REVISED PASS','evidence':f"Original minimum {min(x['radial_clearance_mm'] for x in old):.3f} mm. Revised 2.525 mm nominal chamber radius gives minimum {min(x['radial_clearance_mm'] for x in new):.3f} mm."},{'requirement':'10 Pa hydraulic budget with revised dimensions and input range','outcome':'PASS','evidence':f"Worst {worst['pressure_Pa']:.4f} Pa."},{'requirement':'90 s indicated t99 reference criterion','outcome':'CASE DEPENDENT','evidence':f"Passive probe {summaries[0]['t99_sensor_s']:.2f} s; 30 s sensor lag {late:.2f} s; slow wall/inlet/sensor case {wall['t99_sensor_s']:.2f} s."},{'requirement':'Fluid plus adsorbed inventory conserved','outcome':'PASS','evidence':f"Maximum residual {max(x['maximum_mass_balance_fraction'] for x in summaries):.3g}."},{'requirement':'Measured pressure equals fluid resistance without corrections','outcome':'FAIL','evidence':f"Declared 0.2 mm tap-height error alone adds {hydro:.4f} Pa."}]
csvfile(out/'acceptance-matrix.csv',acceptance)
acceptance.append({'requirement':'Selected 600 s flush on revised geometry at the declared slow-retention condition','outcome':'PASS' if max(max_fraction,fine_fraction)<=.01 else 'FAIL','evidence':f"Maximum old fraction {max(max_fraction,fine_fraction):.6f}; four tolerance corners and finer worst-corner solve."})
csvfile(out/'acceptance-matrix.csv',acceptance)
dump(out/'results.json',{'selected_operating_envelope':decision,'synthetic':True,'tolerance_cases':len(envelope),'reference_response_limit_s':v['response_limit_s'],'response_cases':summaries,'refinement':refinements,'acceptance':acceptance,'checks':checks})
import matplotlib.pyplot as plt
from support import COLORS
fig,axs=plt.subplots(1,2,figsize=(11,4.5),layout='constrained')
for j,name in enumerate(['ideal_passive','sensor_30s','wall_slow']):
rows=all_traces[name];axs[0].semilogy([x['time_s'] for x in rows],[max(x['sensor_fraction'],1e-5) for x in rows],label=name.replace('_',' '),color=COLORS[j]);axs[1].semilogy([x['time_s'] for x in rows],[max(x['total_inventory_fraction'],1e-5) for x in rows],label=name.replace('_',' '),color=COLORS[j])
for ax in axs:ax.axhline(.01,color=COLORS[3],linestyle='--');ax.set(xlabel='Time after switch (s)',ylim=(.001,1.1));ax.grid(alpha=.3);ax.legend()
axs[0].set(ylabel='Indicated old-sample fraction',title='Sensor response includes declared lag');axs[1].set(ylabel='Remaining total tracer / initial tracer',title='Fluid plus reversible wall capacity')
fig.savefig(out/'virtual-flush.png',dpi=160);plt.close(fig)
fig,ax=plt.subplots(figsize=(10,4.4),layout='constrained')
ax.fill_between([0,length],-a,a,color=COLORS[0],alpha=.7)
ax.fill_between([0,length],robust,robust+2,color=COLORS[1],alpha=.8)
ax.fill_between([0,length],-robust-2,-robust,color=COLORS[1],alpha=.8)
ax.plot([0,length],[0,0],color=COLORS[2],linestyle='--',linewidth=.7)
ax.annotate('Chamber bore 5.050 +/- 0.030 mm diameter',(5,robust),xytext=(1,5.7),arrowprops={'arrowstyle':'-','color':COLORS[2]},fontsize=10)
ax.annotate('Probe 4.000 +/- 0.020 mm diameter',(11,a),xytext=(10,6.5),arrowprops={'arrowstyle':'-','color':COLORS[2]},fontsize=10)
ax.annotate('Nominal radial gap 0.525 mm; worst 0.500 mm',(15,-(a+robust)/2),xytext=(1,-6.5),arrowprops={'arrowstyle':'-','color':COLORS[2]},fontsize=10)
ax.set(xlabel='Axial position (mm)',ylabel='Radial coordinate (mm)',xlim=(-1,21),ylim=(-7.2,7.6),title='Selected tolerance-allowance reference: 20 mm active length, 2 mm nominal wall')
ax.grid(alpha=.25);fig.savefig(out/'selected-section.svg');fig.savefig(out/'selected-section.png',dpi=160);plt.close(fig)
methods='''This is an executed virtual flushing bench. All sensor constants, wall parameters, dimensional tolerances and instrument errors are synthetic reference inputs. No physical chamber, tracer or sensor was operated.
Coupled transport and wall model
The existing finite-volume mass form is extended with two wall states per axial strip. A wall capacity W = K A has units of equivalent liquid volume. Transfer J = (W / tau) (cfluid - cwall) removes tracer from the adjacent fluid cell and adds exactly the same amount to that wall state. Equilibrium is cwall = cfluid. This linear reversible retention model describes sensitivity to extra stored tracer; K and tau are declared parameters, not chemistry fitted to a material. Both walls start at concentration 1. The reported inventory includes fluid and wall mass.
The inlet may pass through a perfectly mixed 20 mm3 volume. Its old-sample fraction decays as exp(-Qt/Vin); each numerical step uses the exact average of that function. The measured probe signal follows tau_sensor dy/dt = cprobe - y. It has an exact exponential decay factor and trapezoidal forcing over each step. The zero-lag case reads the probe average directly. Inlet tracer, outlet tracer and retained inventory are recorded at every time step.
The baseline and response variations use the original 2.5 mm chamber radius to isolate wall and instrument effects. A separate tolerance study evaluates the 2.525 mm nominal-radius variant and exports its CAD. This revised geometry is a dimensional proposal supported by hydraulic calculations, not a fabricated part. Its transport response is not substituted into the original baseline curves.
Hydraulic and dimensional calculation
For a straight concentric annulus, deltaP = 8 mu L Q / [pi (b^4 - a^4 - (b^2-a^2)^2/ln(b/a))]. The 162 cases cover both original and allowance variants, independently varying probe/body radii, flow and viscosity. A new FEM solve checks the worst hydraulic corner of the revised geometry. The 10 Pa budget and minimum 0.5 mm clearance remain unchanged.
The 90 s indicated-response limit is an explicit new virtual-bench criterion, not a requirement from a customer. Failure is reported when a slow sensor or retained wall inventory causes that limit to be exceeded. This prevents a fluid-only result from being presented as sensor performance.
Instrument check
Differential pressure contains hydrostatic head rho g delta_h and instrument zero bias. Twenty readings include the declared 0.2 mm elevation error, 0.2 Pa zero offset and 0.08 Pa independent Gaussian noise. Correcting the known injected errors tests the reduction workflow. It does not establish the calibration of a real low-pressure transducer.
'''
page='Sampling chamber virtual benchSampling chamber: virtual flush bench
'+methods+'Selected design and operating decision
Adopt the supplied 2.525 mm chamber-radius CAD with the stated radial tolerances. For the declared slow-retention virtual fixture, use a 600 s flush at 4.5 mm3/s before accepting a reading. Four geometry corners and a finer worst-corner solve put fluid, outlet, sensor and total inventory below 1% at that time. Maximum remaining fraction is '+f'{max(max_fraction,fine_fraction):.6f}'+'.
The 90 s response objective is rejected for that slow-retention condition. The selected wait does not claim that objective was met. Wall and sensor values remain explicit synthetic inputs; a different measured response requires recalculating the wait.
'+table(selected_runs)+'Acceptance results
'+table(acceptance)+'Response cases
'+table(summaries)+'Refinement
'+table(refinements)+'Verification
'+table(checks)+'
Selected reference dimensions
Model boundary
The model remains axisymmetric and straight. Wall retention is linear, inlet dead volume is perfectly mixed, and the sensor is first order. It does not resolve real side ports, trapped gas geometry, nonlinear binding, seals or leakage. These are completed theoretical sensitivity cases with exposed assumptions, not a replacement for physical characterization.
Method context: COMSOL transport and adsorption example describes coupled fluid and wall species balances. This package uses its own explicitly derived linear two-capacity model, not the COMSOL example chemistry or software.
Results JSON | Inputs | Run instructions
'
(out/'report.html').write_text(page)
sources=['virtual_bench.py','virtual-bench-inputs.json','inputs.json','flow.py','transport.py','model.py','run.py']
completion={'schema_version':1,'evidence_type':'synthetic/theoretical','command':'python virtual_bench.py','passed':all(x['passed'] for x in checks),'checks':checks,'source_sha256':{p:hashlib.sha256(Path(p).read_bytes()).hexdigest() for p in sources},'report':'outputs/virtual-bench/report.html','limitations':['No physical tracer, sensor or material adsorption measurements.','Wall capacity and rate are declared linear sensitivities, not fitted chemical properties.','The straight annulus omits asymmetric port geometry, bubbles, leakage and nonlinear chemistry.'],'summary':f"{len(checks)} checks; 162 dimensional/operating cases; six coupled response cases; three slow-wall refinement runs; four selected-design corners plus refined worst corner; 20 synthetic instrument readings."}
dump('completion.json',completion)
delivery=json.loads(Path('delivery.json').read_text());delivery['completion']=completion['summary'];delivery['status']='Flow and virtual flushing study completed';delivery['selected_operating_envelope']=decision
files=[{'path':'outputs/virtual-bench/selected-section.svg','label':'Selected reference section and tolerances'},{'path':'outputs/virtual-bench/selected-operating-envelope.json','label':'Selected design and 600 s flush decision'},{'path':'completion.json','label':'Completed virtual-bench verification'},{'path':'virtual_bench.py','label':'Wall, inlet and sensor sensitivity source'},{'path':'virtual-bench-inputs.json','label':'Synthetic tolerances and instrument inputs'},{'path':'outputs/virtual-bench/report.html','label':'Virtual flushing acceptance report'},{'path':'outputs/virtual-bench/acceptance-matrix.csv','label':'Acceptance outcomes'},{'path':'outputs/virtual-bench/response-summary.csv','label':'Six response cases'},{'path':'outputs/cad/radius-2.525mm-assembly.step','label':'Selected tolerance-allowance assembly STEP'}]
for item in files:delivery['key_files']=[x for x in delivery['key_files'] if x['path']!=item['path']]+[item]
for item in delivery['key_files']:
if item['path']=='outputs/cad/radius-2.5mm-assembly.step':item['label']='Original 2.5 mm-radius baseline assembly STEP'
delivery['summary']='Annular flow, tracer transport and a virtual flushing bench. The selected reference adds a tolerance allowance and a verified 600 s wait for the declared slow-retention case.'
delivery['figures']=[x for x in delivery['figures'] if x['path']!='outputs/virtual-bench/selected-section.png']+[{'path':'outputs/virtual-bench/selected-section.png','caption':'Selected reference dimensions and tolerance allowance. The original three comparison geometries remain in the package.'}]
delivery['assumptions']=[x for x in delivery['assumptions'] if not x.startswith('Final selected variant:')]+['Final selected variant: 5.050 +/- 0.030 mm chamber bore and 4.000 +/- 0.020 mm probe diameter. The declared slow-retention virtual fixture uses a 600 s flush at 4.5 mm3/s, 0.1 mm equivalent wall-capacity length, 120 s exchange, 20 mm3 inlet mixing and 30 s sensor lag.']
delivery['implementation']=[x for x in delivery['implementation'] if not x.startswith('The virtual bench couples')]+['The virtual bench couples conservative fluid transport to reversible wall capacity, a mixed inlet volume and sensor lag. Four tolerance corners and a refined worst-corner run verify the selected 600 s flushing decision; the revised assembly and section drawing are included.']
delivery['verification']['command']='python run.py'
delivery['verification']['summary']=delivery['verification']['summary'].split(' Additional virtual bench:')[0]+' Additional virtual bench: '+completion['summary']
delivery['results']=[{'label':'Selected nominal radial clearance','value':robust-a,'unit':'mm'},{'label':'Minimum gap after declared tolerances','value':round(min(x['radial_clearance_mm'] for x in new),3),'unit':'mm'},{'label':'Worst hydraulic pressure in selected envelope','value':round(worst['pressure_Pa'],4),'unit':'Pa'},{'label':'Selected virtual-fixture flush time','value':600,'unit':'s'},{'label':'Maximum old fraction after selected flush','value':round(max(max_fraction,fine_fraction),6)},{'label':'Original passive-fluid probe t99','value':39.58,'unit':'s'}]
delivery['limitations']=[x for x in delivery['limitations'] if not x.startswith('Uniform annular passages only.')]
if not any(x.startswith('The virtual bench adds') for x in delivery['limitations']):delivery['limitations'].append('The virtual bench adds linear wall retention, mixed inlet volume and first-order sensor response. They are declared sensitivities, not measured chemistry. Asymmetric fittings and moving bubbles remain outside the resolved geometry.')
dump('delivery.json',delivery)
verification=json.loads(Path('verification.json').read_text())
verification['synthetic_completion']={'command':completion['command'],'passed':completion['passed'],'checks':len(checks),'report':completion['report']}
dump('verification.json',verification)
root_report=Path('report.html')
text=root_report.read_text()
import re
text=re.sub(r'', '', text, flags=re.S)
text=text.replace('The original comparison below is retained as baseline evidence. The final selection is the 2.525 mm-radius tolerance-allowance CAD and the 600 s virtual-fixture flush decision in the supplement.
','')
supplement=''
text=re.sub(r'()',r'\1'+supplement+'The original comparison below is retained as baseline evidence. The final selection is the 2.525 mm-radius tolerance-allowance CAD and the 600 s virtual-fixture flush decision in the supplement.
',text,count=1)
root_report.write_text(text)
print(json.dumps({'passed':completion['passed'],'checks':len(checks),'minimum_revised_clearance_mm':min(x['radial_clearance_mm'] for x in new),'maximum_revised_pressure_Pa':worst['pressure_Pa'],'slow_sensor_t99_s':late,'slow_wall_t99_s':wall['t99_sensor_s']},indent=2))
if not completion['passed']:raise SystemExit(1)
if __name__=='__main__':main()