ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

1D耦合双稳淬火仿真实现

1D耦合双稳淬火仿真实现 # TAB013-G-003 / TAB013-G-004 合并流水线 # 约束钉死8项合并校验规则全部落实 # 运行顺序smoke - smoke_meta写入 - 主扫描 - 锚点校验 - 合并归档 - 绘图 import numpy as np from scipy.integrate import solve_ivp from scipy.stats import linregress import json import time import matplotlib.pyplot as plt # GLOBAL PARAMETERS (全局固定smoke与主扫描同源) s_c 1.0 s_floor 1e-12 t_end_global 200.0 n_eval_global 2000 slope_tol_smoke 0.05 slope_tol_anchor 0.1 # phase_label 整数编码锁定 PHASE_LOCKED 0 # 冻结自锁 PHASE_SLOW 1 # 慢清偿 PHASE_CLEARABLE 2 # 有效清偿 TRANSITION_MASK 3 # 过渡带掩码仅绘图轮廓不写入phase_label # ODE 定义 def ode_ddt(t, s, beta, g): rhs -s * (1.0 / ((1.0 s)**beta) g) return rhs def mu_debt(s): return s / (s s_c) # solve_single_run完整状态返回全分支status def solve_single_run(beta, g, s0, t_endt_end_global, n_evaln_eval_global): t_eval np.linspace(0, t_end, n_eval) sol solve_ivp(ode_ddt, t_span(0, t_end), y0[s0], args(beta, g), t_evalt_eval, methodLSODA, atol1e-9, rtol1e-7) out { t: None, s: None, mu: None, T_half_S: np.nan, T_half_m: np.nan, status: None } if not sol.success: out[status] integration_failed return out t sol.t s sol.y[0, :] mu mu_debt(s) out[t] t out[s] s out[mu] mu out[status] ok # 检测是否触及s_floor if np.any(s s_floor): out[status] floor_reached # 计算T_half_S half_S_target s0 / 2 idx_S np.argmin(np.abs(s - half_S_target)) if np.abs(s[idx_S] - half_S_target) 0.05 * half_S_target: out[T_half_S] t[idx_S] else: out[status] t_half_not_reached # 计算T_half_m mu0 mu_debt(s0) half_mu_target mu0 / 2 idx_m np.argmin(np.abs(mu - half_mu_target)) if np.abs(mu[idx_m] - half_mu_target) 0.05 * half_mu_target: out[T_half_m] t[idx_m] else: out[status] t_half_not_reached return out # SMOKE TEST g0 锚点校验 def smoke_g0_powerlaw(): print( START SMOKE TEST g0 ANALYTIC ANCHOR ) beta_list [1.0, 1.2, 1.5, 2.0, 3.0] s0 100.0 smoke_records [] for beta in beta_list: run solve_single_run(beta, g0.0, s0s0, t_endt_end_global, n_evaln_eval_global) t run[t] s run[s] mu run[mu] T_half_S run[T_half_S] T_half_m run[T_half_m] delta_T T_half_m - T_half_S mask_fit (s s0 / 10) (s 10) t_fit t[mask_fit] s_fit s[mask_fit] n_fit len(s_fit) s_fit_min np.min(s_fit) if n_fit0 else np.nan s_fit_max np.max(s_fit) if n_fit0 else np.nan trans_idx np.argmax(s s0/10) t_transition t[trans_idx] if np.any(s s0/10) else np.nan slope None theory_slope None fit_ok False fit_type None print(f --- beta {beta} ---) print(fT_half_S {T_half_S:.4f}, T_half_m {T_half_m:.4f}, ΔT {delta_T:.4f}) print(ft_transition {t_transition:.4f}, fit points n{n_fit}) if beta 1.0: fit_type linear print(β1 branch: linear decay regime (not powerlaw)) if n_fit 10: res linregress(t_fit, s_fit) slope res.slope print(fLinear fit slope s~t: {slope:.4f}, expected ~ -1.0) fit_ok True else: fit_type powerlaw theory_slope -1.0 / (beta - 1.0) if n_fit 10: logt np.log10(t_fit) logs np.log10(s_fit) res linregress(logt, logs) slope res.slope print(fLog-log slope: {slope:.4f}, theoretical slope: {theory_slope:.4f}) if np.abs(slope - theory_slope) slope_tol_smoke: fit_ok True print(f✓ Slope passed tolerance ±{slope_tol_smoke}) else: print(f✗ Slope OUTSIDE tolerance!) else: print(⚠️ Not enough points in power-law fitting window) smoke_records.append({ beta: beta, T_half_S: T_half_S, T_half_m: T_half_m, delta_T: delta_T, t_transition: t_transition, slope: slope, theory_slope: theory_slope, fit_ok: fit_ok, fit_type: fit_type, n_fit: n_fit, s_fit_min: s_fit_min, s_fit_max: s_fit_max }) high_beta_recs [r for r in smoke_records if r[beta]1.0] deltas [r[delta_T] for r in high_beta_recs] beta_vals [r[beta] for r in high_beta_recs] is_monotonic_inc np.all(np.diff(deltas) 0) print( GLOBAL CHECK: ΔT separation trend ) print(fβ sequence: {beta_vals}) print(fΔT sequence: {[f{d:.4f} for d in deltas]}) if is_monotonic_inc: print(✓ ΔT T_half^m - T_half^S is monotonically increasing with β) else: print(✗ ΔT not monotonic increasing with β!) smoke_meta { beta_list: [r[beta] for r in smoke_records], s0: 100.0, t_end: t_end_global, slope_fit: [r[slope] for r in smoke_records], slope_theory: [r[theory_slope] for r in smoke_records], slope_pass: [r[fit_ok] for r in smoke_records], T_half_S: [r[T_half_S] for r in smoke_records], T_half_m: [r[T_half_m] for r in smoke_records], delta_T: [r[delta_T] for r in smoke_records], t_transition: [r[t_transition] for r in smoke_records], n_fit: [r[n_fit] for r in smoke_records], s_fit_min: [r[s_fit_min] for r in smoke_records], s_fit_max: [r[s_fit_max] for r in smoke_records], delta_T_monotonic: bool(is_monotonic_inc), timestamp: time.strftime(%Y-%m-%d %H:%M:%S), script_version: smoke_g0_powerlaw v0.1, slope_tol: slope_tol_smoke } return smoke_records, smoke_meta # Smoke 元数据写入独立smoke_meta.npz规避大文件锁 def write_smoke_meta_npz(smoke_meta, filenamesmoke_meta.npz): import json smoke_meta_json json.dumps(smoke_meta, indent2) np.savez(filename, smoke_meta_jsonsmoke_meta_json) print(f✅ Smoke meta saved into {filename}) def load_smoke_meta_from_npz(filenamesmoke_meta.npz): import json with np.load(filename, allow_pickleTrue) as f: json_str f[smoke_meta_json].item() return json.loads(json_str) # fit_slope_from_trace 与smoke同源拟合 def fit_slope_from_trace(t_arr, s_arr, s0): mask_fit (s_arr s0 / 10) (s_arr 10) t_fit t_arr[mask_fit] s_fit s_arr[mask_fit] n_fit len(s_fit) s_fit_min np.min(s_fit) if n_fit0 else np.nan s_fit_max np.max(s_fit) if n_fit0 else np.nan slope np.nan if n_fit 10: logt np.log10(t_fit) logs np.log10(s_fit) res linregress(logt, logs) slope res.slope return { slope: slope, n_fit: n_fit, s_fit_min: s_fit_min, s_fit_max: s_fit_max } # Anchor Consistency Check 返回可序列化JSON字符串 def anchor_consistency_check(smoke_meta, beta_list, g_list, s0_list): import json anchor_report {mismatch_flags: [], records: []} g_min g_list[0] s0_target 100.0 idx_s0 np.argmin(np.abs(s0_list - s0_target)) for ib, beta in enumerate(beta_list): run solve_single_run(beta, g_min, s0_target, t_endt_end_global, n_evaln_eval_global) if run[status] integration_failed: anchor_report[records].append({ beta: beta, g_min: g_min, s0_target: s0_target, status: integration_failed }) continue t_trace run[t] s_trace run[s] T_half_main run[T_half_S] fit_out fit_slope_from_trace(t_trace, s_trace, s0_target) slope_main fit_out[slope] sm_idx smoke_meta[beta_list].index(beta) T_half_smoke smoke_meta[T_half_S][sm_idx] slope_smoke smoke_meta[slope_fit][sm_idx] theory_slope smoke_meta[slope_theory][sm_idx] mismatch_flag False if not np.isnan(slope_main) and not np.isnan(slope_smoke): if np.abs(slope_main - slope_smoke) slope_tol_anchor: mismatch_flag True anchor_report[records].append({ beta: beta, g_min: g_min, s0_target: s0_target, T_half_main: T_half_main, T_half_smoke: T_half_smoke, slope_main: slope_main, slope_smoke: slope_smoke, theory_slope: theory_slope, n_fit_main: fit_out[n_fit], s_fit_min_main: fit_out[s_fit_min], s_fit_max_main: fit_out[s_fit_max], mismatch_flag: mismatch_flag }) if mismatch_flag: anchor_report[mismatch_flags].append(beta) if len(anchor_report[mismatch_flags]) 0: print(f⚠️ SMOKE_ANCHOR_MISMATCH detected at beta {anchor_report[mismatch_flags]}) else: print(✅ Anchor consistency check passed) anchor_report_json json.dumps(anchor_report, indent2) return anchor_report_json # 主扫描 scan_phase_diagram修复g_c收集bug def scan_phase_diagram(beta_list, g_list, s0_list): n_beta len(beta_list) n_g len(g_list) n_s0 len(s0_list) phase_label np.zeros((n_beta, n_g, n_s0), dtypeint) T_half_S_grid np.full((n_beta, n_g, n_s0), np.nan) T_half_m_grid np.full((n_beta, n_g, n_s0), np.nan) D_dec_grid np.full((n_beta, n_g, n_s0), np.nan) transition_mask np.zeros((n_beta, n_g, n_s0), dtypebool) # 阈值定义沙盘分类判据 t_clear_thresh 50.0 t_lock_thresh 150.0 for ib, beta in enumerate(beta_list): for ig, g in enumerate(g_list): r_end_list [] for is0, s0 in enumerate(s0_list): run solve_single_run(beta, g, s0) T_half_S run[T_half_S] T_half_m run[T_half_m] T_half_S_grid[ib,ig,is0] T_half_S T_half_m_grid[ib,ig,is0] T_half_m D_dec_grid[ib,ig,is0] T_half_m - T_half_S # 分类 if np.isnan(T_half_S) or T_half_S t_lock_thresh: pl PHASE_LOCKED elif T_half_S t_clear_thresh: pl PHASE_CLEARABLE else: pl PHASE_SLOW phase_label[ib,ig,is0] pl # 过渡带掩码 if 40 T_half_S 60: transition_mask[ib,ig,is0] True # 收集终态r用于g_c临界梯度 s_end run[s][-1] if run[t] is not None else np.nan r_end_list.append(s_end) # 循环内收集完成再求临界g_c r_series np.array(r_end_list) # 临界梯度判定示例r_end从高值跌落点 pass return phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask # 绘图plot_phase_slices不阻塞仅保存图片 def plot_phase_slices(beta_list, g_list, s0_list, phase_label, transition_mask): plt.rcParams[font.size] 9 for ib, beta in enumerate(beta_list): fig, ax plt.subplots(figsize(7,5)) im ax.pcolormesh(s0_list, g_list, phase_label[ib,:,:], vmin0, vmax2, shadingauto) cs ax.contour(s0_list, g_list, transition_mask[ib,:,:], colorswhite, linestyles--) ax.set_xlabel(r$s_0$) ax.set_ylabel(r$g$) ax.set_title(fβ {beta:.2f}) cbar plt.colorbar(im, axax) cbar.set_ticks([0,1,2]) cbar.set_ticklabels([Locked,Slow,Clearable]) plt.tight_layout() plt.savefig(fphase_slice_beta_{beta:.2f}.png, dpi150) plt.close(fig) print(✅ Phase slices saved to image files) # 合并归档写入phase_scan_result.npz def merge_all_to_phase_npz(phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask, beta_list, g_list, s0_list, anchor_json, smoke_meta_pathsmoke_meta.npz): with np.load(smoke_meta_path, allow_pickleTrue) as f: smoke_meta_json f[smoke_meta_json] np.savez(phase_scan_result.npz, phase_labelphase_label, T_half_S_gridT_half_S_grid, T_half_m_gridT_half_m_grid, D_dec_gridD_dec_grid, transition_masktransition_mask, beta_listnp.array(beta_list), g_listnp.array(g_list), s0_listnp.array(s0_list), smoke_meta_jsonsmoke_meta_json, anchor_check_jsonanchor_json) print(✅ Merged all data into phase_scan_result.npz) # MAIN ENTRY POINT if __name__ __main__: # 1. Run smoke smoke_records, smoke_meta smoke_g0_powerlaw() # Hard gate: smoke fail - abort main scan powerlaw_records [r for r in smoke_records if r[beta] 1.0] if not all(r[fit_ok] for r in powerlaw_records): raise RuntimeError(SMOKE FAILED: power-law anchor mismatch, abort main scan) # Write smoke meta standalone npz write_smoke_meta_npz(smoke_meta) # 2. Define scan grid beta_list [1.0, 1.2, 1.5, 2.0, 3.0] g_list np.logspace(-4, 1, 40) s0_list np.logspace(0, 3, 30) # 3. Main scan print( START MAIN PHASE DIAGRAM SCAN ) phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask scan_phase_diagram(beta_list, g_list, s0_list) # 4. Anchor consistency check anchor_json anchor_consistency_check(smoke_meta, beta_list, g_list, s0_list) # 5. Merge and save big npz merge_all_to_phase_npz(phase_label, T_half_S_grid, T_half_m_grid, D_dec_grid, transition_mask, beta_list, g_list, s0_list, anchor_json) # 6. Plot slices (non-blocking, save only) plot_phase_slices(beta_list, g_list, s0_list, phase_label, transition_mask) print( PIPELINE COMPLETED )参考来源激光淬火技术comsol相变模拟的实践与应用陶瓷淬火时“啪“一声裂开的瞬间背后藏着相场模型里的连续损伤演化。今天咱们用Matlab玩个热应力场相场断裂的耦合计算看看脆性材料怎么被温度场玩坏磨削区仿真技术深度解析单颗磨粒作用与计算建模实战重型自卸汽车离合器系统多物理场耦合设计与工程验证——基于热-机-流协同调控的集成解决方案电磁场理论在多物理场耦合仿真中的工程化解读与实战指南
返回列表