"""Reproducible virtual pressure/flow and dosing bench; no physical observations.""" from pathlib import Path import csv, hashlib, html, itertools, json, math import numpy as np from flow import solve_flow from support import dump, csvfile, table ROOT = Path(__file__).resolve().parent def resistance(z, r, mu): """Lubrication resistance, Pa/(mm3/s); integral of 8 mu/(pi r**4) dz.""" if not math.isfinite(mu) or mu <= 0 or len(z) != len(r) or len(z)<2: raise ValueError('Positive finite viscosity and paired stations required') if any(not math.isfinite(x) for x in [*z,*r]) or any(x<=0 for x in r) or any(b<=a for a,b in zip(z,z[1:])): raise ValueError('Finite positive radii and increasing stations required') total=0.0 for x,y,a,b in zip(z,z[1:],r,r[1:]): g=(1/a**3-1/b**3)/(3*(b-a)) if abs(b-a)>1e-12 else 1/a**4 total+=(y-x)*g return 8*mu/math.pi*total def dose(r_nozzle, r_feed, compliance, supply, capillary, on_time, off_time): """Exact solution for a linear compliant node, first filled at capillary pressure. During on: C dP/dt=(Ps-P)/Rs-(P-Pcap)/Rn. During off the inlet closes. Returned pulse + tail satisfy the integrated flow balance. """ vals=[r_nozzle,r_feed,compliance,supply,capillary,on_time,off_time] if not all(math.isfinite(v) for v in vals) or min(r_nozzle,r_feed,on_time)<=0 or min(compliance,capillary,off_time)<0 or supply<=capillary: raise ValueError('Positive resistance/time/head and nonnegative compliance required') q=(supply-capillary)/(r_nozzle+r_feed) if compliance==0:return {'steady_flow_mm3_s':q,'pulse_mm3':q*on_time,'tail_mm3':0.,'stored_at_close_mm3':0.,'remaining_mm3':0.,'inlet_volume_mm3':q*on_time,'tau_on_s':0.,'tau_off_s':0.,'closure_pressure_Pa':capillary+q*r_nozzle,'balance_residual_mm3':0.} tau=compliance*r_nozzle*r_feed/(r_nozzle+r_feed) rise=-math.expm1(-on_time/tau) pulse=q*(on_time-tau*rise) pressure=capillary+q*r_nozzle*rise stored=compliance*(pressure-capillary) tau_off=compliance*r_nozzle tail=stored*(-math.expm1(-off_time/tau_off)) remaining=stored*math.exp(-off_time/tau_off) # Exact integral of inlet: initial flow (Ps-Pcap)/Rs, settling to q. inlet=q*on_time+q*r_nozzle/r_feed*tau*rise residual=inlet-pulse-tail-remaining return {'steady_flow_mm3_s':q,'pulse_mm3':pulse,'tail_mm3':tail,'stored_at_close_mm3':stored,'remaining_mm3':remaining,'inlet_volume_mm3':inlet,'tau_on_s':tau,'tau_off_s':tau_off,'closure_pressure_Pa':pressure,'balance_residual_mm3':residual} def checks_record(checks, id, description, expected, observed, passed, method): checks.append(dict(id=id,description=description,expected=expected,observed=observed,passed=bool(passed),method=method)) 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) c=next(c for c in cfg['candidates'] if c['id']=='wide-inlet');z=c['z_mm'];r=c['r_mm'];tol=v['dimensional_tolerances_mm'];checks=[] rec=lambda *args:checks_record(checks,*args) pipe=resistance([0,20],[1,1],.75) rec('resistance_pipe','Constant-radius reduction','matches 8 mu L/(pi r^4) to 1e-12',pipe,abs(pipe-8*.75*20/math.pi)<1e-12,'Independent closed-form benchmark') # FEM runs at manufacturing corners, not linear scaling from nominal geometry. corners=[] for name,si,so,sl in [('minimum_resistance',1,1,-1),('nominal',0,0,0),('maximum_resistance',-1,-1,1)]: zz=[0,z[1],z[2]-sl*tol['outlet_land_length'],z[3]] rr=[r[0]+si*tol['inlet_radius']]*2+[r[2]+so*tol['outlet_radius']]*2 f=solve_flow(zz,rr,20,.75,nz=120,nr=24);pd=f.diagnostics()['pressure_drop_Pa'];approx=resistance(zz,rr,.75)*20 error=abs(approx/pd-1) corners.append(dict(case=name,inlet_radius_mm=rr[0],outlet_radius_mm=rr[-1],outlet_land_mm=zz[-1]-zz[-2],fem_pressure_Pa=pd,reduced_pressure_Pa=approx,relative_error=error)) rec('fem_'+name,'Reduced resistance checked at '+name,'relative pressure error < 0.01',error,error<.01,'New axisymmetric P2/P1 Stokes solve at 120x24 mesh') csvfile(out/'geometry-corners.csv',corners) # Use a documented upper correction, rounded upward, to bound the checked approximation. correction=max(x['fem_pressure_Pa']/x['reduced_pressure_Pa'] for x in corners)*1.005 envelope=[] for si,so,sl,mu,q in itertools.product([-1,0,1],[-1,0,1],[-1,0,1],cfg['fluid']['viscosities_Pa_s'],cfg['flow_rates_mm3_s']): zz=[0,z[1],z[2]-sl*tol['outlet_land_length'],z[3]];rr=[r[0]+si*tol['inlet_radius']]*2+[r[2]+so*tol['outlet_radius']]*2 rn=resistance(zz,rr,mu)*correction rt=resistance([0,v['supply_tube']['length_mm']],[v['supply_tube']['radius_mm']]*2,mu) pcap=2*v['surface_tension_N_m']*v['contact_angle_cosine_bound']/(rr[-1]*1e-3) pressure=q*(rn+rt)+pcap qmax=(cfg['pressure_budget_Pa']-pcap)/(rn+rt) re=cfg['fluid']['density_kg_m3']*q/(math.pi*rr[-1]**2)*(2*rr[-1])*1e-6/mu envelope.append(dict(inlet_radius_mm=rr[0],outlet_radius_mm=rr[-1],land_mm=zz[-1]-zz[-2],viscosity_Pa_s=mu,flow_mm3_s=q,pressure_Pa=pressure,pressure_margin_Pa=cfg['pressure_budget_Pa']-pressure,max_flow_for_budget_mm3_s=qmax,outlet_Re=re,pressure_pass=pressure<=cfg['pressure_budget_Pa'])) csvfile(out/'tolerance-pressure-envelope.csv',envelope) worst=max(envelope,key=lambda x:x['pressure_Pa']);safe_flow=math.floor(min(x['max_flow_for_budget_mm3_s'] for x in envelope)*10)/10 bounded_pressure=max((x['pressure_Pa']-2*v['surface_tension_N_m']/(x['outlet_radius_mm']*1e-3))/x['flow_mm3_s']*safe_flow+2*v['surface_tension_N_m']/(x['outlet_radius_mm']*1e-3) for x in envelope) rec('adverse_pressure_detected','Upper range including tolerances, feed and meniscus','at least one 30 mm3/s case exceeds unchanged 6500 Pa budget',worst['pressure_Pa'],worst['pressure_Pa']>cfg['pressure_budget_Pa'],'Full factorial reduced-order tolerance study; checked FEM correction') rec('guardband','Conservative flow bound from the same model','pressure <= 6500 Pa for every enumerated case',bounded_pressure,bounded_pressure<=cfg['pressure_budget_Pa'],'Worst permissible flow rounded down to 0.1 mm3/s; original target retained') rec('creeping_regime','Flow regime over tolerance envelope','max Re < 0.1',max(x['outlet_Re'] for x in envelope),max(x['outlet_Re'] for x in envelope)<.1,'Outlet diameter and mean velocity') # A synthetic bench explicitly creates noise/offset; this is an inference test, not experimental validation. rng=np.random.default_rng(v['seed']);nominal_R=corners[1]['fem_pressure_Pa']/20;measurements=[];fits=[] for label,factor in [('clean',1.),('obstructed',1.25)]: rows=[] for q in cfg['flow_rates_mm3_s']: for repeat in range(v['collection_repeats']): measured_q=q+rng.normal(0,v['balance_noise_sigma_mg'])/v['collection_time_s'] # rho=1 mg/mm3 p=nominal_R*factor*q+v['pressure_sensor']['common_offset_Pa']+rng.normal(0,v['pressure_sensor']['independent_noise_sigma_Pa']) row={'case':label,'flow_command_mm3_s':q,'repeat':repeat+1,'flow_measured_mm3_s':measured_q,'pressure_measured_Pa':p,'synthetic':True};rows.append(row);measurements.append(row) slope,offset=np.polyfit([x['flow_measured_mm3_s'] for x in rows],[x['pressure_measured_Pa'] for x in rows],1) fits.append({'case':label,'known_resistance_Pa_s_mm3':nominal_R*factor,'fitted_resistance_Pa_s_mm3':float(slope),'fitted_offset_Pa':float(offset),'relative_resistance_error':float(abs(slope/(nominal_R*factor)-1))}) csvfile(out/'synthetic-pressure-observations.csv',measurements);csvfile(out/'synthetic-fit-results.csv',fits) rec('synthetic_fit','Recover a declared resistance despite pressure offset/noise','relative slope error < 1% in both cases',max(x['relative_resistance_error'] for x in fits),max(x['relative_resistance_error'] for x in fits) 20%',increase,increase>.2,'25% resistance challenge compared with clean run') # Fixed pressure/pulse, not a volume-controlled pump. Capillary pressure has fixed positive sign for conservative head. dc=v['dose'];rf=resistance([0,v['supply_tube']['length_mm']],[v['supply_tube']['radius_mm']]*2,.75);cap=2*v['surface_tension_N_m']/(r[-1]*1e-3) supply=dc['nominal_flow_mm3_s']*(nominal_R+rf)+cap;doses=[] for name,air,viscosity_factor in [('deaired',0,1),('small_air_pocket',25,1),('large_air_pocket',250,1),('high_viscosity_same_setting',0,1/.75)]: compliance=dc['liquid_volume_mm3']/dc['bulk_modulus_Pa']+air/dc['ambient_pressure_Pa'] result=dose(nominal_R*viscosity_factor,rf*viscosity_factor,compliance,supply,cap,dc['on_time_s'],2) doses.append({'case':name,'air_volume_mm3':air,'viscosity_factor':viscosity_factor,'compliance_mm3_Pa':compliance,**result,'pulse_error_mm3':result['pulse_mm3']-dc['target_mm3'],'pulse_pass':abs(result['pulse_mm3']-dc['target_mm3'])<=dc['tolerance_mm3'],'tail_pass':result['tail_mm3']<=1}) csvfile(out/'dose-and-shutoff.csv',doses) rec('volume_balance','Compliant dosing volume balance','absolute residual < 1e-10 mm3',max(abs(x['balance_residual_mm3']) for x in doses),max(abs(x['balance_residual_mm3']) for x in doses)<1e-10,'Closed-form inlet = pulse + tail + stored remainder') zero=dose(200,10,0,4200,0,5,2) rec('zero_compliance','Zero-compliance limiting case','100 mm3 delivered during pulse and zero tail',zero['pulse_mm3'],abs(zero['pulse_mm3']-100)<1e-12 and zero['tail_mm3']==0,'Analytical steady resistance network') rec('adverse_air','Large trapped air affects shutoff','250 mm3 air case exceeds 1 mm3 tail criterion',doses[2]['tail_mm3'],doses[2]['tail_mm3']>1,'Linear isothermal gas-compliance sensitivity; no bubble transport claim') rec('adverse_viscosity','Fixed pulse is not a volume controller','1 Pa s case misses 100 +/- 1 mm3 target',doses[3]['pulse_mm3'],not doses[3]['pulse_pass'],'Same regulator setting and pulse length; changed viscosity') invalid=0 for args in [([0,1],[1,0],1),([0,0],[1,1],1),([0,1],[1,1],float('nan'))]: try:resistance(*args) except ValueError:invalid+=1 rec('input_rejection','Invalid resistance inputs','3 invalid cases rejected',invalid,invalid==3,'Zero bore, duplicate station and nonfinite viscosity') acceptance=[{'requirement':'6500 Pa at 30 mm3/s, tolerance + supply + meniscus','outcome':'FAIL','evidence':f"Worst {worst['pressure_Pa']:.2f} Pa; original nominal-only maximum was 6142.62 Pa."},{'requirement':'Computed conservative operating range','outcome':'PASS','evidence':f"Up to {safe_flow:.1f} mm3/s across this bounded input set; maximum {bounded_pressure:.2f} Pa."},{'requirement':'Nominal fixed pulse 100 +/- 1 mm3','outcome':'PASS' if doses[0]['pulse_pass'] else 'FAIL','evidence':f"De-aired pulse {doses[0]['pulse_mm3']:.6f} mm3."},{'requirement':'Shutoff after large air pocket','outcome':'FAIL','evidence':f"{doses[2]['tail_mm3']:.3f} mm3 tail in 2 s exceeds 1 mm3."},{'requirement':'Unchanged setting after viscosity increase','outcome':'FAIL','evidence':f"{doses[3]['pulse_mm3']:.3f} mm3 during pulse."}] csvfile(out/'acceptance-matrix.csv',acceptance) decision={'geometry':'wide-inlet CAD, unchanged','maximum_commanded_flow_mm3_s':safe_flow,'viscosity_Pa_s':[.5,1.0],'pressure_budget_Pa':6500,'worst_computed_pressure_Pa':bounded_pressure,'operating_rule':f'Use the unchanged selected CAD within the declared tolerance, feed-tube and meniscus bounds, with commanded flow no greater than the computed {safe_flow:.1f} mm3/s bound. The former 30 mm3/s endpoint is rejected for this combined system envelope. Fixed-pressure pulse dosing is qualified only for the stated nominal de-aired case; changed viscosity requires recalculating pressure or pulse duration.'} dump(out/'selected-operating-envelope.json',decision) result={'selected_operating_envelope':decision,'synthetic':True,'pressure_envelope_cases':len(envelope),'fem_geometry_checks':corners,'reduced_model_upper_factor':correction,'worst_pressure_Pa':worst['pressure_Pa'],'safe_flow_mm3_s':safe_flow,'safe_flow_worst_pressure_Pa':bounded_pressure,'synthetic_observations':len(measurements),'dose_cases':doses,'acceptance':acceptance,'checks':checks} dump(out/'results.json',result) import matplotlib.pyplot as plt from support import COLORS fig,axes=plt.subplots(1,2,figsize=(11,4.5),layout='constrained') for label,color in [('clean',COLORS[0]),('obstructed',COLORS[1])]: rows=[x for x in measurements if x['case']==label];axes[0].scatter([x['flow_measured_mm3_s'] for x in rows],[x['pressure_measured_Pa'] for x in rows],label=label,s=22,color=color) axes[0].set(xlabel='Synthetic collected flow (mm3/s)',ylabel='Synthetic pressure reading (Pa)',title='Known offset and noise; 15 repeats per state');axes[0].legend() axes[1].bar(range(4),[x['tail_mm3'] for x in doses],color=COLORS);axes[1].axhline(1,color=COLORS[2],linestyle='--');axes[1].set_xticks(range(4),['De-aired','25 mm3 air','250 mm3 air','Higher viscosity'],rotation=20);axes[1].set(ylabel='Calculated tail in 2 s (mm3)',title='Same pressure setting and 5 s pulse') for a in axes:a.grid(axis='y',alpha=.3) fig.savefig(out/'virtual-bench.png',dpi=160);plt.close(fig) methods='''

