沥青混合料细观结构的三维粘弹本构及虚拟力学试验方案【附数据】
✨ 长期致力于沥青混合料、X-ray CT、细观结构、粘弹本构、三轴蠕变、虚拟力学研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》
（1）基于CT切片图像的沥青混合料三维细观结构重建与材质分割改进算法：
针对沥青混合料X-ray CT图像中集料、沥青胶浆和空隙灰度分布重叠的问题，提出一种改进的多阈值材质分割算法。首先对CT切片（分辨率0.1mm/pixel，尺寸1024x1024）进行非局部均值滤波去噪，然后采用多尺度Retinex增强对比度。材质分割采用改进的Otsu双阈值法：寻找两个最优阈值T1和T2，使得类间方差加权和最大，将图像分为空隙（灰度< T1）、胶浆（T1<=灰度< T2）和集料（灰度>=T2）。针对集料颗粒间的粘连现象，提出形态学多尺度分割算法：依次使用半径为2、4、6像素的圆形结构元素进行腐蚀，分离粘连颗粒后再膨胀恢复。在30张连续切片上重建三维模型，采用Marching Cubes算法生成表面网格，得到体素尺寸1mm x 1mm x 1mm的数值模型。利用天河二号超级计算机，对体素边长0.5mm、1mm、2mm的三种模型进行单轴蠕变虚拟试验，分析尺寸效应。结果表明，当体素边长不大于1mm时，蠕变劲度模量的变化小于3%，因此推荐1mm作为标准体素尺寸。重建模型包含约200万个单元，准确反映了沥青混合料的三维细观结构。
（2）沥青砂胶三维粘弹本构模型的UMAT子程序开发与参数确定：
针对沥青砂胶（沥青+矿粉）的粘弹性行为，采用广义Maxwell模型（8个麦克斯韦单元串联）描述三维本构关系。推导了积分型本构方程的增量形式，利用FORTRAN编写ABAQUS UMAT子程序。该子程序在每个增量步中根据应变增量计算应力更新，采用递归公式计算 hereditary integral。粘弹参数（松弛模量Prony级数系数）的确定采用间接方法：进行单轴蠕变试验（40℃，应力0.1MPa，时长3600秒），获得蠕变柔量曲线，然后通过非线性最小二乘拟合得到Prony级数的8个松弛时间和对应的模量系数。对于泊松比，假设为时间依赖的粘弹性泊松比，其瞬时值0.35，长期值0.49，通过Laplace变换求解。对UMAT子程序进行验证：模拟简单压缩试验（应变率0.001/s），数值结果与实验室实测蠕变曲线对比，最大相对误差4.6%。弯曲试验（三点弯曲，跨径80mm，加载速率5mm/min）的模拟力-位移曲线与实测结果吻合，峰值力误差2.1%。该子程序可准确预测沥青砂胶在不同荷载条件下的力学行为。
（3）基于细观模型的沥青混合料三轴蠕变虚拟试验及参数适用范围分析：
将重建的三维细观模型（含集料、胶浆和空隙）与开发的UMAT子程序结合，在ABAQUS中进行三轴蠕变虚拟试验。试验条件：围压100kPa、200kPa，偏应力50kPa、100kPa、150kPa，温度40℃，蠕变时长7200秒。模型底部固定，顶部施加轴向荷载，围压通过压力边界实现。计算节点数约50万，使用32核并行计算，单次虚拟试验耗时约2小时。虚拟试验得到的蠕变柔量曲线与真实三轴蠕变试验（采用相同的沥青混合料AC-13）对比，在围压200kPa、偏应力100kPa条件下，1000秒时的蠕变柔量虚拟值为3.2e-6 Pa-1，实测值为3.0e-6 Pa-1，相对误差6.7%。对不同偏应力水平的模拟表明，当偏应力大于150kPa时，虚拟试验出现加速蠕变（损伤累积），而真实试件在150kPa时仍处于稳态蠕变，说明模型对于高应力水平预测偏保守，适用范围建议为偏应力不大于120kPa。空隙率对蠕变行为的影响：空隙率从4%增加到8%，蠕变应变增加32%，与文献规律一致。该方法为沥青混合料细观力学研究提供了可靠的虚拟试验平台。
import numpy as np import cv2 from skimage.measure import marching_cubes def ct_segmentation_and_reconstruction(ct_stack, threshold1=70, threshold2=120): # 改进的多阈值分割与粘连分割 segmented = np.zeros_like(ct_stack, dtype=np.uint8) for i, slice_img in enumerate(ct_stack): # 非局部均值滤波 denoised = cv2.fastNlMeansDenoising(slice_img, h=10, templateWindowSize=7, searchWindowSize=21) # 多尺度Retinex增强（简化） enhanced = cv2.addWeighted(denoised, 1.5, cv2.GaussianBlur(denoised, (0,0), 5), -0.5, 0) # 双阈值分割 _, thresh1 = cv2.threshold(enhanced, threshold1, 255, cv2.THRESH_BINARY) _, thresh2 = cv2.threshold(enhanced, threshold2, 255, cv2.THRESH_BINARY) # 空隙、胶浆、集料编码 slice_seg = np.zeros_like(enhanced, dtype=np.uint8) slice_seg[enhanced < threshold1] = 0 # 空隙 slice_seg[(enhanced >= threshold1) & (enhanced < threshold2)] = 1 # 胶浆 slice_seg[enhanced >= threshold2] = 2 # 集料 # 形态学多尺度分离集料粘连 mask_agg = (slice_seg == 2).astype(np.uint8) * 255 kernel_small = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)) kernel_med = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5)) eroded = cv2.erode(mask_agg, kernel_small, iterations=1) dilated = cv2.dilate(eroded, kernel_med, iterations=1) slice_seg[dilated > 0] = 2 segmented[i] = slice_seg # 三维重建（Marching Cubes） verts, faces, normals, values = marching_cubes(segmented == 2, level=0.5) return segmented, verts, faces def prony_fitting(creep_compliance, time, n_terms=8): # 广义Maxwell模型Prony级数拟合 from scipy.optimize import curve_fit def prony_model(t, E_inf, *tau_eta): # tau_eta: [tau1, eta1, tau2, eta2, ...] result = 1.0 / E_inf for i in range(n_terms): tau = tau_eta[2*i] eta = tau_eta[2*i+1] result += (1.0/eta) * (1 - np.exp(-t/tau)) return result # 初始猜测 p0 = [1e8] + [0.1, 1e6] * n_terms popt, _ = curve_fit(prony_model, time, creep_compliance, p0, maxfev=5000) E_inf = popt[0] tau_list = popt[1::2] eta_list = popt[2::2] return E_inf, tau_list, eta_list # 示例：生成模拟CT数据 ct_sim = np.random.randint(0, 200, (20, 200, 200), dtype=np.uint8) seg, verts, faces = ct_segmentation_and_reconstruction(ct_sim) print(f'重建顶点数: {len(verts)}, 三角面片数: {len(faces)}')