一种基于MATLAB进行P‑Y曲线的自动化计算方法与流程

本发明涉及一种计算方法，尤其涉及一种基于MATLAB进行P-Y曲线的自动化计算方法。

背景技术：

海上风机基础结构机理较为复杂，影响因素较多，虽然研究人员关注这一课题已有三十来年，但迄今为止仍未形成完善的分析理论和方法，很多国家的规范中也还没有相关的内容和规定。目前，我国正在编撰海上风机基础结构方面的行业标准。

在进行海上风电机组基础结构设计分析时，必须要考虑土对桩的水平向约束作用。目前计算水平载荷作用下桩土相互作用常用的还是地基反力法，即忽略土的连续性，认为某点的抗力只与该点的位移有关。地基反力法应用Winkler地基模型，把桩周土离散为一个个单独作用的弹簧。在地基反力法中应用最广泛的是P-Y曲线法。根据《港口工程桩基规范》，P-Y曲线法是由复杂的多个公式和图表查询共同手算得到，完成一个海上风机基础结构的P-Y曲线数据对的计算工作量相当大。而且，目前还没有任何一个软件或者相似方法可以实现直接提供P-Y曲线数据对。

目前阶段，海上风机基础结构形式主要包括单桩基础、高桩承台基础、多脚架基础、导管架基础、重力式基础以及其他概念型基础等。常用的基础型式中导管架基础在海上石油平台施工中有一定经验，也涉及到桩基础的建设；重力式基础体积大、重量大、安装、运输不方便，施工周期长，对海床条件要求较高；而桩基基础具有结构构造简单，受力明确、易于建造和海上施工方便而被广泛应用于海上风机基础结构中。

桩土水平相互作用属于海上风机桩基基础设计中最为关键性的技术，适用于所有海上风机桩基基础的建模计算。但是，目前P-Y曲线参数计算和其他作用一样仅靠手算进行，而且由于插值图表的存在，其手算计算难度相当大，工作量大，成本高，错误率高。

技术实现要素：

为了解决上述问题中的不足之处，本发明提供了一种基于MATLAB进行P-Y曲线的自动化计算方法。

为了解决以上技术问题，本发明采用的技术方案是：一种基于MATLAB进行P-Y曲线的自动化计算方法，它包括以下步骤：

(一)：首先，根据地质勘察单位提供的地勘资料，填写六个参数并输入到MATLAB中调用：

(1a)软黏土包括直径、土层号、三轴仪试验中最大主应力差一半时的应变值、土的重度、原状黏土不排水抗剪强度的标准值、泥面以下桩的任一深度；

依据地质勘查报告获得土的物理性质参数：不排水抗剪强度标准值C_u、土的重度γ、泥面以下桩的任一深度Z、系数ζ、桩径D、相关系数ρ、泥面以下Z深度处桩的侧向水平变形Y、三轴仪试验中最大主应力差一半时的应变值ε₅₀；

(1b)砂土包括直径、土层号、土层位于水上/水下、土的重度、内摩擦角、泥面以下桩的任一土层深度；

依据地质勘查报告获得土的物理性质参数：内摩擦角φ、泥面以下桩的任一深度Z、桩径D、土的重度γ、计算系数Ψ、泥面以下Z深度处桩的侧向水平变形Y；

(二)：根据输入的参数对黏土或砂土进行数据处理：

(2a)软黏土的数据处理计算过程如下：

不排水抗剪强度标准值C_u小于等于96kPa的软黏土，在非反复荷载作用下P－Y曲线按下列规定确定：

(2a.1)极限水平土抗力转折点的深度Z_r，可按下式计算：

(2a.2)桩侧单位面积的极限水平土抗力标准值P_u可按下列公式计算：

当Z<Z_r：

当Z≥Z_r：P_u＝9C_u

(2a.3)桩周土达极限水平土抗力之半时，相应桩的侧向水平变形Y₅₀可按下列公式计算：

Y₅₀＝ρε₅₀d

(2a.4)软黏土中桩的P-Y曲线，可按下列公式确定：

当Y/Y₅₀<8时：

当Y/Y₅₀≥8时：P＝P_u

