Comput Geotech: 用ChatGPT写的MATLAB代码能做CFD？

ChatGPT（由OpenAI开发）是基于Transformer架构的最先进的大语言模型（LLM），用于生成类似人类的文本响应。ChatGPT能够接收对话提示，进行处理，并生成相关的文本响应。虽然还有其他的大语言模型，如BERT和T5（Devlin等人，2018年；Raffel等人，2020年），但ChatGPT通过其对上下文的细腻理解、对各种查询的适应性以及生成连贯响应的卓越能力脱颖而出。正因为如此，ChatGPT开创了人工智能领域的新篇章，并成为现代大语言模型（LLMs）所能实现的目标和对社会带来影响的代表。自2019年11月发布完整的GPT-2模型以来，关于ChatGPT在各个领域的应用已有广泛的报道和讨论，这些内容在表1中进行了简要总结。随着对ChatGPT的优缺点、潜在风险以及伦理问题的讨论，ChatGPT的应用领域似乎在迅速扩展，尽管其生成响应的可信度和真实性尚未被正式报道和公布。

表1. ChatGPT的已报道应用

人们认为，现在是时候设想、讨论并分享如何在科学和工程领域采用ChatGPT的强大功能，尤其是在这些领域中已有的文献参考数量有限的情况下。Bahrini等人（2023年）通过让ChatGPT生成系统与信息工程专业的初级考试问卷，评估了其在教育中的有效性。Tsai等人（2023年）提出了利用ChatGPT辅助的解决问题程序，用于开发化学工程教育课程，尤其是在建立蒸汽动力装置、相反应和传质模型的情况下。Kashefi和Mukerji（2023年）探讨了ChatGPT在生成和改进不同编程语言的数值算法代码方面的能力，应用于诸如泊松方程、扩散方程、不可压缩Navier-Stokes方程和可压缩无粘流动等各种数学问题。在建筑领域，ChatGPT不仅能够生成满足项目要求的施工进度表（Prieto等人，2023年），还能够促进危险识别和安全培训（Uddin等人，2023年）。

对先前研究的回顾得出的结论是，进一步探索ChatGPT在数值建模中的潜在应用是值得的。因此，本研究旨在调查ChatGPT在编写编程代码方面的能力，这是ChatGPT提供的最可靠功能之一，用于解决诸如渗流和边坡稳定性等著名且基础的岩土工程问题。同时，还尝试了对部分饱和砂的X射线CT图像进行图像处理。对于每个示例，首先定义了给定问题的概述，并通过逻辑提示向ChatGPT请求生成MATLAB编程代码。通过改进提示反复检查结果，直到获得正确的结果。通过这些用于自动化编写编程代码的指定示例，展示了将ChatGPT集成到解决传统岩土工程问题中的潜在应用。

选择了岩土工程中的三个基础问题，以探讨ChatGPT-4.0在通过提示生成MATLAB代码时的性能。以下是问题定义及后续优化过程的描述。

2.1. 二维渗流分析 

图1a中渗流的配置显示，土层高度为10米，宽度为25米，一道宽度为1米的板桩插入土层中，深度为5米。初始水位设置在板桩的左侧，以在两侧施加水力梯度。假设土壤是均质且各向同性的，土层底部被定义为不透水边界。板桩作为不透水屏障，允许渗流通过土层。

(a) 板桩周围的二维渗流分析

(b) 使用Fellenius切片法的二维边坡稳定性分析

(c) 通过图像处理对部分饱和砂的XCT图像进行相位分割

图1. 岩土工程中问题场景的示意图

拉普拉斯方程用于描述土壤域内水力头的分布。可以应用有限差分法（FDM）来获得近似解。为了实施FDM，将土层等效离散为在x方向上间隔为dx、y方向上间隔为dy的正方形网格，如图1a右侧网格化框中所示。通过中心差分法来近似拉普拉斯算子，第i行j列节点的水力头可以通过以下公式获得：

基于水位、初始条件和边界条件，可以迭代计算每个节点的水力头，直到收敛为止。最终水力头值的分布可以用于描绘等势线和流线。使用ChatGPT生成的代码所得结果与使用商业软件GeoStudio SEEP/W得到的结果进行了比较。

2.2. 边坡稳定性分析

边坡稳定性分析的关键参数包括边坡几何形状，如高度（h）和边坡倾角（β），以及土壤性质，包括黏聚力（c）、单位重度（γ）和摩擦角（ϕ）。边坡的稳定性通过安全系数（FS）来量化，FS定义为可用抗剪强度与作用在潜在破坏面上的剪应力之比。需要注意的是，最小化FS的潜在破坏弧及其位置指示了最脆弱的区域，称为“临界破坏面”。