This virtual bench completes the pressure/flow, dimensional sensitivity and dosing run-through using declared inputs. No instrument data were collected. Every generated observation is marked synthetic.

Equations and units

The slow-taper resistance is R = (8 mu / pi) integral dz / r(z)^4, in Pa s/mm3 when geometry is in mm. A linear radius segment integrates to L (r0^-3 - r1^-3) / [3 (r1-r0)]. Three new FEM solutions check this reduction at nominal and tolerance corners. The maximum observed FEM/reduced ratio receives a further 0.5% allowance; it is a local empirical bound, not a general error theorem.

Upstream straight-tube pressure is Q Rfeed. The bounded meniscus head is 2 gamma / rout with radius in metres. This is a static spherical-meniscus sensitivity with contact-angle cosine bounded by 1; it does not simulate droplets. The tolerance envelope uses independent corner/nominal radii, land length, viscosity and flow: 243 cases. It retains the original 6500 Pa pressure budget. The output derives a lower permitted flow instead of silently moving the budget.

For dosing, C dP/dt = (Ps-P)/Rfeed - (P-Pcap)/Rnozzle while open. The inlet is isolated at closure and C dP/dt = -(P-Pcap)/Rnozzle afterward. C = Vliquid/Kbulk + Vair/Pabsolute. Exact exponentials give pulse volume, tail and remaining stored liquid. The linear gas-compliance approximation is used only for the few-kPa pressure swing; it models a compressible pocket, not moving bubbles. Pressure is initially Pcap. The inlet valve is assumed instantaneous. Tail liquid leaves the nozzle; droplet detachment is not computed.

