船舶作为海上运输、作业等的重要工具，其安全性受到广泛关注．作为影响船舶安全性主要因素之一的船体结构疲劳破坏问题，研究者在船体结构的设计中给予了特别的关注．而为兼顾经济性和安全性的要求，船体结构的优化设计是必不可少的．文献[1]提出一种新的基于机器学习的船舶设计优化方法，并对某散货船结构进行了优化设计；文献[2]提出一种结合遗传算法和优化方法的混合遗传算法(GA)，实现了减轻结构质量的目的；文献[3]使用改进的果蝇算法(IFOA)和IFOA-BP(反向传播)方法对影响浮式生产储油卸油装置(FPSO)船体梁疲劳强度的主要设计参数进行优化；文献[4]基于参数化子模型法及疲劳强度评估的设计波法，提出了一套结构抗疲劳设计的优化方法；文献[5]在某大型SWATH船波浪载荷和结构强度直接计算的基础上，提出了一种子模型法结合参数化模型的优化计算方法；文献[6]为研究三体船结构细节的疲劳行为，提出了简化疲劳分析方法．事实上，采用船舶相关规范要求、基于经验公式得到的计算载荷[1，7]，虽然计算简便，但其与船体结构实际遭受的载荷会有一定的出入；采用设计波载荷[4，8]，又涉及到设计波的选取情况会影响船体结构优化结果的问题．而船体结构疲劳评估直接计算的谱分析法[6，9-11]，具有较高的精度、能够反映船舶具体结构细节等不可替代的优点，但其计算过程较复杂，计算工作量十分繁重，目前仅限于应用在能够确定船体结构的准确优化方案的情况；而若想将其应用于须要寻求船体结构优化方案的情况，则还有待进一步研究．本研究以船体结构应力响应函数曲线为研究对象，在考虑计算结果误差控制的前提下，提出了船体结构计算频率的筛选方法；分析了船体结构应力响应函数曲线形状的快速识别和确定的可行性，并总结给出了确定的方法；同时结合船体结构浪向疲劳累积损伤度的计算方法，提出了船体结构疲劳评估直接计算的简化方法．1 浪向疲劳累积损伤度的计算方法考虑Nload个装载状态，船体结构的疲劳累积总损伤度D可表示为[12-13]D=δTLa¯Γ1+m2∑n=1Nloadpn∑i=1nS∑j=1nHpipjvijn(22m0ijn)m,式中：δ为在航率系数；TL为船舶疲劳计算的回复期；a¯，m为所用疲劳曲线(S-N曲线)的两个参数；Γ(1+m/2)为伽玛函数；pn为第n个装载的时间分配系数；nS为选定的海况散布图所包含的短期海况的总数；nH为计算浪向的总数；pi为第i个海况出现的概率；pj为第j个浪向出现的概率；υijn为第n个装载状态、海况i和浪向j的应力响应过零率，νijn=[1/(2π)]m2ijn/m0ijn；m0ijn和m2ijn为第n个装载状态及海况i和浪向j下应力响应谱的零阶矩和二阶矩，且m0,2ijn=∫0+∞ω0,2GXX(ω)dω．从而在某一装载状态下，船体结构的浪向疲劳累积损伤度可表示为Dj=pjδTLa¯Γ1+m2pn∑i=1nSpivijn(22m0ijn)m.特别地，结合m0ijn和m2ijn具体的计算方法，船体结构的浪向疲劳累积损伤度可进一步表示为Dj=[pj(δTL/a¯)]Γ(1+m/2)pn[1/(2π)]∙(22)m∑i=1nSpi∑k=1lωek2σA2(ωk)Gηη(ωk)Δωk1/2∙∑k=1lσA2(ωk)Gηη(ωk)Δωk(m-1)/2,式中：ωek为船舶的遭遇频率，ωek=ωk(1+(ωkU/g)cos θ)；θ为航向角；U为船舶的航速；σA(ωk)为与计算频率ωk对应的船体结构应力响应；Gηη(ωk)为用有义波高Hs和平均跨零周期Tz来表达的波浪谱，且Gηη(ωk)=Hs24π2πTz4ωk-5exp-1π2πTz4ωk-4;Δωk=ωk+1-ωk    (k=1),(ωk+1-ωk-1)/2    (1kl),ωk-ωk-1    (k=l),其中l为计算频率的总个数．令Rj=pjδTLa¯Γ1+m2pn12π(22)m，则船体结构的浪向疲劳累积损伤度可进一步表示为Dj=Rj∑i=1nSpi∑k=1lωek2σA2(ωk)Gηη(ωk)Δωk1/2∙∑k=1lσA2(ωk)Gηη(ωk)Δωk(m-1)/2. (1)2 疲劳评估直接计算的简化方法考虑到在采用疲劳谱分析方法开展船体结构的优化设计研究中，须要重复进行“获得优化方案→获得优化模型→进行结构分析”这样耗费时间和精力的过程，因此本研究对该过程进行合理简化．2.1　船体结构计算频率的筛选方法A. 由波浪谱曲线确定的研究频率根据选定的海况散布图，取有义波高Hs为最小值，获得由Hs与各平均跨零周期Tz组成的全部参数组合(Hs，Tz)．结合波浪谱的计算公式，得到与各参数组合(Hs，Tz)对应的短期海况的波浪谱曲线；统计各波浪谱曲线峰值对应的计算频率，获得由波浪谱曲线确定的研究频率．B. 由基准应力响应函数曲线确定的研究频率采用下式处理由船舶原设计方案得到的应力响应函数曲线，获得与之对应的船体结构基准应力响应函数曲线，记为基准Qk0(ω)曲线，Qk(ωk)=σA(ωk)/σA(ω0), (2)式中：σA(ω0)和ω0分别为应力响应函数曲线最高峰值对应的结构应力响应和计算频率，且将ω0记为最危险频率；σA(ωk)为船体结构应力响应函数曲线中与各计算频率ωk对应的结构应力响应．然后，研究与各计算频率对应的Qk0(ωk)值，对船体结构的基准Qk0(ω)曲线进行区域划分，令Qk0(ωk)≥0.5的部分为高应力区，将剩余部分中的低频部分记为低频低应力区，高频部分记为高频低应力区，且规定在应力响应函数曲线最危险频率左侧的计算频率为低频频率，在其右侧的计算频率为高频频率，进而可获得对应的研究频率．C. 由Qk02(ω)Gηη(ω)曲线确定的研究频率结合前述内容，可获得与各短期海况对应的Qk02(ω)Gηη(ω)曲线；统计各曲线峰值对应的计算频率，获得由Qk02(ω)Gηη(ω)曲线确定的研究频率．D. 船体结构计算频率的筛选方法对于由前述内容确定的研究频率，若其满足如下条件，则应将其筛选出来，并称其为船体结构的主要研究频率．a. 由Qk02(ω)Gηη(ω)曲线确定的研究频率；b. 既是由船体结构基准Qk0(ω)曲线确定的高应力区的研究频率，又是由波浪谱曲线确定的研究频率；c. 经过条件b筛选后剩余的船体结构基准Qk0(ω)曲线高应力区的研究频率，且在由波浪谱曲线确定的研究频率的最值所构成范围内；d. 船体结构基准Qk0(ω)曲线低应力区的研究频率，同时包含在由波浪谱曲线确定的研究频率内，且与其对应的Qk0(ωk)值满足Qk0(ωk)β，规定β的取值范围为0.1≤β≤0.4．2.2　应力响应函数曲线形状的识别和确定方法确定船体结构的主要研究频率后，可对应地获得船体结构优化方案在这些主要研究频率的结构应力响应σA'(ωk)；然后利用式(2)处理上述σA'(ωk)值，得到与各σA'(ωk)值对应的Qk(ωk)值；用上述Qk(ωk)值对应地更新Qk0(ωk)值，便获得了船体结构优化方案在对应计算浪向的Qk(ω)曲线．由公式ΔQk(ωk)=Qk(ωk)-Qk0(ωk) (3)可分析各Qk(ωk)值和其基准Qk0(ωk)值之间的变化情况，若判定各Qk(ωk)值均不变，则说明船体结构优化方案在对应计算浪向的应力响应函数曲线的形状未发生变化．根据下式可计算得到船体结构优化方案在对应计算浪向各计算频率的结构应力响应，进而可获得与之对应的应力响应函数曲线，σA'(ωk)=Qk(ωk)σA'(ω0), (4)式中：σA'(ωk)为船体结构优化方案在计算频率ωk的结构应力响应；σA'(ω0)为在最危险频率ω0下，根据船体结构优化方案计算得到的结构应力响应．3 简化方法的应用对象本研究选择某个极具代表性、结构疲劳损伤很大的疲劳热点结构作为本文提出的简化方法的应用对象．该应用对象主要由主甲板、上层建筑后端壁板、主甲板下的连续纵桁、主甲板下的横舱壁板及在后端壁板处中断的纵壁板组成，见图1．10.13245/j.hust.211208.F001图1应用对象结构示意图在该应用对象的结构分析中，选取12个计算浪向(0°～330°，间隔为30°)，18个计算频率(0.1～1.8 Hz，间隔为0.1 Hz)，计算海况取为全球海况．首先，由确定的计算频率，结合选择海况的海况散布图，根据由波浪谱曲线确定研究频率的方法，获得对应的研究频率{0.3，0.4，0.5，0.6，0.7，0.8，1.0，1.3}．以0°计算浪向(j=1，Dj0=D10=7.688)为例，应用对象的应力响应函数曲线最高峰值对应的结构应力响应值σA(ω0)=76.09 MPa，最危险频率ω0=0.6 Hz，与各计算频率对应的σA(ωk)值见表1．10.13245/j.hust.211208.T001表1应用对象的结构应力响应汇总ωk/HzσA(ωk)/MPaωk/HzσA(ωk)/MPa0.10.501.014.130.22.381.117.150.310.781.215.320.430.541.315.200.557.791.410.730.676.091.510.690.773.861.69.550.864.051.76.400.946.101.84.28采用式(2)处理表1中的各σA(ωk)值，得到应用对象的基准Qk0(ω)曲线．然后根据由基准Qk0(ω)曲线确定研究频率的方法，得到高应力区的研究频率{0.5，0.6，0.7，0.8，0.9}、低频低应力区的研究频率{0.1，0.2，0.3，0.4}及高频低应力区的研究频率{1.0，1.1，1.2，1.3，1.4，1.5，1.6，1.7，1.8}．之后根据由Qk02(ω)Gηη(ω)曲线确定研究频率的方法，获得对应的研究频率{0.4，0.5，0.6，0.7，0.8，0.9，1.3}．最后，结合船体结构计算频率的筛选方法，分别取β=0.1，0.2，相应地确定应用对象的主要研究频率Yfβ依次为：Yf 0.1={0.3，0.4，0.5，0.6，0.7，0.8，0.9，1.0，1.3}；Yf 0.2={0.4，0.5，0.6，0.7，0.8，0.9，1.3}．4 简化方法的验证与单一计算频率对应的Qk(ωk)值改变对Dj的影响见图2．特别地，规定ΔDj=Dj0–Dj (下同)．10.13245/j.hust.211208.F002图2单一计算频率对Dj的影响由图2可知：对于高应力区的研究频率ωk=0.5，0.6，0.7，0.8，0.9和紧邻高应力区的计算频率ωk=0.4 (Qk(ωk)=0.4)，与之对应的Qk(ωk)值改变对Dj有很大的影响．同样地，与多个计算频率对应的Qk(ωk)值改变对Dj的影响见图3．由图3可知：对于计算频率ωk=0.1，0.2，与之对应的Qk(ωk)值改变对Dj几乎无影响；对于计算频率ωk=1.4，1.5，1.6，1.7，1.8，与之对应的Qk(ωk)值的改变对Dj的影响不大；而对于高应力区的计算频率，与之对应的Qk(ωk)值改变会对Dj造成相当大的影响．由上述分析可知：本研究提出的船体结构计算频率的筛选方法，不仅可以筛选出对Dj贡献大的计算频率，而且能够将对Dj贡献很小的计算频率筛选掉． 10.13245/j.hust.211208.F003 图3 多个计算频率对Dj的影响 1—ωk=0.3，0.4；2—ωk=0.2，0.3，0.4；3—ωk=0.1，0.2，0.3，0.4；4—ωk=0.4；5—ωk=0.7，0.8；6—ωk=0.5，0.7，0.8； 7—ωk=0.5，0.7，0.8，0.9；8—ωk=0.7；9—ωk=1.0；10—ωk=1.0，1.1；11—ωk=1.0，1.1，1.2；12—ωk=1.0，1.1，1.2，1.3；13—ωk=1.0，1.1，1.2，1.3，1.4；14—ωk=1.0，1.1，1.2，1.3，1.4，1.5；15—ωk=1.0，1.1，1.2，1.3，1.4，1.5，1.6；16—ωk=1.0，1.1，1.2，1.3，1.4，1.5，1.6，1.7；17—ωk=1.0，1.1，1.2，1.3，1.4，1.5，1.6，1.7，1.8． 进一步，考虑是否存在合适的ΔQ0值，当Qk(ωk)值的变化在±ΔQ0值确定的范围内时，就可以认为Qk(ωk)值无变化．本研究特别关注ΔQ0=±0.01和ΔQ0=±0.05的情况．根据已确定的主要研究频率Yf 0.1和Yf 0.2，本研究分析了与Yf 0.1和Yf 0.2对应的Qk(ωk)值改变对Dj的影响，具体见图4．由图4可以看出：当|ΔQk(ωk)|=0.01时，有|ΔDj/Dj0|=2.6%；而当|ΔQk(ωk)|=0.05时，有|ΔDj/Dj0|12%．10.13245/j.hust.211208.F004图4主要研究频率对Dj的影响为不失一般性且能说明问题，本研究规定与主要研究频率Yfβ对应的ΔQk(ωk)值取值见表2，并认为这是最危险的状态，继续分析Qk(ωk)值改变对Dj的影响．10.13245/j.hust.211208.T002表2与Yfβ对应的ΔQk(ωk)取值情况ωk/HzΔQk(ωk)Yfβ-1Yfβ-20.3-0.050.050.4-0.050.050.5-0.010.010.6-0.010.010.7-0.010.010.8-0.010.010.9-0.010.011.0-0.050.051.3-0.050.05对应地，与各主要研究频率对应的Qk(ωk)值改变对Dj的影响见表3．由表3知：对于与Yfβ (β=0.1,0.2)对应的ΔQk(ωk)取不同值时，Dj的误差基本在4%左右．10.13245/j.hust.211208.T003表3与各主要研究频率对应的Qk(ωk)值改变对Dj的影响参数Yf 0.1Yf 0.2Yf 0.1-1Yf 0.1-2Yf 0.2-1Yf 0.2-2Dj7.3918.0157.4167.980ΔDj0.297-0.3270.271-0.292ΔDj/Dj00.039-0.0430.035-0.038结合图4和表3的结果，本研究确定了合适的ΔQ0，具体如下：a. 对于高应力区中的主要研究频率，取ΔQ0=0.01；b. 对于低应力区中的主要研究频率，若对计算结果的精度要求高，取ΔQ0=0.01；若对计算结果的精度要求不是特别高，取ΔQ0=0.05．5 简化方法的应用实例针对前文确定的应用对象，采用本研究提出的简化方法对其进行结构优化设计，根据获得的结构优化方案，得到其在0°计算浪向( j=1)各主要研究频率的结构应力响应σA'(ωk)，具体见表4．10.13245/j.hust.211208.T004表4与主要研究频率Yf 0.1对应的结构应力响应( j=1)ωk/Hz0.30.40.50.60.70.80.91.01.3σA'(ωk)/MPa10.2523.9945.4359.8658.1450.4036.2218.8512.51采用式(2)处理表5中的各结构应力σA'(ωk)，得10.13245/j.hust.211208.T005表5与主要研究频率Yf 0.1对应的ΔQk(ωk)值ωk/Hz0.30.40.50.60.70.80.91.01.3|ΔQk(ωk)|0.029 50.000 60.000 60.000 00.000 60.000 20.000 80.129 20.009 2到与其对应的Qk(ωk)值；用各Qk(ωk)值对应地更新基准Qk0(ωk)值，获得该结构优化方案的Qk(ω)曲线，见图5．10.13245/j.hust.211208.F005图5该结构优化方案的Qk(ω)曲线由式(3)判断该结构优化方案应力响应函数曲线的形状是否改变，并规定ΔQ0=0.01，与主要研究频率Yf 0.1对应的ΔQk(ωk)值见表5．根据式(4)，由该结构优化方案的Qk(ω)曲线，可对应地得到其应力响应函数曲线，具体见图6．10.13245/j.hust.211208.F006图6该结构优化方案的应力响应函数曲线结合已获得的该结构优化方案的应力响应函数曲线，由本研究提出的简化方法可计算得到0°计算浪向( j=1)的疲劳累积损伤度D1'(Yf 0.1)=3.786 8，D1'(Yf 0.2)=3.743 4；同时，由全频率的谱分析方法可计算得到0°计算浪向( j=1)的疲劳累积损伤度D1=3.786 4；比较两种方法的计算结果，且规定ΔD1=D1-D1'，从而有ΔD1(Yf 0.1)/D1=-0.000 1，ΔD1(Yf 0.2)/D1=0.011 4．6 结语a. 采用疲劳谱分析法对船体复杂结构进行优化研究时，须要重复进行“获得优化方案→获得优化模型→进行结构分析”的过程，十分耗费时间和精力．本研究提出了一种船体结构疲劳评估直接计算的简化方法，可在很大程度上对结构的优化计算过程进行简化，且能达到很好的误差控制．b. 不同的计算频率对船体结构浪向疲劳累积损伤度的贡献不同，本研究提出了一种船体结构计算频率的筛选方法，通过将贡献大的计算频率确定为主要研究频率，实现了对计算频率的简化．c. 应力响应函数曲线是船体结构疲劳评估中十分重要的内容，结合船体结构的主要研究频率，本研究提出了应力响应函数曲线形状的快速识别和确定方法，实现了对船体结构优化方案的应力响应函数曲线的快速确定，节省了研究时间和人力．本研究提出的船体结构疲劳评估直接计算的简化方法提高了船体结构优化研究效率，减少了人力，同时也保证了研究结果的准确性，可具体应用于复杂工程结构的结构优化研究．