Fellenius切片法，又称普通切片法，被用于给定边坡几何形状的边坡稳定性分析（Fellenius, 1936）。该方法假设潜在破坏面为圆形滑动面，并将可能滑动的土体分割为多个垂直切片以评估FS。它还假设切片间的作用力平行于每个切片的基底，且作用在侧面的力在破坏面上的法向为零。因此，FS可以仅基于力矩平衡推导出来，从而进行确定性的边坡稳定性分析。在此分析中，假设孔隙水压力为零。潜在破坏面假设为如图1b中红色弧线所示的圆弧，其中心坐标为C (x, y)，半径为r。作用在每个切片基底（如蓝色阴影梯形所示）上的法向力（Nn）和剪切力（Tn）为：

其中，Wn 是第 n 个切片的重量，αn 是第 n 个切片基底的倾角。基于力矩平衡的安全系数（FS）定义为相对于中心 C 的每个切片的驱动力矩（Md）之和与抗力矩（Mr）之和的比值。第 n 个切片的驱动力矩、抗力矩和 FS 可以通过以下公式确定：

通过改变潜在破坏弧的中心和半径，可以计算出不同情况下的安全系数（FS），最终确定最小的FS及其对应的几何形状。ChatGPT实现了上述过程，并使用商业软件GeoStudio SLOPE/W对结果进行了验证。

2.3. 相位分割的图像分析  

部分饱和土壤包括三种相态，即固体颗粒、水和空隙（空气），其切片通过二维X射线计算机断层成像在8位灰度下可视化，如图1c-i所示。A-B段的CT值轮廓线如图1c-ii所示，可以大致勾勒出相邻相态之间逐渐变化的边界。图1c-iii中的强度水平（即CT值）直方图清楚地表明，空气-水界面和水-固体界面的边界值分别位于大约91和113附近。使用传统的分割方法时，由于水的CT值与界面处逐渐变化的CT值相似，水不可避免地会在固-气界面附近被误分割。

因此，从灰度图像中分割多个相位不仅需要采用传统的阈值分割方法，还需要使用先进的图像处理技术，如迭代膨胀-侵蚀法（Adams和Bischof, 1994; Serra, 1982）。在本研究的范围内，目标是仅基于可感知的提示逐步开发由ChatGPT驱动的数值代码，而无需明确提供图像处理算法。

2.4. 通过使用ChatGPT生成代码解决问题的过程  

图2概述了问题解决的工作流程。对于每个案例，通过描述几何形状、初始条件、边界条件和计算方法来清晰地定义问题。需要注意的是，提示由叙述性描述组成，不包含任何数学公式，从而允许ChatGPT根据其理解自由生成数值代码。在请求ChatGPT生成MATLAB代码后，生成的代码被执行。一旦代码运行无误，用户会仔细审查代码和后续的结果以进行验证。根据此审查，用户会提供额外的提示来纠正代码，直到它完全没有错误。这些提示会不断被细化和优化，以确保ChatGPT生成的结果与预期的正确答案一致。

图2. 通过使用ChatGPT进行自动代码生成的解决问题流程图

3.1. 二维渗流分析  

图3展示了首次尝试为二维渗流分析生成MATLAB代码的过程。图3的上框中包含了全面的问题定义，如几何细节、初始和边界条件、空间域网格配置的规格以及迭代计算的收敛标准。需要注意的是，虽然提示明确要求使用有限差分法（FDM），但并未详细说明该方法的具体方程或操作原理。图3下部分显示了由ChatGPT生成的结构合理的MATLAB代码。为了清晰起见，生成的代码中每个部分基于用户提供的提示被用不同颜色标注。

图3. 第一次向ChatGPT提供的渗流分析提示（上框）及ChatGPT生成的相应MATLAB代码（下框）

在第一次尝试中，用户定义的特定尺寸变量正确地反映在第21行。第24到29行（以绿色高亮）建立了土层的分析域、初始条件以及有限差分法（FDM）应用的收敛标准。土层根据用户定义的间隔dx和dy被划分为网格，每个网格点被分配给矩阵变量H，表示总水头值。该矩阵H根据初始条件初始化，反映了板桩两侧的平均水位。随后的代码部分（第32-34行，以黄色高亮）定义了板桩在土壤域中的位置，并通过为两侧的水位分配常量水头值H1和H2来建立边界条件。利用用户指定的阈值作为FDM的收敛标准，ChatGPT自主构建了第37-52行（以蓝色高亮）的FDM迭代计算代码，该代码集成了while循环、for循环和if语句。值得注意的是，在这个迭代部分中，第43-45行包含了一个排除板桩内计算的条件语句。为了简洁起见，绘制流网的部分被省略。第一次尝试生成的MATLAB代码运行无误，ChatGPT生成的流网如图4a所示。GeoStudio SEEP/W生成的流网结果也显示在图4b中。对比这两个图可以发现，关键差异源于第32行中与水位相关的不正确边界条件语句。在矩阵H的第一行中，板桩左侧应指定为高水位H1 = 10米，而右侧应为低水位H2 = 0米，作为常量水头值。