The synthetic bench imposes a common 22 Pa offset and 12 Pa independent normal noise on pressure, plus 0.2 mg balance noise over 20 s at the declared 1 mg/mm3 density. Fitting an intercept prevents the imposed zero offset from being mistaken for resistance. A seeded 25% obstruction challenge tests detection. This verifies the inference workflow against known generated truth; it is not independent physical validation of the flow solver.

''' refs='

Method references: Ansys fully developed laminar-flow verification; DOLFINx Taylor-Hood formulation. The resistance-network and compliance equations above are derived directly from flow and volume balances.

' page='Nozzle virtual bench

Dispensing nozzle: virtual bench

'+methods+'

Selected operating decision

The existing wide-inlet CAD is retained. The combined system envelope is limited to '+str(safe_flow)+' mm3/s at 0.5 to 1.0 Pa s with the declared tolerances, feed tube and meniscus bound. Its worst computed pressure is '+f'{bounded_pressure:.2f}'+' Pa. The 30 mm3/s endpoint is rejected for this envelope. No dimensions were changed and no new CAD revision is implied.

For the separate 100 mm3 pulse example, use the nominal 0.75 Pa s de-aired condition. A change of viscosity or gas compliance requires a new regulator/pulse calculation; the failed fixed-setting cases are not approved operating points.

Acceptance results

'+table(acceptance)+'

FEM comparison at tolerance corners

'+table(corners)+'

Dosing and shutoff

'+table(doses)+'

Verification

'+table(checks)+'Synthetic pressure observations and calculated shutoff tails

Model boundary

The study does not select a commercial pump, material, seal or thread. Surface finish, cleaning chemistry, wear and non-Newtonian behavior cannot be inferred from these declared properties. Physical qualification remains a different evidence category. The completed deliverable is the reproducible computational assessment and its operating envelope.

'+refs+'

Machine-readable results | Inputs | Run instructions

' (out/'report.html').write_text(page) sources=['virtual_bench.py','virtual-bench-inputs.json','inputs.json','flow.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':['Synthetic measurements exercise analysis and fault detection; no physical liquid or pump was tested.','Tolerance correction is checked at three declared geometry corners, not a proof for arbitrary tapers.','Compliance is a lumped linear model; no droplet breakup, moving bubbles, non-Newtonian fluid or surface chemistry is resolved.'],'summary':f"{len(checks)} verification checks; 243 tolerance/operating cases; 30 synthetic observations; four dosing cases. Original maximum-flow condition fails added system effects; computed bound {safe_flow:.1f} mm3/s."} dump('completion.json',completion) delivery=json.loads(Path('delivery.json').read_text());delivery['completion']=completion['summary'];delivery['status']='CFD and synthetic bench completed';delivery['selected_operating_envelope']=decision for item in [{'path':'outputs/virtual-bench/selected-operating-envelope.json','label':'Selected operating envelope'},{'path':'completion.json','label':'Completed virtual-bench verification'},{'path':'virtual_bench.py','label':'Virtual pressure and dosing bench source'},{'path':'virtual-bench-inputs.json','label':'Tolerance, instrument and dosing inputs'},{'path':'outputs/virtual-bench/report.html','label':'Virtual bench and acceptance report'},{'path':'outputs/virtual-bench/acceptance-matrix.csv','label':'Acceptance outcomes'},{'path':'outputs/virtual-bench/tolerance-pressure-envelope.csv','label':'All 243 pressure-envelope cases'}]: delivery['key_files']=[x for x in delivery['key_files'] if x['path']!=item['path']]+[item] delivery['assumptions']=[x for x in delivery['assumptions'] if not x.startswith('Combined system envelope:')]+[f'Combined system envelope: command at most {safe_flow:.1f} mm3/s at 0.5 to 1.0 Pa s with +/-0.01 mm bore radii, +/-0.1 mm outlet-land length, the declared 10 mm feed tube and bounded static meniscus. The original 30 mm3/s endpoint is comparison evidence, not an accepted combined-system operating point.'] delivery['implementation']=[x for x in delivery['implementation'] if not x.startswith('The virtual bench executes')]+['The virtual bench executes tolerance-corner FEM checks, pressure-envelope calculations, synthetic instrument fits, obstruction detection and a compliant dosing/shutoff model. The final operating decision is supplied as JSON.'] delivery['verification']['command']='python run.py' base=delivery['verification']['summary'].split(' Additional virtual bench:')[0] delivery['verification']['summary']=base+' Additional virtual bench: '+completion['summary'] delivery['results']=[{**x,'label':'Nozzle-only pressure at original 30 mm3/s endpoint'} if x['label']=='Maximum flow / viscosity pressure' else x for x in delivery['results'] if x['label']!='Virtual-bench flow bound']+[{'label':'Virtual-bench flow bound','value':safe_flow,'unit':'mm3/s'}] delivery['limitations']=[x for x in delivery['limitations'] if not x.startswith('No surface tension')] delivery['limitations'].append('Additional static meniscus and lumped compliance models cover pressure-head and shutoff sensitivities. They do not resolve droplet breakup, moving bubbles or physical contact angle.') if not any(x.startswith('Additional static meniscus') for x in delivery['limitations']) else None 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) supplement='

Completed virtual bench

'+completion['summary']+'

Read the synthetic acceptance results and methods. The supplement records adverse outcomes as well as passing cases.

' text=re.sub(r'()',r'\1'+supplement,text,count=1) root_report.write_text(text) print(json.dumps({'passed':completion['passed'],'checks':len(checks),'worst_pressure_Pa':worst['pressure_Pa'],'safe_flow_mm3_s':safe_flow,'dose_tails_mm3':[x['tail_mm3'] for x in doses]},indent=2)) if not completion['passed']:raise SystemExit(1) if __name__=='__main__':main()