其中，P_u：泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，单位：kPa；C_u：原状黏土不排水抗剪强度的标准值，单位：kPa；γ：土的重度，单位(kN/m³)；Z：泥面以下桩的任一深度，单位：m；ζ：系数，取0.5；D：桩径，单位：m；Z_r：极限水平土抗力转折点的深度；P：泥面以下Z深度处作用于桩上的水平土抗力标准值，单位：kPa；Y：泥面以下Z深度处桩的侧向水平变形，单位：mm；Y₅₀：桩周土达极限水平土抗力之半时，相应桩的侧向水平变形，单位：mm；ρ：相关系数，取2.5；ε₅₀：三轴仪试验中最大主应力差一半时的应变值；对饱和度较大的软黏土，也可取无侧线抗压强度q_u一半时的应变值；

(2b)砂土的数据处理计算过程如下：

(2b.1)依据《港口工程桩基规范》获得系数C₁、C₂、C₃以及土抗力的初始模量K；

(2b.2)砂土单位桩长的极限水平土抗力标准值P_u可按下列公式计算:

当Z<Z_r：P_u＝(C₁Z+C₂d)γZ

当Z≥Z_r：P_u＝C₃dγZ

(2b.3)计算系数Ψ，按下列公式计算:

(2b.4)泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，可按下列公式计算:

其中，P_u：泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，单位：kPa；C₁、C₂、C₃：系数，按《港口工程桩基规范》图D.2.3确定；Z：泥面以下桩的任一深度，单位：m；D：桩径，单位：m；γ：土的重度，单位(kN/m³)；P：泥面以下Z深度处作用于桩上的水平土抗力标准值，单位：kPa；Ψ：计算系数；K：土抗力的初始模量，按《港口工程桩基规范》图D.2.4确定，单位：kN/m³；Y：泥面以下Z深度处桩的侧向水平变形，单位：mm；φ：内摩擦角；

(三)：将步骤(二)得到的数据处理结果经过输出处理后即得到最终计算结果：P-Y文本文件和P-Y图形文件，供海上风机桩基基础设计人员使用。

本发明在砂土的输入参数中，内摩擦角的输入范围为[20,40]。

本发明实现P-Y曲线的自动化计算，计算速度快，并经过严密的验证程序，计算完整无误的同时计算速度得到了质的提高，人工成本节约1000倍以上；此外，本发明的计算方法简单快捷，操作方便，使用户易于学习和掌握，其高效的数值计算及符号计算功能，具有完备的图形输出功能，实现了计算结果和编程的可视化，能使用户从繁杂的数学运算分析中解脱出来。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明软黏土的数据处理计算流程图。

图2是本发明砂土的数据处理计算流程图

图3是无量纲化的软黏土P-Y曲线输出图。

图4是有量纲化的软黏土P-Y曲线输出图。

图5是砂性土的P-Y曲线输出图。

具体实施方式

如图1-图5所示，本发明包括以下步骤：

(一)：首先，根据地质勘察单位提供的地勘资料，填写六个参数并输入到MATLAB中调用：

(1a)软黏土包括直径、土层号、三轴仪试验中最大主应力差一半时的应变值、土的重度、原状黏土不排水抗剪强度的标准值、泥面以下桩的任一深度；

依据地质勘查报告获得土的物理性质参数：不排水抗剪强度标准值C_u、土的重度γ、泥面以下桩的任一深度Z、系数ζ、桩径D、相关系数ρ、泥面以下Z深度处桩的侧向水平变形Y、三轴仪试验中最大主应力差一半时的应变值ε₅₀；

(1b)砂土包括直径、土层号、土层位于水上/水下、土的重度、内摩擦角输入范围为[20,40]、泥面以下桩的任一土层深度；

依据地质勘查报告获得土的物理性质参数：内摩擦角φ、泥面以下桩的任一深度Z、桩径D、土的重度γ、计算系数Ψ、泥面以下Z深度处桩的侧向水平变形Y；

(二)：根据输入的参数对黏土或砂土进行数据处理：

(2a)软黏土的数据处理计算过程如下：

不排水抗剪强度标准值C_u小于等于96kPa的软黏土，在非反复荷载作用下P－Y曲线按下列规定确定：

(2a.1)极限水平土抗力转折点的深度Z_r，可按下式计算：

(2a.2)桩侧单位面积的极限水平土抗力标准值P_u可按下列公式计算：

当Z<Z_r：

当Z≥Z_r：P_u＝9C_u

(2a.3)桩周土达极限水平土抗力之半时，相应桩的侧向水平变形Y₅₀可按下列公式计算：

Y₅₀＝ρε₅₀d

