"""Support/geometry uncertainty, torsion screening, and synthetic fixture calibration."""
import json,csv,math
import numpy as np
import matplotlib;matplotlib.use('Agg')
import matplotlib.pyplot as plt
from completion_util import ROOT,OUT,csv_file,check,table,finish
def fit_compliance(rows):
x=np.array([1/float(r['support_N_m']) for r in rows]);y=np.array([float(r['load_displacement_mm']) for r in rows]);A=np.column_stack([np.ones(len(x)),x]);coef=np.linalg.lstsq(A,y,rcond=None)[0];return coef,float(np.max(abs(A@coef-y)))
def predicted(coef,force,E_GPa,thickness_mm,k_MN_m):
if min(force,E_GPa,thickness_mm,k_MN_m)<=0:raise ValueError('Positive finite mechanical inputs required')
return (coef[0]*69/E_GPa*20/thickness_mm+coef[1]/(k_MN_m*1e6))*force/500
def rectangle_J(a,b):
a,b=max(a,b),min(a,b)
if b<=0:raise ValueError('Positive dimensions required')
return a*b**3/3*(1-.63*b/a+.052*(b/a)**5)
def torsion(force,offset_mm,rotation_Nm_rad,E_GPa=69,nu=.33,thickness_mm=20,center_mm=120):
# Two elastic rotational load paths, including actual lengths of solid/window sections.
G=E_GPa*1000/(2*(1+nu));Jsolid=rectangle_J(40,thickness_mm);Jwindow=2*rectangle_J(thickness_mm,12)
def compliance(lo,hi):
window=sum(max(0.,min(hi,b)-max(lo,a)) for a,b in [(40,100),(140,200)])
return (window/Jwindow+(hi-lo-window)/Jsolid)/G
left=compliance(20,center_mm);right=compliance(center_mm,220);spring=1/(rotation_Nm_rad*1000)
torque=force*offset_mm;structural_angle=torque*left*right/(left+right)
total_angle=torque*(left+spring)*(right+spring)/(left+right+2*spring)
return {'offset_z_mm':offset_mm,'rotational_support_Nm_rad':rotation_Nm_rad,'structural_twist_rad':structural_angle,'support_twist_rad':total_angle-structural_angle,'added_loadpoint_displacement_mm':offset_mm*total_angle}
def main():
c=json.loads((ROOT/'completion-inputs.json').read_text());g=json.loads((ROOT/'inputs.json').read_text());rows=list(csv.DictReader((ROOT/'outputs/load-cases.csv').open()));selected=[r for r in rows if r['profile']==c['profile']];fits={};calibration=[];checks=[]
for center in [.12,.15]:
rr=[r for r in selected if float(r['load_center_m'])==center];coef,err=fit_compliance(rr);fits[center]=coef
for r in rr:calibration.append({'center_m':center,'support_MN_m':float(r['support_N_m'])/1e6,'FEA_displacement_mm':float(r['load_displacement_mm']),'reduced_model_mm':float(coef[0]+coef[1]/float(r['support_N_m'])),'absolute_error_mm':float(abs(coef[0]+coef[1]/float(r['support_N_m'])-float(r['load_displacement_mm'])))})
checks.append(check(f'surrogate-{center}','Support-compliance surrogate matches executed FEA at all three supports','max error < 0.0001 mm',err,err<.0001,'Least-squares fit delta = elastic term + support coefficient / stiffness'))
csv_file('response-model-benchmark.csv',calibration)
rng=np.random.default_rng(c['seed']);sample=[];summaries=[]
for center,coef in fits.items():
values={name:rng.uniform(*c[key],c['samples_per_load_position']) for name,key in [('load_N','load_range_N'),('E_GPa','E_range_GPa'),('thickness_mm','thickness_range_mm'),('support_MN_m','support_range_MN_m')]}
delta=(coef[0]*69/values['E_GPa']*20/values['thickness_mm']+coef[1]/(values['support_MN_m']*1e6))*values['load_N']/500
for i in range(len(delta)):sample.append({'center_m':center,'draw':i+1,**{k:float(v[i]) for k,v in values.items()},'loadpoint_displacement_mm':float(delta[i]),'meets_0p05mm':bool(delta[i]<=.05)})
worst=predicted(coef,max(c['load_range_N']),min(c['E_range_GPa']),min(c['thickness_range_mm']),min(c['support_range_MN_m']))
denom=.05/(max(c['load_range_N'])/500)-coef[0]*69/min(c['E_range_GPa'])*20/min(c['thickness_range_mm']);required=float(coef[1]/denom/1e6) if denom>0 else None
summaries.append({'load_center_m':center,'synthetic_draws':len(delta),'sample_min_mm':float(min(delta)),'sample_max_mm':float(max(delta)),'sample_exceedances':int(sum(delta>.05)),'worst_corner_mm':worst,'minimum_support_for_corner_MN_m':required})
csv_file('synthetic-envelope.csv',sample);csv_file('envelope-summary.csv',summaries)
checks.append(check('uncertainty-count','Two load positions each have 10000 synthetic combinations','20000 rows',len(sample),len(sample)==20000,'Deterministic seeded uniform input sweep; no claim of field probability'))
checks.append(check('nominal','Selected nominal centered case meets deflection requirement','< 0.05 mm',predicted(fits[.12],500,69,20,20),predicted(fits[.12],500,69,20,20)<.05,'Reduced model anchored to converged FEA'))
checks.append(check('envelope-rejection','Combined compliant/thin/high-load corner is identified','at least one corner > 0.05 mm',[r['worst_corner_mm'] for r in summaries],any(r['worst_corner_mm']>.05 for r in summaries),'Bounded monotone worst-case corner, separate from random sampling'))
# Increasing the support design target to the rounded-up derived requirement closes the declared envelope.
design_k=math.ceil(max(r['minimum_support_for_corner_MN_m'] for r in summaries));closed=[predicted(coef,550,66,19.9,design_k) for coef in fits.values()]
checks.append(check('support-requirement','Derived minimum support target closes translational envelope','all corner values <= 0.05 mm',{'support_MN_m':design_k,'corner_mm':closed},max(closed)<=.05,'Analytical inversion of calibrated compliance law'))
torsions=[];baseline=predicted(fits[.12],500,69,20,20)
for offset in c['out_of_plane_offsets_mm']:
for rot in c['rotational_support_Nm_rad']:
r=torsion(500,offset,rot);r['combined_screen_mm']=baseline+r['added_loadpoint_displacement_mm'];r['meets_target']=r['combined_screen_mm']<=.05;torsions.append(r)
csv_file('out-of-plane-screen.csv',torsions)
selected_envelope=[]
for center,coef in fits.items():
t=torsion(550,2,500,E_GPa=66,thickness_mm=19.9,center_mm=center*1000)
value=predicted(coef,550,66,19.9,30)+t['added_loadpoint_displacement_mm']
selected_envelope.append({'center_m':center,'load_N':550,'E_GPa':66,'thickness_mm':19.9,'translation_support_MN_m':30,'rotation_support_Nm_rad':500,'offset_mm':2,'combined_displacement_mm':value,'including_fit_error_mm':value+.0001,'target_mm':.05,'accepted':value+.0001<=.05})
csv_file('selected-operating-envelope.csv',selected_envelope)
checks.append(check('selected-envelope','Final selected supports include translation, torsion and fit residual','all <= 0.05 mm at declared worst corner', [r['including_fit_error_mm'] for r in selected_envelope],all(r['accepted'] for r in selected_envelope),'30 MN/m translation and 500 N m/rad rotation per support; 550 N, 66 GPa, 19.9 mm, 2 mm offset'))
checks.append(check('zero-offset','Zero out-of-plane offset produces zero extra movement','0 mm',torsion(500,0,100)['added_loadpoint_displacement_mm'],torsion(500,0,100)['added_loadpoint_displacement_mm']==0,'Torque equilibrium and Saint-Venant rectangular-section approximation'))
weak=next(r for r in torsions if r['offset_z_mm']==2 and r['rotational_support_Nm_rad']==100)
strong=next(r for r in torsions if r['offset_z_mm']==2 and r['rotational_support_Nm_rad']==500)
checks.append(check('rotation-screen','Weak rotational mounts can invalidate the 2D-only result','100 N m/rad rejected; 500 N m/rad accepted at 2 mm offset',{'weak_mm':weak['combined_screen_mm'],'strong_mm':strong['combined_screen_mm']},not weak['meets_target'] and strong['meets_target'],'Two equal rotational springs plus integrated cross-section torsion compliance'))
records=[];slope=baseline/500
for cycle in range(c['cycles']):
last=0
for j,F in enumerate(c['load_cycle_N']):
direction=1 if F>=last else -1;noise=float(rng.normal(0,c['synthetic_sensor_noise_mm']));y=slope*F+direction*c['synthetic_hysteresis_halfwidth_mm']+noise
records.append({'cycle':cycle+1,'step':j,'load_N':F,'direction':direction,'ideal_fixture_mm':slope*F,'simulated_sensor_mm':y});last=F
csv_file('synthetic-load-cycles.csv',records);X=np.column_stack([np.ones(len(records)),[r['load_N'] for r in records]]);coef=np.linalg.lstsq(X,[r['simulated_sensor_mm'] for r in records],rcond=None)[0];error=abs(coef[1]/slope-1)
checks.append(check('synthetic-calibration','Regression recovers generating compliance from noisy load cycles','slope error < 5%',error,error<.05,'Ordinary least squares with intercept; three load/unload cycles'))
stress=max(float(r['stress_away_from_corners_MPa']) for r in selected);screens=[{'assumed_stress_factor':k,'screening_stress_MPa':k*stress,'allowable_MPa':c['nominal_yield_MPa']/c['elastic_screen_factor']} for k in c['stress_concentration_screens']];csv_file('stress-factor-screen.csv',screens)
checks.append(check('stress-screen','Finite stress multipliers remain below assumed yield/2','all <= 120 MPa',max(r['screening_stress_MPa'] for r in screens),all(r['screening_stress_MPa']<=r['allowable_MPa'] for r in screens),'Sensitivity only; factors are not resolved notch stresses'))
fig,axes=plt.subplots(1,2,figsize=(11,4));
for center,color in [(.12,'#ad6349'),(.15,'#436b55')]:axes[0].hist([r['loadpoint_displacement_mm'] for r in sample if r['center_m']==center],bins=45,alpha=.65,label=f'x = {center*1000:.0f} mm',color=color)
axes[0].axvline(.05,color='black',ls='--');axes[0].set(xlabel='Synthetic displacement (mm)',ylabel='Input combinations');axes[0].legend();axes[1].scatter([r['load_N'] for r in records],[r['simulated_sensor_mm'] for r in records],s=13,color='#66576d');axes[1].plot([0,500],[coef[0],coef[0]+coef[1]*500],color='#ad6349');axes[1].set(xlabel='Declared load (N)',ylabel='Synthetic sensor displacement (mm)');fig.tight_layout();fig.savefig(OUT/'commissioning.png',dpi=140);plt.close(fig)
limits=['The input draws are uniformly distributed scenarios, not measured manufacturing or installation statistics. Exceedance counts cannot be interpreted as the probability a real fixture will fail.','The response model is calibrated to 2D FEA. Matching its source data and synthetic calibration checks verify the workflow; they do not independently validate real hardware.','The torsion screen uses an approximate rectangular-section Saint-Venant constant and prescribed rotational end springs. It omits warping restraint, local load introduction, 3D contact and bolt preload.','Stress multipliers are declared sensitivity factors. They do not replace a filleted 3D notch model, fatigue assessment or a load rating.','The 0.05 mm target is met only inside the stated force, material, translation and rotation envelope. Mount stiffness values are requirements to be achieved, not measured fixture properties.']
sections=[('Selected operating envelope','The final reference selection is the 12 mm flange crosshead with at least 30 MN/m vertical stiffness and 500 N m/rad rotational stiffness per support. The modeled envelope is 0-550 N, modulus at least 66 GPa, thickness at least 19.9 mm, load at x = 120 or 150 mm, and through-thickness offset no greater than 2 mm. The combined calculation includes a 0.0001 mm response-fit allowance and remains below 0.05 mm. Larger offsets and weaker supports are rejected conditions.
'+table(selected_envelope)),('Decision',f'The selected 12 mm flange profile meets 0.05 mm at the nominal 500 N/20 MN/m support condition. The combined 550 N, 66 GPa, 19.9 mm, 16 MN/m corner does not. The derived translational support requirement rounds up to {design_k} MN/m per support for that declared envelope. A separate torsion screen shows why off-plane load placement and rotational mount stiffness also matter.
'),('Inputs and response model',table([{'parameter':k,'value':v} for k,v in c.items()])+'At each load position, fit delta = A + B/k to the three converged FEA support cases. Scale the elastic A term by force, inverse modulus and inverse thickness. Scale B/k by force only. This separates material/section compliance from the mount springs. The fit residual is retained below.
'+table(calibration)),('Bounded sensitivity',table(summaries)+'
The random sweep explores combinations; deterministic corners govern the stated support requirement. The required stiffness is obtained by solving the compliance equation for k, not by changing the displacement target.
'),('Out-of-plane screening','A through-thickness eccentricity e creates torque P e. Centered loads split torque equally; off-center loads split it according to each path compliance. Torsion compliance integrates the actual solid and windowed lengths on both sides of the load. Each load path includes its end rotational spring; torque divides inversely with their total compliances. The centered case reduces to equal half-torque paths. The window section uses two 20 x 12 mm rectangles; solid sections use 40 x 20 mm. The rectangular torsion approximation is J = a b cubed / 3 [1 - 0.63 b/a + 0.052 (b/a) to the fifth]. Mount rotation adds P e / (2 Krot). Extra load-point vertical movement is e times the total twist.
'+table(torsions)),('Synthetic instrument exercise',f'Three load/unload cycles contain {len(records)} observations, 0.0005 mm Gaussian sensor noise and a declared plus/minus 0.001 mm hysteresis offset. The recovered compliance differs from its generating value by {100*error:.2f}%. This tests data reduction; it is not an experimental FEA correlation.
'),('Stress sensitivity',table(screens)+'The peak reported finite stress excludes ideal corners. The factor sweep is a sensitivity screen; no sharp-corner stress convergence or fatigue life is implied.
')]
finish('Crosshead installation and measurement study',checks,sections,limits,['complete_validation.py','completion_util.py','completion-inputs.json','inputs.json','solver.py','outputs/load-cases.csv'],{'synthetic_combinations':len(sample),'derived_translational_only_requirement_MN_m':design_k,'selected_support_MN_m':30,'selected_rotational_support_Nm_rad':500,'selected_worst_with_fit_error_mm':max(r['including_fit_error_mm'] for r in selected_envelope),'calibration_observations':len(records),'synthetic_slope_error_fraction':error})
if __name__=='__main__':main()