图4. 通过ChatGPT生成的MATLAB代码在第一次尝试中分析的流网(a)与在相同几何和边界条件下由GeoStudio SEEP/W分析的基准结果(b)的比较

第二次尝试的提示如图5所示，旨在纠正矩阵H的边界条件，代码也在第14-15行通过正确分配常量水头值进行了修正。然而，这一修正并未产生正确的流网，如图6a所示。原因是“越界”情况下缺少了用于特定节点（如AF、FE、ED、CH、HG和GB线上）的有限差分法（FDM）计算所需的四个相邻值。因此，第三次尝试的提示详细解释了常用的“镜像”技术，以用于那些缺少相邻值的节点的FDM计算，这部分内容在图5中以蓝色高亮显示。首先在板桩区域分配了“NaN”值（第18-19行），然后检查“NaN”值并正确实施镜像技术（第30-33行）。

图5. 为解决错误向ChatGPT提供的第二次和第三次尝试提示（上方两个框），以及ChatGPT生成的相应修订后的MATLAB代码（下方框）

(a) 第二次尝试结果                      (b) 第三次尝试结果

图6. 通过ChatGPT生成的MATLAB代码分析的流网

通过迭代改进提示，最终生成的流网（图6b）使用ChatGPT自动生成的MATLAB代码，与基准流网准确匹配。尽管SEEP/W基于有限元法（FEM），但两种方法生成的矩阵H中的水头值比较后显示出差异极小。

3.2. 边坡稳定性分析

第一次尝试的提示指定使用Fellenius方法来计算安全系数（FS），并包括变量的声明，如边坡几何形状、土壤性质、坐标和表示潜在破坏面的弧的半径（图7）。需要注意的是，虽然建议使用Fellenius方法，但其基本原理并未解释。ChatGPT生成的代码以用户自定义函数“Fellenius”的形式显示在图7的下部分。从第6行到第16行，生成了滑动体的坐标、滑动体切片的代码，以及每个切片重量的计算。第19行到第21行定义了驱动力和抗力的基本方程，以计算FS。

图7. 第一次向ChatGPT提供的用于创建计算FS的MATLAB函数的提示（上框）及ChatGPT生成的相应MATLAB代码（下框）

初始代码虽然没有语法错误，但包含逻辑错误，如图8所示。计算边坡的安全系数（FS）需要识别滑动体、将其分为N个切片，然后基于每个切片的重量计算驱动力和抗力。错误情况1表明滑动体识别不正确，具体是由于边坡长度范围设置错误（红色箭头显示的错误定义长度和黑色箭头显示的正确长度）。在错误情况2中，分成N个部分的切片被视为矩形，且每个切片的高度错误地设置为第n个y坐标值（红色箭头所示）。这应通过考虑边坡的y坐标与切片基底的y坐标之间的差异进行修正，如黑色箭头所示。错误情况3中，用于计算FS评估的切片基底倾角的变量错误地设置为与边坡相关的倾角。

图8. (a) 三种不同破坏类型的分析域示意图，以及使用ChatGPT在第一次尝试中生成的MATLAB代码计算的FS与GeoStudio SLOPE/W结果的准确性评估比较。(b) ChatGPT在第一次尝试中创建的MATLAB函数中的三个错误导致了不正确的FS计算结果。

为了纠正这些错误，向ChatGPT提供了额外的提示以指导代码改进，修订后的MATLAB函数代码如图9所示。错误1被修正为确保所有部分都包含在边坡长度内（第9行）。为纠正错误2，ChatGPT添加了一个for循环，根据每个切片在边坡长度内的位置不同来计算其高度（第16-24行），并考虑到梯形的形状。对于错误3，替代了在确定边坡几何形状时提供的倾角变量（beta），在第29行引入了一个新定义的角度（alpha），以考虑未定义的第n个切片的倾斜度。这个变量描述了从破坏圆心到第n个切片基底的线与垂直线之间形成的角度。

图9. 为解决三种不同错误而向ChatGPT提供的第二次尝试提示（上框）及ChatGPT生成的相应修订后的MATLAB代码（下框）

该代码在第30-32行使用角度变量alpha计算驱动力和抗力，从而得出FS值。为了简洁起见，主函数（使用for循环和MATLAB函数代码来计算潜在破坏圆每个坐标的FS值）被省略。

对于具有不同半径的潜在破坏圆心范围，ChatGPT得到了最小的FS值为1.630（图10a）。热图表示了给定坐标下计算得到的FS值的分布。这个结果与GeoStudio SLOPE/W的结果相同（图10b）。在这两种情况下，趾部破坏类型的临界破坏圆的坐标x、y和半径r分别确定为-0.30、14.78和14.77，对应的最小FS值计算为1.630。