(2a.4)软黏土中桩的P-Y曲线，可按下列公式确定：

当Y/Y₅₀<8时：

当Y/Y₅₀≥8时：P＝P_u

其中，P_u：泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，单位：kPa；Cu：原状黏土不排水抗剪强度的标准值，单位：kPa；γ：土的重度，单位(kN/m³)；Z：泥面以下桩的任一深度，单位：m；ζ：系数，取0.5；D：桩径，单位：m；Z_r：极限水平土抗力转折点的深度；P：泥面以下Z深度处作用于桩上的水平土抗力标准值，单位：kPa；Y：泥面以下Z深度处桩的侧向水平变形，单位：mm；Y₅₀：桩周土达极限水平土抗力之半时，相应桩的侧向水平变形，单位：mm；ρ：相关系数，取2.5；ε₅₀：三轴仪试验中最大主应力差一半时的应变值；对饱和度较大的软黏土，也可取无侧线抗压强度q_u一半时的应变值；

(2b)砂土的数据处理计算过程如下：

(2b.1)依据《港口工程桩基规范》获得系数C₁、C₂、C₃以及土抗力的初始模量K；

(2b.2)砂土单位桩长的极限水平土抗力标准值P_u可按下列公式计算:

当Z<Z_r：P_u＝(C₁Z+C₂d)γZ

当Z≥Z_r：P_u＝C₃dγZ

(2b.3)计算系数Ψ，按下列公式计算:

(2b.4)泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，可按下列公式计算:

其中，P_u：泥面以下Z深度处桩侧单位面积极限水平土抗力标准值，单位：kPa；C₁、C₂、C₃：系数，按《港口工程桩基规范》图D.2.3确定；Z：泥面以下桩的任一深度，单位：m；D：桩径，单位：m；γ：土的重度，单位(kN/m³)；P：泥面以下Z深度处作用于桩上的水平土抗力标准值，单位：kPa；Ψ：计算系数；K：土抗力的初始模量，按《港口工程桩基规范》图D.2.4确定，单位：kN/m³；Y：泥面以下Z深度处桩的侧向水平变形，单位：mm；φ：内摩擦角；

(三)：将步骤(二)得到的数据处理结果经过输出处理后即得到最终计算结果：P-Y文本文件和P-Y图形文件，供海上风机桩基基础设计人员使用。

本发明主要是通过MATLAB编程技术，节省海上风机桩基基础结构设计时间、节约成本、增加结构设计的可靠性。其自动动化实现是海上风电行业乃至整个建筑行业发展的必然趋势。本发明是基于MATLAB实现结构参数的处理功能及计算结果的后处理。它通过编写好的复杂应用程序，实现结构参数的计算，计算结果可直接应用于其它相关软件(包括有限元软件)。本发明只需简单提供几个参数，即可完成一次风机基础桩土水平相互作用P-Y曲线的自动化计算，即使在设计过程中，各单位间多次反复提资，只要调节这六个参数，可以很快就完成设计修改。

本发明总体思路为：依据软黏土和砂土的P-Y曲线计算原理，编制两个MATLAB程序命令流主程序和两个MATALB程序命令流子程序，其中，MATLAB程序命令流流程为：输入少数几个结构设计或土层参数。将图形计算参数(砂土)提取成数据对参数，并编制成文本文件。然后对所有计算参数，包括设计参数、图层参数和图形提取的数据对参数，按照《港口工程桩基规范》第D.2节规定，进行数据处理。运行两个MATLAB程序命令流主程序，输出处理后得到的计算结果——P-Y文本文件和P-Y图形文件，供海上风机桩基基础设计人员使用。与现属的现有技术相比，本发明有以下几方面优点：

(1)本发明直接以文本形式提供了P-Y曲线数据对，该项发明完成了前所未有的便于所有工程设计人员使用的工作。(2)手算完成一个6～8个土层的桩土水平相互作用P-Y曲线数据对的计算工作，原本需要一周的时间，使用本项发明可以在五分钟之内完成该项工作。(3)对于所有海上风电桩基基础工程，均需考虑桩土水平相互作用，该项发明在海上风电桩基基础工程的设计中发挥了关键性作用。上述步骤中，本发明通过MATLAB编程语言来实现此项技术，此外，本发明还可以采用其它编程语言来实现此项技术，比如C、C++、Fortran、Python等编程语言，均可实现桩土水平相互作用P-Y曲线的自动化计算功能。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。