Python遍历双列表实现层合板剪切应力计算的问题排查

我要实现层合板剪应力的计算公式，其中zc是从z_[0]到z_[-1]的等距序列（示例：t_array = np.arange(z_[0],z_[-1],0.01)，当前z_[0] = -15.336，z_[-1] = 15.336）。需要遍历t_array，当zc到达z_的当前分界值时，切换到下一个层的参数计算。

现有代码生成的结果下半部分偏移，预期图表x轴首尾都是0，但实际不符合；同时需要代码适配任意层数，满足len(b) = len(nk)且len(z_) = len(b)+1。

现有代码片段：

import numpy as np

z_ = np.array([-15.336,-14.5,-14,14,14.5,15.336])  #层合板中面到各面的距离(zk)
nk = np.array([1,0.58,0.004126,0.58,1])          #归一化系数
b = np.array([300,300,300,300,300])             #层合板宽度(mm)

t_array = np.arange(z_[0],z_[-1],0.001)
t_array = np.append(t_array, z_[-1])
Itot = 146790.02
Qhat = np.zeros(len(t_array))

r = 0
# 遍历t_array
for i, t in enumerate(t_array):
    # 根据公式计算Qhat值
    if r == 1 :
        Qhat[i] = b[r] * nk[r] / 2 * (t ** 2 - z_[r+1] ** 2)
    elif r == len(z_) - 1:
        Qhat[i] = b[r] * nk[r] / 2 * (t ** 2 - z_[r-1] ** 2)
    else:
        Qhat[i] = b[r] * nk[r] / 2 * (t ** 2 - z_[r-1] ** 2)
    
    # 当t_array到达当前z_值时，切换到下一个z_索引
    if r < len(z_) - 1 and t >= z_[r + 1]:
        r += 1

Qhat = Qhat+Qhat_const
tau = -(9.81/Itot)*Qhat
plt.plot(tau,t_array)

问题根源 #

1. 初始层计算错误：r=0时进入else分支，误用z_[r-1]（即z_[-1]），导致初始段计算逻辑完全错误，引发下半部分偏移。
2. 分界判断精度问题：t >= z_[r+1]的浮点数判断容易出现提前/延迟切换层的情况。
3. 公式分支冗余混乱：r==1的分支无公式依据，属于错误逻辑。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# 支持任意层数：满足len(z_) = len(b)+1，len(b)=len(nk)
z_ = np.array([-15.336, -14.5, -14, 14, 14.5, 15.336])  # 层合板中面到各面的距离
nk = np.array([1, 0.58, 0.004126, 0.58, 1])          # 归一化系数
b = np.array([300, 300, 300, 300, 300])             # 层合板宽度(mm)
Itot = 146790.02
Qhat_const = 0  # 根据实际需求赋值

# 生成等距序列，确保包含首尾端点，避免浮点数精度问题
t_array = np.linspace(z_[0], z_[-1], num=int((z_[-1]-z_[0])/0.001)+1)
Qhat = np.zeros(len(t_array))

# 按层计算，每层对应区间[z_[r], z_[r+1]]，r为b/nk的索引
for r in range(len(b)):
    # 筛选当前层的t值索引
    layer_mask = (t_array >= z_[r]) & (t_array <= z_[r+1])
    t_layer = t_array[layer_mask]
    
    if r == 0:
        # 第一层：保证z_[0]处Qhat为0
        Qhat[layer_mask] = b[r] * nk[r] / 2 * (t_layer ** 2 - z_[r] ** 2)
    elif r == len(b)-1:
        # 最后一层：保证z_[-1]处Qhat为0
        Qhat[layer_mask] = b[r] * nk[r] / 2 * (z_[r+1] ** 2 - t_layer ** 2)
    else:
        # 中间层：累加前序层的常量值，再计算当前层变量部分
        prev_total = 0
        for prev_r in range(r):
            prev_total += b[prev_r] * nk[prev_r] / 2 * (z_[prev_r+1] ** 2 - z_[prev_r] ** 2)
        Qhat[layer_mask] = prev_total + b[r] * nk[r] / 2 * (t_layer ** 2 - z_[r] ** 2)

# 计算剪应力并绘图
Qhat += Qhat_const
tau = -(9.81 / Itot) * Qhat

plt.plot(tau, t_array)
plt.xlabel('剪应力τ')
plt.ylabel('zc')
plt.grid(True)
plt.show()

关键修正说明 #

• 分层逻辑重构：用mask直接筛选每层对应的t_array区间，避免循环中索引切换的精度问题，天然支持任意层数。
• 公式边界修正：第一层和最后一层的公式分别保证首尾端点Qhat为0，匹配预期图表要求。
• 序列生成优化：用np.linspace替代np.arange+append，彻底解决浮点数精度导致的端点缺失问题。
• 中间层累加逻辑：正确累加前序层的常量值，确保层间计算连续。

内容的提问来源于stack exchange，提问作者Justincomp