图10. 使用热图可视化潜在破坏面的FS计算结果，并展示临界破坏面。结果分别通过(a) ChatGPT创建的MATLAB函数和(b) GeoStudio SLOPE/W获得。（边坡尺寸：h = 10米，β = 40°，土壤性质：Ф = 35°，c = 12 kN/m²，γ = 17 kN/m³）。

3.3. 相位分割的图像分析  

图11展示了对部分饱和砂的X射线CT图像中三种相态（即颗粒、水和空气）进行分割的第一次尝试。每个相态的CT值范围被大致提供，并在提示中指定了分割相态的颜色。通过使用名为Link Reader的ChatGPT插件，上传了8位灰度图像。ChatGPT生成的代码与提示中的输入几乎相同，首次输出的图像显示在图11底部。需要注意的是，为了简洁起见，加载图像的代码被省略。由于部分体积效应（PVE），使用第一次尝试中编码的传统阈值方法无法避免在颗粒和空气边界周围捕获到的不必要的环状水相。由于空间分辨率的限制，PVE往往会导致相位边界处不同物体的模糊和融合（Glover和Pelc, 1980; Herman, 1980）。

图11. 向ChatGPT提供的用于生成执行相位分割的MATLAB代码的第一次尝试提示（上框）及ChatGPT生成的相应MATLAB代码（下框）

因此，提供了第二次提示，告知ChatGPT部分体积效应（PVE）的起源以及所需的输出图像，而未指定相关算法。图12展示了第二次提示、ChatGPT生成的MATLAB代码以及输出图像。ChatGPT自愿提出了腐蚀和膨胀的方法，这确实已被广泛应用于克服PVE（Haralick等，1987；Serra，1982）以去除伪影，正如ChatGPT在代码之前的回答中准确解释的那样。ChatGPT最初设置的 `strel` 函数和跨度值 `3` 可以由用户调整以产生所需的结果。然而，这种方法往往过度去除了环状水相，导致饱和度从59.03%降低到51.66%，如图12底部的输出图像所示。第一次尝试中捕获的被去除水相的轮廓在图12底部接触区的高亮图像中以红色呈现。

图12. 向ChatGPT提供的用于生成执行相位分割的MATLAB代码的第二次尝试提示（上框）及ChatGPT生成的相应MATLAB代码（下框）

在第三次尝试中，提示建议使用替代策略来对抗PVE，例如“边缘保留平滑”或“平滑约束”，这些方法已被用于在减少噪声的同时保留PVE引入伪影的重要结构细节（Rudin等, 1992；Smith和Brady, 1997；Tomasi和Manduchi, 1998）。图13展示了第三次提示、ChatGPT生成的代码以及结果。这导致ChatGPT使用了MATLAB的 `imbilatfilt` 函数（第15行），该函数基于像素的接近度和强度相似性执行非线性和边缘保留平滑。由于这种方法特别补偿了粒子接触区域附近的水相，饱和度略微回升至52.83%。

图13. 向ChatGPT提供的用于生成执行相位分割的MATLAB代码的第三次尝试提示（上框）及ChatGPT生成的相应MATLAB代码（下框）

本文探讨了在岩土工程中利用ChatGPT辅助编程来解决渗流、边坡稳定性和X射线图像处理等基本但重要的问题的可行性。对于每个案例，通过对话式提示向ChatGPT提供了问题定义、初始条件和边界条件，而无需直接输入数值或数学定义与方程。对初始提示生成的MATLAB代码进行了校对并执行，以验证结果的准确性。通过为需要进一步修改的代码提供附加提示，代码被逐步更新，直到能够生成准确且可接受的结果。提出的用于自动化编程任务的ChatGPT-MATLAB框架引出了以下几点观察：

研究结果进一步表明，ChatGPT并不能替代人类的判断力和专业知识，而是可以加速和简化编程过程的工具。用户必须具备基础的学科知识和编程原理，才能有效地利用ChatGPT。本质上，ChatGPT能够加速代码生成的复杂过程，但并不能完全替代它，仍然需要人类的监督和伦理考量，以确保结果的准确性和相关性。

在本研究中展示的借助ChatGPT进行数值编码的案例，仅是生成式人工智能在工程、研究和教育领域广泛潜在应用中的一个例子。预计这一技术的能力将扩展到处理更复杂的工程问题，例如多过程耦合建模、设计任务以及大数据分析。生成式人工智能的潜在应用及其在科学和工程领域的作用，是应该由我们社会共同讨论的问题。更广泛的社区应参与这些可能性的讨论，因为集体的经验和知识将提升生成式人工智能在研究和教育中的应用效果。

翻译转载自：《Computers and Geotechnics》: "A ChatGPT-MATLAB framework for numerical modeling in geotechnical engineering applications"