基于轴对称接触问题的有限单元法及其在桩基工程中的应用

'西南交通大学硕士学位论文轴对称接触问题的有限单元法及其在桩基工程中的应用姓名：王新申请学位级别：硕士专业：岩土工程指导教师：毛坚强20031001 西南交通大学硕士研究生学位论文第1页摘要桩土之间的接触面对桩基的受力特性有着重要的影响，在计算中不应无视它的存在。有限元等数值方法虽然已在桩基特性分析中得到了广泛的应用，但计算中对桩士接触面的处理却仍然是一个并未很好解决的问题。本文的主要目的就是将该类问题作为接触问题进行处理，然后建立一种桩土相互作用的计算方法，并应用于桩基的工程实践之中。在有限元中计算中，常用的方法是利用节理单元来模拟接触面。本论文对这种方法经过分析之后，发现这种处理方式对接触面的描述并不合理，也无法反映出接触面的主要力学特征。同时指出，将桩土相互作用按接触问题来处理，可以更好地反映出该接触面的主要力学特征，即：(1)接触面具有相应的“抗剪强度”，即当接触面上的某一点处的切向应力小于该点处的抗剪强度时，该处桩土共同变形；达到其抗剪强度时，则沿切向发生滑动；(2)在整个变形过程中，桩与土不能互相侵入。针对于接触问题的求解，目前已经发展了很多计算方法，如数学规划法、罚函数法等。本文采用的是[2]中提出的“自由度法”，并根据其计算原理建立了轴对称等参数单元的有限元计算公式，并详细介绍了弹塑性接触问题的有限元求解实施过程。“自由度法”涉及到的概念与普通的接触有限元相比，既有同也有异，因此作者在阐述本法相关概念的同时，还与普通接触有限元的处理方式作了对比。传统的有限元程序采用结构化的语言编写，由于结构化语言本身的缺陷，使有限元程序的更新、扩展功能受到很大的限制，80年代出现的面向对象编程思想为解决这个问题提供了转机。目前，已有学者将面向对象编程思想引入到有限元编程当中，但在接触问题中，却研究较少。本论文中，作者将面向对象编程的思想应用于接触问题有限元程序的编制，并建立起了轴对称接触问题的对象类框架。完成上述工作后，作者对西昌某工地现场的两根试桩进行了计算分析。通过计算，得到相应的荷载一沉降曲线、轴力图、桩侧阻力图、塑性区、土的竖向位移等值线、桩的滑移过程等结果，并与现场试验结果吻合较好，说明用本文方法计算分析桩基的工作特性是可行和有效的。关键词接触问题；有限元；桩基；面向对象编程 西南交通大学硕士研究生学位论文第ll页Absl：ractAstheinterfacebetweenpileandgeoteclmicalmediae腩ctsthemechanicalbehaviorsofthepilefoundationsignificantly，areasonablemodeliscriticalinthecalculationandanalysisofpile．Sofar,thoughfiniteelementmethodhasbeenwidelyusedinpilefoundationengineering，itstillremainsanintractableproblemtosimulatetheinterfaceideallyinnumericalanalysis．Theauthordeemsthatitisapracticalwaytotreatthisproblemascontactprobleminsolidmechanics．Themainintensionofthispaperistodevelopafiniteelementapproachforaxisymmetriccontactproblemsandapplyittopilefoundation．Sofar,intheFEManalysis，thecontactinterfacesusuallyaretakenasthejointelement，whichwasintroducedbyGoodmanearlyin1968．Nevertheless，inthemodel，thecontactinterfacesarecharacteredsimply诵t11stiffnessKnandKs，whickisnotperfectandrigorousintheviewofsolidmechanicsandratherroughforthedescriptionofthebehavioroftheinterfaces．Incontrast，takenasacontactproblem，theprincipalcharactersoftheinterfacescanbereflected：(1)thecontactinterfaceisofshearstrengthwhichimpliesthatagenericpairofcontactpointsontheinterfacewilldisplacecoincidentallyifthetangentialforceatthispointhasnotreachedthelimitresistance，else，relativeslidebetweenthepointswillOccuralongthetangentialdirection．(2)NomattertherelativeslideOccursornot,thecontactbodiescannotinvadeeachotherinanywise．Thereareconsiderablenumbersofresearchconcerningfiniteelementmethodoncontactproblemandalotoftechniquessuchasmathematicsprogramming，penaltymethodandsoon．Themethoddependsonthe‘FreedomDegree’technique，whichwasdevelopedin【2]，isaneffcctivemethodtosolvecontactproblems．InthisPaper，onbasisofthetheory,theanisymmetricisoparametricFEMformulaewereconstructedwiththeprincirlleofvirtualwork．Moreover,theimplementationoftheseformulaeonthecomputerisalsodiscussedindetail．Also，thesimilarityanddifferencebetweenthemethodandconventionalmethodalecompared．Mostofthepresentfiniteelementprogramsarcwritteninthestructurallanguage，suchasFORTRAN，Candsoon．Asthesizeofprogramsgrows，thedefectsofthestructurallanguage，whicharchardtomaintainandupdate，areexposed．Theobject—orientedprogramming(OOP)thatappearedin1980’scansolve 西南交通大学硕士研究生学位论文第11I页thisproblemwell．Sofar,thoughOOPhasbeenappliedtothefieldofFEMprogramming，itrarelybeusedtosolvecontactproblems．Inthepaper,theprogrammingofcontactFEMCOdewithtechniqueofOOPiSdiscussedandtheprogramframewi也object-classisestablished．Lastly,twoaxiallyloadedboredpilesinsomebuildingsitenearXichangcitywascalculatedwiththemethodandsomeusefuldataandcurvesincludingload—seRlementcurvesanddistributionoffrictionresistanceetc．arepresented．Theaccordancebetweentheresultsofcalculationandin—situtestverifiesthevalidityofthemethoddevelopedinthePapel-．Keywordscontactproblem；finiteelement；pile；object—orientedprogramming 西南交通大学硕士研究生学位论文第l页第1章绪论1．1接触问题有限元法的国内外研究现状桩基由于具有承载力高、沉降量小而均匀、能适应各类地层条件等优点，在工程中得到了广泛的应用⋯，然而目前桩基的计算理论和方法却落后于桩基的工程实践。一般来说，各类工程结构的形式虽然不同，但其计算过程主要是两项内容：(1)针对所研究的问题建立相应的力学计算模型；(2)应用适当的数学方法进行求解。桩基的计算也不例外，为了得到满意的计算结果，桩基的计算模型应尽可能反映出桩基的主要力学特点。例如，在计算分析桩的承载特性时，对桩一土间接触面的力学特征的准确描述就十分重要，这也是建立桩一土作用计算模型的一个关键问题。计算模型建立后的工作就是将建立的模型在计算中实现。我们知道，数值方法特别是有限元法的发展和应用，使得很多复杂力学问题的求解成为可能。目前，若不考虑桩土间的接触面而将桩和土作为连续体，其计算与一般的固体力学问题并无本质的区别，相应的弹性、塑性、粘性等计算方法已相当成熟。但实际上，对于桩基而言，在荷载传递过程中，桩土之间会沿接触界面发生相对滑移，如果还用连续体的方法，则与桩的实际受力特点相差较远，计算结果也会相差较大。在有限元中常用的处理方法是用接触面单元来模拟桩土间的接触面，但由于对接触面力学性质的过度简化，使该类单元存在着难以弥补的缺陷H1(详见2．2．1节)。从本质上讲，对桩土间的接触面的模拟应属于固体力学中的接触问题。虽然有限元等数值方法的应用使求解复杂问题的能力大大提高，但对接触问题的求解至今仍是一个难点，所发展的种种求解方法，迄今为止并没有那种非常令人满意。但从目前计算力学发展水平看，将桩土间相互作用作为接触问题来求解是一个比较合理的途径。在文献【2]中提出了一种接触问题的计算方法并成功地用于岩土工程问题的求解，本文将在此基础上，发展一种计算桩一土相互作用的有限元计算方法。1．1．1非线性有限元的研究概况有艰元的基本恩想，最早可以追溯到1943年，Courant[41应用定义在三角形区域上的分片连续函数和最小位能原理相结合，求解了st．Venant扭转问题。但 西南交通大学硕士研究生学位论文第2页此后，由于各种原因，有限元的概念没有得到大家的重视。直到20世纪60年代以后，随着电子计算机的广泛应用和发展，有限元的发展才显著加快。现代有限元的第一个成功尝试，是将刚架位移法推广应用到弹性力学问题，这是Turner、CloughpJ等人在1956年在分析飞机结构时得到的成果。他们首次给出了用三角形单元求解平面应力I"ol题的正确解答。1960年Clough[61在平面应力问题的基础上，进一步处理了平面弹性问题，并第一次提出了“有限单元法”的名称。轴对称有限单元法首先由E．L．Wilsonl7『1在1965年发表，并且很快就有了相应的计算程序。最先将有限元应用于土力学的学者是Clough，1966年，他与Woodward用有限元方法做了土坝的非线性分析【81。近几十年来，有限元法的理论和应用都得到了迅速、持续不断的发展。首先，有限元的理论基础在进一步的深入。Turner、Clough等人开始提出有限单元法时利用的是直接刚度法，单元的特性分析也是由弹性理论方程来确定的。直接刚度法来源于结构分析的刚度法，力学概念清楚，容易理解，但只能处理一些比较简单的实际问题。1963～1964年，Besseling、Melosh和Jones等人证明了有限单元法是基于变分原理的Ritz法的另一种形式，从而可用Ritz法分析的所有理论基础都适用于有限单元法，确认了有限单元法是处理连续介质问题的一种普遍方法pJ。利用变分原理建立有限元方程和经典的Ritz法的主要区别是有限单元法假设的近似函数不是在全部求解域而是在单元上分片连续的。从20世纪60年代后期，人们进一步利用加权余量法来确定单元特性和建立有限元求解方程。目前，有限单元中所利用的主要是Galerkin法，它可以用于已知微分方程和边界条件、但泛函尚未找到或者根本不存在的微分方程的求解，因而进一步扩大了有限元的应用领域。有关有限元的数学基础，详见文献[10～11]。非线性问题一直是有限元分析中的一个热点和难点，对固体力学问题，三种非线性类型为：。(1)材料非线性。包括非线性弹性、弹塑性、粘弹塑性等。(2)几何非线性。即计算时需计入几何方程中的高阶微量部分。(3)边界非线性。如果边界条件随应力、应变等量的变化而变化，则构成边界非线性问题。接触问题就是典型的边界非线性问题。实际问题往往可能是上述三种基本非线性问题的组合。本论文所研究的桩土相互作用问题就同时考虑了材料非线性和边界非线性。组合后的问题，要比单一的非线性问题更为复杂。接触问题与经典问题的主要区别在于存在着接触边界。接触边界对整体求 西南交通大学硕士研究生学位论文第3页解结果影响较大，因此如何将接触边界代入到整体平衡方程当中又不增加很大的计算求解难度，成为接触问题求解的难点。本文将在下一节详细介绍目前主要的求解方法。1．1．2接触问题求解的研究概况固体力学主要研究变形体对外部作用的响应。事实上，除重力、电磁力等场力外，对一个物体施加荷载的过程，就是施加者、被施加者之间相互接触、作用的过程，只不过当满足一定条件时，我们将其中的一方简化为明确的、不变的作用，成为应力或位移边界，而研究另一物体的受力和变形，此即为经典固体力学的主要研究内容；当一方不能或不宜化简为应力或位移边界条件而必须同时考虑各物体的受力和变形时，即为我们所说的接触问题【2】。如上节所述，接触问题是典型的边界非线性问题。在接触问题中边界条件通常不能直接给出，它随接触过程的发展而变化。简单的弹性接触问题自上世纪末就由Hertz研究过，即著名的Hertz问题。Hertz求解接触问题时，对接触面的性质作了如下简化：(1)接触面是光滑的，即接触面上只有法向作用力，而无切向力。(2)预先假定法向应力具有某种特定的分布形式。由此可得到某些简单问题的解析解，且应用比较方便，因此一直受到重视，直到今天，仍有广泛的研究和应用i14’1”。但Hertz问题的研究对象只是诸如球体、圆柱体等简单形状物体的接触，且对接触面的性质作了很大的简化，因此，虽然该法至今仍在发展，但也远不能解决实际工程中存在的各类复杂的接触问题。实际上，即使对一般的弹性力学问题，能够利用解析法求解的问题也十分有限。有限元法等数值方法及计算机的出现才使这一问题的研究有了长足的进展。求解接触问题，一般包括以下四个方面的内容：(1)物理模型：采用什么样模型来描述两接触体之间力的传递以及在不同荷载作用下接触状态的变化。(2)几何运动规律：在接触面上两物体位移必须满足的条件。(3)本构规律：在接触面上力与位移或压力与切向力之间的关系。(4)建立方程与求解的方法：如何建立数学方程描述以上各规律并使系统处于平衡以及如何求解该方程。以上研究内容中，建立方程与求解方法是当前的研究热点，总结起来目前主要有以下几种：(1)罚函数法。它的基本思想是将两接触体在接触时材料不得相互侵入的 西南交通大学硕士研究生学位论文第4页条件作为惩罚项引入接触方程。有限元计算时对已知位移边界条件采用的“主元赋大值法”，实际上就是应用了罚函数的技巧。但这种解法的精度强烈依赖于所用罚因子的大小，过小则起不到作用，过大有可能导致求解方程的病态。(2)Lagrange乘子法。该法将接触边界条件用Lagrange乘子相乘并与总势能一起构成修正的势能，再求驻值得到最后的求解方程。若将Lagrange乘子取为无穷大，则该法变为罚函数法。该法得到的有限元计算公式中，与Lagrange乘子对应的刚度矩阵对角元素为零，由此可能造成总刚的奇异性而导致无法求解，此外还有增广的Lagrange乘子法¨⋯。(3)数学规划法。接触问题的数学规划法是基于势能或余能原理，利用变分不等式等现代数学方法导出的，理论上比较严格和直观。最初该方法是针对无摩擦接触问题提出的09]，它利用了无摩擦接触问题的非穿透条件和互补条件：g。=U。一u。’一《≥0(1-la)(U—go’一钟)Po=0(1-lb)只≤0(1-lc)上式中的U，、u。。、d?分别为物体力及D’的接触边界上，任一接触结点对的法向位移及其初始间距，只为其法向作用力(切向力为O)，其中，第一式表示接触体不能互相侵入，第二、三式表示法向作用力只能为压力(接触时)或0(分离时)。数学规划法解有摩擦接触问题的关键在于将摩擦条件表示成互补形式。一种方法是利用凸分析理论，把摩擦条件可以写成如下带导数的互补形式【20】；另～种处理方法是引进惩罚因子，然后仿照塑性力学将摩擦接触条件表示成有惩罚因子的互补形式。有了上述摩擦接触条件的互补关系，就可以利用参变量变分原理或虚功原理建立摩擦接触问题的有限元二次规划(线性互补)模型。数学规划方法在弹塑性接触问题的应用，通常用迭代法反映材料非线性特征，在每次迭代中用数学规划方法求解接触问题口”。然而，钟万勰等【22】利用参变量变分原理将接触问题和弹塑性问题表示成具有相同形式的有限元参数二次规划问题，很方便地实现了弹塑性接触问题的数学规划解法。(4)杂交或混合有限元法。PianTHH等将接触面上连续条件作为一个约束并把接触力处理成附加场变量的方法引入余能泛函中构成修正余能原理，从而导出应力杂交接触有限单元，但最后可以导出仅以位移为未知变量的单元矩阵【23甾]j将接触压力与位移作为未知数用摄动Lagrange乘子法(PerturbedLagrangianMethod)导出混合接触有限单元也是这一方法的另一种途径。但该类有限单元中，单元多为简单单元，该法的成功与否关键取决予与假定位移多项式联系的压力函数的选择是否恰当，而这点常常是比较困难的。 西南交通大学硕士研究生学位论文第5页1．2基桩竖向工作特性分析方法简介基桩的竖向工作特性是桩基础的重要研究内容之一。其主要研究方法可分为以下四类：1．试验方法包括现场试验和模型试验两种。现场静载试验是传统的也是最可靠的确定桩承载特性的方法，它不仅可确定桩的极限承载力，而且可以通过埋设各类测试元件获得桩身轴力、桩侧阻力、桩端阻力等诸多信息。由于试验费用、工期、设备等原因，往往只能对部分工程的少量桩进行试验。桩动测法是近年来发展起来的一种快速测桩方法，但目前该技术的成熟程度尚不足以取代静载试验。与现场试验相比，模型试验更易控制试验条件，但难以完全真实地反映桩的工作特性，而且岩土的相似准则也不好确定。2弹性理论法以弹性理论为基础的一些计算方法，各种解法略有不同，但一般都是基于桩的位移与邻近土位移之间的协调条件，只是假定桩长方向桩侧剪应力分布形式的不同而己，本文在第四章有详细介绍。弹性理论法不可能解决桩一土体系的的非线性问题，且难以考虑成层土问题，故该法在应用上受到较大的限制，工程实践上也受到一定的局限性。3．简化模型法该法将桩周土和桩端土简化为一系列的等效弹簧，其力学系统为桩弹簧系统。桩身一点的平衡微分方程为：EAd2u／d2x一曲(Ⅳ)=0(I-2)E一混凝土的弹性模量；4一桩横截面面积；“一桩在深度z处的位移；q俐一桩侧摩阻力的分布规律(荷载传递函数)，s为桩身截面周长。另外根据实测资料，桩底土可用一个非线性弹簧来代替，即可假定桩底土反力p倒为：P=f(u、(1—3)边界条件为：EAdu／djdr=，=一P(1-4)由方程(1-2)、(1-3)、(1-4)可确定桩在荷载作用下的荷载一沉降曲线。简化模型法最主要的两种简化模型为：传递函数法和剪切变形传递法。两种方法的不同之处在于对弹簧特征、描述方法的不同以及求解(1-2)的方法上。 西南交通大学硕士研究生学位论文第6页求解(1-2)必须知道q(u)的表达形式，cl(u)曲线不仅与土的自身特征有关，还和桩的埋置深度、桩径有关。在传递函数法中，各个学者利用不同的拟合曲线性式来对试验数据拟合，从而得到q(u)的形式，如viiayergiya法、佐藤悟法、陈竹昌法、双折线法等，因此从本质上传递函数法是一种经验方法。而剪切变形传递法则是根据模型的假设，从理论上得到了桩表面位移与桩土表面剪切应力的关系，该关系式实际上是q(u)的隐式表达，从理论上讲，我们可以从该关系导出q(u)的表达。4．数值解法数值解法是随着计算机的发展而发展起来的，目前应用最多的是有限单元法。通过有限元计算，不但可以得到荷载．沉降曲线，而且可以确定桩周岩土的应力场、位移场、塑性区等其它方法难以得到的结果。采用有限元法，不同的研究者出于对不同的土性、桩型、研究内容及其它考虑，因此采用了不同的具体的方法，但总的看来，主要反映在两个方面：(1)单元类型和计算方法的不同；(2)岩土的本构关系。Potts“”对置于正常固结粘性土中的单桩，考虑各种排水条件，分析了土中的有效应力，计算按轴对称问题考虑，采用的是八结点等参数单元，土的本构关系采用修正的剑桥模型：E1lison071等研究伦敦粘土中钻孔灌注桩的荷载传递机理时，按轴对称问题考虑，采用的是三角形单元、矩形单元及连接点的弹簧组成的组合模型，土的本构关系采用的是实测的三折线模型：Trochanis“”在应用三维有限元分析矩形截面单桩和群桩时，桩和土采用的是27结点三维矩形等参数单元，桩土界面采用的是18结点二维接触面单元，土体根据位移条件按线形弹簧或Drucker—Prager模型考虑；Desai汹1在研究砂土中桩土体系性状时，采用矩形等参单元，土的本构关系采用Duncan-Chang模型，接触面应力应变关系采用的是双曲线模型。但上述分析都在桩土之间设置了接触单元，通过对接触单元来模拟桩土之间的接触面，但接触单元由于自身的缺陷，使其在应用上受到了一定的限制。近些年来，也有学者利用接触问题的处理方法来研究桩基，但大都想是通过建立一种特殊的接触面单元，通过这些单元来模拟接触面，但构造出来的单元大都有一定的限制，通用性较差。本文采用文献[2]中的“自由度法”，它没有直接去构造接触单元，而是直接对接触点对自由度的处理来实现接触问题的求解。1．3接触问题常用计算方法简介1．3．1常用接触面单元简介 西南交通大学硕士研究生学位论文第7页岩土与结构共同作用时，接触面上的受力和变形可能会呈现两种状态：一是岩土与结构之间只有力的传递，没有相对位移，也就是没有错位或拉开。计算时可以看成由两种材料组成的连续体，进行有限元计算时不会产生任何困难；另一种情况是岩土与结构之间发生相对位移，从整体上来说是不连续的，为了进行有限元计算，就要设置所谓接触面单元来处理这种不连续性。目前，岩土工程中常用的接触单元有两种：Goodman单元和薄单元。两种单元对接触面的力学行为考虑不同，Goodman单元考虑了两相介质接触界面间的位移不连续性；而薄层界面单元则相反，通过结点位移对单元内部位移场的插值构造，认为两相介质接触界面间的位移有某种连续性。在本文的公式中，以黑体字表示矢量、张量、张量运算符号等，以后出现时不再一一说明。1．Goodman单元Goodman等人1968年提出的岩石节理单元[3们，被广泛用作接触面单元。这种单元为无厚度的四结点单元。如图1—1所示，上、下接触面之间由无数的法向和切向的微小弹簧相联系。相应的应力为盯制o-¨s，【。j“‘J’幽l-1Goodman单元图1-2Desai单兀两片接触面上各点在受力后产生f|匀j每封位移为舻1％』“-6)假定接触面上的法向应力和剪应力与法向相对位移和切向相对位移之间无交叉影响。则应力与相对位移的关系式为盯；D∞(I-?)式中。=心三l㈨s， 西南交通大学硕士研究生学位论文第8页t和矗分别为接触面七刀向和法向劲度系数(kN／m3)。取线性位移模式，可将接触面上任一点的位移用结点位移表示出来。r1“顶=I互f1‰2B(1—9)两式之差为切向相对位移m；=“顶一“麻a同理可写出法向相对位移∞。=”碰一”底。把它们用矩阵表示出来为：棚：I绌}-B占。(1-10)【曲J其中嚣=i：：：：苫三苫三I上式中口=j1一主，6=j1+·至由虚功原理可得单元平衡方程如下：Fc：f7．B7DBdx6‘：足cJ‘(1-11)‘号单元刚度矩阵为：k。：三62k。O2“(1．12)劲度系数的取值与应力变形状态有关。对‰，在接触面拉开时，取一个很小的数，否则取～个很大的数。‰可由直剪试验确定。Goodman单元能较好地模拟接触面上的滑移或张开，若(1-7)式中的D采用非线性的本构关系，还可以考虑接触面变形的非线性特性。然而，它也有两个大的缺点：一是单元无厚度，在受压时就会使两侧的二维单元相重叠，二是矗任取一个大值时，若法向相对位移有微小的误差，则会使盯。=岛l∞。有较大误差，因此导致计算的盯。有时是不合理的。有关节理单元本构以及求解方法的“Ⅳ、●11．／＼●11．／苎上x一￡一+，一21—2，，●●●●．．＼／，●●●●．＼+“、，ll／、、，●／x一￡x一￡+一批猕。称k，b{!蕈％。b执。缸o‰。％。咄执。珧。地。。k。“。‰缸。地。讹。 西南交通大学硕士研究生学位论文第9页详细论述，还可参考文献[31～331。2薄单元为了避免无厚度单元存在的两侧单元的重叠，许多学者主张用薄层单元。Zienkiewicz曾用等参单元作为接触单元，Ghabouss和Wilson也提出过另～种形式的有厚度单元，但计算较麻烦。这里主要介绍Desai提出的薄单元【341。如图1．2所示，接触单元长B，厚度t，在单元刚度矩阵的形成等方面与其它固体单元一样，可以取为四结点，八结点等参数单元。但在本构矩阵中，将法向与切向的分量分开来考虑。D：f仉孤1(1_13)D=Il(1—13)＼pnsDmj上式中JD。为剪切分量，眈。为法向分量，而Dm和五k为考虑耦合效应的分量。由于目前很难用试验测定法向和切向的耦合影响，丑k和田h可不予考虑，即取为0。法向分量D。。与接触面t厚度区域内材料的性质有关，也与两侧材料性质有关。可以将其表示为D一=丑(Dn)，+如(Dnn)。+也(Dm)。(1-14)上式中下标i，g和st分别表示接触区材料，岩土材料和结构材料，而九l、沁和b都是变化于0～1之间的系数。上式表示眈一是三种材料法向分量的某种加权平均值，九l、沁和b是权。一种简单的处理是令Xl=l、X2矗o=O。在动力问题中取九l-O．75、k=O．25、沁zO是合适的。剪切分量D。，是由接触面剪切模量G构成。其中G：堡f(1．15)a弛对一般材料来说，剪切模量G与弹性模量E和泊松比有关，G=j÷1。，⋯⋯但在接触面上，当t接近抗剪强度tl时，G会变得很小，而弹性模量不一定变小，压缩性不一定很大。令G≠志就便于考虑接触面的这种特性。DesaiZll+vJ在算例中取接触面材料的E=69000kP。，1)=O．3而G=I．38kP。，远小于由E，u算得的G值，这样就能使接触面两侧在受剪时的错动滑移得到模拟。由(1-15)式可见，单元厚度t的选择对G的数值有直接影响。当t取得太大，与B处于同一数量级时，接触单元与普通单元就没有什么差别了；当t取得太小了时，也会使相对剪切位移的计算产生误差。研究表明直取r=(o．01～o．1归(1-16) 西南交通大学硕士研究生学位论文第lo页以上的推导都是建立在四结点单元之上。当岩土或结构本分用八结点等参数单元时，就需要六结点接触面单元，更详细的推导请参照文献131。上述两种接触面模型在轴对称的计算中遇到了一些问题。问题主要出在环向劲度系数的确定上。段文峰等人推导了两种单元的轴对称格式，但同时他们也指出，推导过程中只是力学过程上的严密推导，至于环向应变的定义是否存在更合理的方式，还需要探讨。具体讨论请见文献[35]。上述的两种接触面模型都侧重于描述不同材料界面上的力学特性尤其是切向的摩擦特性，若将它们用于描述不同材料界面上的位移不连续现象尤其是大规模滑移和脱开时，常常难以得到合理的计算结果【l“。尽管不少学者对上述两种单元又作出了不少改进，但上述单元都是试图以连续的模型来模拟不连续的接触面，虽简单好用，但却没有触及问题实质，且这一缺点不是可以通过改进、完善所能弥补的。1．3．2接触问题对接触面的力学描述如前所述，本文将按接触问题计算分析桩一土相互作用，研究桩的竖向工作特性，以下将接触问题中对接触面力学特性的描述作简要介绍，既便于下文的应用，又可与上述计算方法进行比较。图1-3接触问题对接触面的力学描述如图1．3所示，两单元Q。及Q”之间为接触面，假设受力变形前后单元之间始终保持接触而不分离，则无论怎样受力和变形，在此过程中，接触面上任意一点处的作用力及位移应满足下列条件：(1)上、下面的作用力大小相等，方向相反儿=-p’。P。=一P’。(1-17)(2)Q及Qt不能相互侵入，即上、下面的法向位移满足“。=U’。(1-18)L 西南交通大学硕士研究生学位论文第11页的相容条件。(3)当接触面上的切向力小于其极限抵抗力(抗剪强度)时，沿切向不发生相对位移；切向力达到抵抗强度时，则沿切向发生相对位移。即当Jp，J0(塑性加载)d孑7：否7一孑“：{0(中性变载)(3-19)d仃=盯一盯=‘L中1，土义飘，、‘l<0(弹性卸载)因为d孑’只有在r次迭代完成并且应力调整后才能实现，而调整的时候又要用到d孑‘来判断塑性的加卸载，所以(3—19)中的d盯‘在程序中用d盯：来代替。设在第r一1次迭代后高斯点的应力为盯’‘，在第r次迭代后单元结点位移增量为彩’，则应力调整步骤为：(1)由单元结点位移增量计算应变增量幽‘=嚣彩7；(2)由应变增量计算弹性应力增量咖：=Dde’；(3)迭加各高斯点应力仃：=盯’1十打：；(4)计算各高斯点的有效应力；’1、；’和换算屈服应力盯；‘1(r)；(5)高斯点可能的发展状态有四种，按(3．17)、(3一18)判断高斯点所处状态；(6)根据高斯点的发展状态，计算调整因子R(0≤R≤1)，如果R>O，————————————————————__———————————————_—————————————————_————————————_————————_——一 西南交通大学硕士研究生学位论文第34页则对R如：进行消除，使按弹性计算的应力回到屈服面上。步；R的具体计算方法见图3-5。由d’=盯’1+出7及(B-13)、(3-14)得do’=d仃：一Dad,z上式中d五为一个正值的比例因子(又称塑性乘数)。以=——二；一aT尼I口：A+a’Da‘否则，转入第(7)(3-20)(3．21)图3-5调整因子R的计算通过(3．15)、(3．5)、(3．4)可以得到第r次迭代后的高斯点应力盯’。(7)对于弹性的高斯点，不用调整，应力为盯’=盯：=盯’1+d仃：。3．3接触问题娜蝴燃图3-6接触状态不愿幽在式(3．1)中，只所反映的是接触边界上的结点荷载，显然，它取决于该结点所处的接触状态。在计算中，结点状态所处的状态可能有三种：粘结、滑动、张开，如图3-6所示。另外，当结点张开后，在下一级荷载作用下有可能发生嵌入。“嵌入”表明接触体之间发生了相互侵入，这显然是不可能存在的，．一————————————————————————————————————————————————————————————一 西南交通大学硕士研究生学位论文第35页它只是由于不能预先判断物体受力变形后的状态，而使得点对的相对法向位移超过它们之间的距离时产生的，最终将通过迭代使其处于粘结或滑动状态。基桩在竖向荷载下的计算中一般不会发生“嵌入”，因此本文就不再讨论“嵌入”时的处理方式。表3—2所列是迭代过程中，由一种接触状态可能变为的另一种接触状态。表3-2计算过程中接触状态的变化初始状态粘结滑动张开嵌入可能变为粘结的状态张开滑动张开滑动粘结嵌入滑动表3—3给出了不同接触状态时，接触结点作用力的计算方法或计算公式。各种状态的判断准则及其由于接触状态变化而引起的接触不平衡荷载如表3．4所列。表中，接触结点的法向位移均以接触体力的外法线方向为正，此外，为避免混淆，接触结点的位移以““”而非前文的“d”表示。袭3-3接触结点作用力的计算初始粘结滑动张开状态接触只(P-。)由式(2-50)计结点只(，’。)、只(P，)直接算而得只=一P。。=0的作由式(2．50)计算而得只=sign(只)∥J只J只；一P’，=0用力(P’，=sign(P’，)ale。。J)袁3-4各种接触状态的判断准则及接触状态变化引起的不平衡荷载m-1次的接触状态m次的接触判断准则m次产生的接触不平衡荷载只<0粘(P’。<0)△P．=一△尸’。=01clcc+川￡1配=一AP’；=0结(1尸’；I(c+卢I尸。。I)滑只<0(P’。<0)叱=一AP。。=0只≥C+卢kl△e=只-sign(P．),ulP。1A以=一△只动(，。；≥c+∥I尸’。I)张只≥0△只=一AP．。=只(P’。≥0)廿。=一世：=P；开公式中C代表接触面的初始抵抗强度，即岩土和结构间粘聚力的等效结点 西南交通大学硕士研究生学位论文第36页值，如图3—7所示，其计算公式为C=2％，f，rN,出=2ZCq。p妻嘭(3-22)广。"3．4局部坐标系下的变形矩阵在本章第三节中，接触结点的接触不平衡荷载都是建立在各个结点的局部坐标系下的。其中，局部坐标系的：+轴平行于接触面外法线方向，：’轴顺时针旋转90度得到r。轴，局部坐标系轴，。方向与整体坐标r轴的夹角cl通过公式(3—33)求得，具体推导见3．6节。图3-8坐标变换图3-9接触结点的局部坐标系由于在接触单元平衡方程的推导中要用到坐标变换矩阵，下面简单介绍坐标变换矩阵计算。如图3-8所示，局部坐标转为整体坐标的转换关系为；㈧篡兰尝埘(3-23a)令上：l。?8口一81n口l，则上式变为：{s111口c．os口l／一二 西南交通大学硕士研究生学位论文第37页fr]f，P]1：}“1开}门-23”上式中三称为变换矩阵，同一结点的坐标变换、位移交换、结点力变换均用此矩阵。在整体坐标系下，单元结点的平衡方程为：F=BilDBi6(3-24)给定局部坐标下的≤，若利用上式计算，需要先将∥转化为整体坐标下的6．，转化公式为：6，=‘∥(3-25)然后将上式代如到(3．24)中，可求得整体坐标系下的结点力。由于接触判断时需要的是局部坐标系下的结点力，因此还需要将其转换到局部坐标系下。其转换公式为：曩+=口E=口骂TDB．厶d：(3—26a)令(F，)e：k’《巧⋯r，(∥)。=慨赴蕊⋯r，(F)。=k，￡：工3⋯jI，就可得到局部坐标系下单元的平衡方程：(F’)。=(丑’)1D嚣。(∥)。(3-26b)其中B’=B(L。)。(3-27)3．5接触结点处于计算域边界时的处理方法豳3．10接触结点上有位移约束图3-11接触结点上有集中力有时候接触结点会处于计算域的边界，此时接触结点有可能存在集中力、位移约束，接触结点所在单元边上可能有分布力等，此时对接触结点的处理与模型内部的接触结点处理方式稍有不同，当然与边界上的非接触结点也不相同。下面将分别进行讨论。 西南交通大学硕士研究生学位论文第38页(1)接触结点上有位移约束如图3—10所示，接触结点上，在与整体坐标Z轴交角为a的方向上位移受到约束。处理此类边界的方法是：在该接触结点上建立局部坐标系osn，其中，S轴垂直于约束方向，n轴平行与约束方向。在进行接触单元平衡方程计算时，用上述的局部坐标系的坐标转换矩阵代入，而不用用接触边求导的方法建立的坐标变换矩阵。余下的计算同模型内部接触结点。对接触结点i而言，可能的滑动方向只有一个，就是垂直于位移约束的方向，如果按接触边界计算的方向与f的可能滑动方向不垂直，那么，肯定处于粘结状态。因此在引入边界时对边界上接触结点先行判断可能滑动方向，看是否与i的计算滑动方向想符合，如果不符合，则可以将该接触结点变为普通结点，从而减少了计算量。当a等于0或90度的时候，该约束就成为整体坐标系下的竖向或水平约束。(2)接触结点上有集中力存在如图3．11所示，先将集中力F在整体坐标系下分解成为Fr、R，然后在以结点i上建立局部坐标系：S轴沿结点所在边的切线方向(假设与整体坐标r轴夹角为a)，n轴由s轴逆时针旋转90度得到(n轴的方向与结点i的法线方向平行)。最后将R、R转换到局部坐标系下即可。转换公式为：rF1rF1{～}：上『{。}(3—28)【FoJlF：J上式中，三为坐标变换矩阵。半》乎呐 西南交通大学硕士研究生学位论文第39页集中力最后都要分解成为与自由度对应的各个分量。在图3．12中，自由度方向中的实线表示初始接触结点的初始自由度，虚线表示可能增加的自由度。集中力F属于单元l，尸属于单元2，各自沿单元的自由度方向分解为日、咒、|P3、P4。对接触结点i而言‘，当处于粘结状态时，只有两个自由度(1、2)，此时计算时将乃、P4分别迭加到Fl、F2中；当处于滑动状态时，分离出自由度(4)，此时计算时将凡和PI迭加，而B和P4则独立计算。在增量计算中，有了上面的规则，各自由度之I剞的力的分配就不会造成混乱。如果没有告诉集中力属于那一个单元，显然在自由度的分离中新增的自由度对应的力的分配将出现混乱。(3)分布力和初应力对体力和分布力，先利用公式(2．25)和公式(2．26)将其转化为结点上的集中力，然后利用(2)中介绍的办法来处理。3．6接触结点局部坐标系的确定如前所述，接触结点的局部坐标系定义为：n方向与接触结点的外法线方向平行，s方向通过n顺时针旋转90。得到。因此，确定局部坐标系首先要确定接触结点的外法线方向。如图3一13所示，设过接触结点的外法线单位向量为n，，与s方向平行的单位向量为s，接触单元结点的排序旋转向量为珊。60根据右手法则确定，r始终顺着结点的输入方向，在图3．13中，若采用1⋯357．8的结点输入顺序，则得到的09和f如图中所示。平时我们所说的顺时针或逆时图3，13接触结点局部坐标系的确定针结点输入顺序，是指从垂直于单元平面向外的方向来观察得到的。由于观察着所处的位置不同，有可能造成不同的结点排序，因此本文并没有采用逆时针或j顿时针的概念来描述它，而采用旋转向量∞来描述。由几何关系知：；：；x吾(3—29a)由}=cosa’i+sinCr‘J，00=k，代入到上式得到：l玎=ICOS口’lO(3-29b)．．●，aS0C一口nS=_￡01．，眦oS 西南交通大学硕士研究生学位论文第40页盘1为f与整体坐标r轴夹角，由导数的几何意义可知：口’=Arctan(-竿)(3-30a)dr由于单元边上的整体坐标r和Z都是局部坐标{的函数。由复合函数求导法则：de：墨武．dr：竺d{0亭。a舌’将上式代入到式(3-30a)中得拈姗an(毫蝈嗉瑚)=Arctan(嘉筹(3-。ob)如果出／凿=3r／鹫=0，则代表单元边界上有结点重合，显然这种情况不可能发生：如果ar／a{=0，az／a掌>O，贝0t2"=万／2；如果静／af=0，拓／af<0，则12"=一7／"／2；将(3-30b)中计算的口’代入到(3．29b)中，则得到接触结点的外法线方向gt。然后将nj嘎时针旋转90。得到S。由于n和s处于同一个平面内，因此可以采用复平面向量旋转的公式来推导。记n=sinoe’一COS‘，，则s=n(-j、=一since’J—COSOf’(3-31)表示成向量的形式为：J=--COS0f’i—sin口’，=COS(石十a’)f十sin(；r+口’)．，(3-32)从上式可以得出S与整体坐标r轴的夹角d为：口=厅+口’(3．33a)上式是通过1—3．5，7．8的结点输入顺序得到的，如果采用1-8．7．5．3的结点输入顺序，则∞=一k，通过同样的推理可知：口=口-(3-33b)将(3．33)式计算的口即可代到坐标变换矩阵中进行运算。但要注意的是，实际结点的滑动方向有可能与t同向，也有可能与t反向，因此具体是什么方向，应结合结点切向滑动力(奠)的正负来判断。因此，滑动方向与整体坐标轴的夹角B可以表示为：口=Sign(P,")a(3-34)若结点同属于几个接触相关单元，则每个接触单元内计算出来的a有可能不同，因此此时应结合实际情况来进行判断。一般情况，可以采用取结点所在接触相关单元各自计算出来的夹角的平均值a即：口=(口。+口，)／2(3-35) 西南交通大学硕士研究生学位论文第41页3．7面向对象编程思想在本文方法中的应用3．7．1面向对象编程思想的特点常用的有限元程序都是用结构化的思想来设计的。结构化程序设计思想是一种自顶向下(top．down)的设计方法，即将大问题分解成更小的子问题，直到子问题一可以用一个函数或过程(子程序)来解决为至。这种传统的结构化编程模式也被称作是“面向过程”的编程方法。结构化的程序思想认为：程序=数据结构+算法。这种方法的优点是：程序具有良好的结构，易于设计、易于理解、易于调试，从而提高了程序设计和维护的效率。但由于采用了全局性的数据结构，随着有限元软件的增大，这种全局性的数据结构便大大降低了程序的灵活性，因而很难修改和扩展已有的代码以适应新的算法。这种不灵活性主要表现在：(1)算法和数据结构分离，修改和扩展代码需要开发者非常熟悉程序的整个数据结构；(2)子程序都是面向特定的问题而设计的，因此从其它程序获得代码的能力有限。(3)数据结构的微小变化将引起整个程序的波动：(4)各子程序之间的相互依赖关系不易确定，当子程序比较多时，对子程序管理非常复杂，而且子程序的功能也有可能产生耦合；(5)对有限元程序，变量众多，变量的命名和管理不便。图3．14面向过程与面向对象的比较面向对象编程(ObjectOrientationProgram，简称“OOP”)方法与结构化编程方法正好相反，它采用自底而上(down—top)的设计方式，提倡从构造相互独立的小模块(对象)做起，然后把小模块组合成更大的模块，直到最终建立所需要的应用程序为至。这些独立的模块不存在相互依赖性，因此是可以重甚b熏 西南交通大学硕士研究生学位论文第42页——————————————————_-_______-●，—____-____-—^_________-——-__--__-——___--___——-______—-—_____-一复使用的。为了使模块具有可组合性，必须提供模块之间相互作用的接口(消息)，而模块的具体实现则是可以隐藏的。只要模块的接口不变，修改模块内部的实现方式并不影晌其它的模块。软件的维护只是模块的更新和重新组装的过程。面向对象编程思想认为：程序=对象+对象+⋯，而对象=算法+数据结构。如图3—14所示。面向对象的模式与人们对客观世界的认识方式是一致的，该思想中的许多概念都来源于我们实际生活中，如继承、动作等，因此用它设计的程序可理解性要远远强于有结构化设计的程序。只是现在很多科技人员习惯了面向结构的设计方法，反而认为面向对象设计的程序不好理解，如果一开始就使用面向对象设计，那就不会出现这样的错觉了。面向对缘的思想通过面向对象的语言(Object．orientedlanguage)来实现。面向对象的语言必须具有如下四个特征：(1)抽象(Abstraction)抽象是指用程序语言来模拟问题的现实特征。抽象的过程是从问题中抽取最令人感兴趣的特点。由于计算机程序就是用来模拟现实世界，所以所有的计算机语言都具有不同程度的抽象能力。利用Fortan、C等面向过程的语言抽象出来的对象数据和过程是分离的，而利用面向对象的语言的可以将数据和过程统一为一个整体。(2)封装(Encapsulation)封装是把相关的数据结构和算法封装成一个独立的模块，即对象。因为对象是独立的模块，所以它不依赖于其它对象。封装实现了信息隐藏，对对象内部信息的访问只能够通过接口来进行。封装提供了程序保护对象内部数据的能力，可以减少直接访问内部数据而造成的错误。在OOP中，一般都有几种控制方法，如Private、Protected、Public等，各种控制方法提供了不同级别的保护。(3)继承(Inheritance)继承刻画了对象的一般性和特殊性，它是指一个对象可以继承祖先对蒙的特点，同时有可以拥有自己的特点。继承性可以大大加强代码的可重用性。(4)多态(Polymorphism)多态是指同一类的对象对相同命令的不同反应。多态使程序的灵活性大大增强。在OOP中，多态是通过对方法的“动态联编”(DynamicBinding)来实现的，它是在程序运行阶段根据引用所指的实际对象来确定所调用的方法。相对应，在设计阶段就确定了对对象方法的调用称之为“静态联编”(StaticBinding)。只有完全具有上述四个特征的语言才是面向对象的语言。面向过程的语言(如Fortran、c语言等)也可以用来设计对象，但却没有对象的继承、封装、 西南交通大学硕士研究生学位论文第43页多念等主要优点。通常，我们还会遇到基于对象的语言(Object．basedlanguage)，如Microsoft公司的VisualBasic，在这种语言中，我们只能使用系统已经定义好的对象，而不能自己创建对象，而且该语言中的对象也不支持继承。本论文采用是Inprise(原Borland)公司的Delphi平台，使用的是ObjectPasacal语言，浚语言是一个真正的面向对象的语言，开发一个程序的过程就象创建一个对象的过程。OOP中，还有一个重要概念：类。类是用户定义的数据类型，通常用来描述一些非常相似的对象所具有的共同特征和行为。对象是类的实例，也就是通过类这个数据类型定义的变量。显然，先有类后有对象，所以更严格的讲，上述面向对象语言的特点是对类而不是对对象。对象是实体，在程序运行中，对象要为它内部的数据和方法占用一些内存，而类则作为一种数据类型，不会在运行期间占用内存。每一个类都是从另一个类继承下来的，被继承的类称为“父类”。在自定义类时，父类名并不是必需的，当在类的声明中没有指定父类时，Delphi会默认这个类的父类是TObject类。TObject类是所有类的父类。每一个类中，必须有构造函数(Constructor)和析构函数(destructor)。类利用构造函数来创建对象，利用析构函数来释放对象所占用的内存。面向对象编程的第一步就是对各个类的抽象过程，也是最关键的一步。与面向过程编程相比，它先难后易，而面向过程则是先易后难，完成某一个特定任务的工作量差不多，但对于以后的维护和更新，显然前者要远远强于后者。3．7．2面向对象编程思想在本文方法中的应用“自由度法”的特点就是所有的分析都是基于自由度基础之上，而不是基于结点之上。结合“自由度法”的特点，在弹塑性有限元的程序基础上，作者设计了面向对象的接触有限元程序。在“自由度法”中，和接触直接相关的类有：结点类(TNode)、单元类(TElement)和自由度类(TFreedom)。由于所有关于接触的计算都是发生在接触结点上的，只需对接触结点进行处理即可，处理后的结果反馈给单元的过程与普通结点并没有什么本质的区别，因此对于单元类，作者只进行了简单的扩充就将起应用要接触相关单元。不过，还里需要在单元类里定义一个接触单元指示变量(Boolean类型)，该值为True时表明该单元为接触相关单元，为False时表明为普通单元。下面将详细介绍～下接触结点类和自由度类的设计和实现。接触结点类的定义如下：TContactNode=Class(TNode) 西南交通大学硕士研究生学位论文第44页lnheriated；，，继承父类Freedom：TContactFreedom；／／重新定义结点的自由度类RelateElementNum：Array[1—2】ofInteger；／／接触相荚单元编号CEdge：Array[1．．2】ofInteger；／／接触相关单元中用作接触面的边的编号C，V：Array[1．．2】ofSingle；∥对应接触边的力学参数NodeFreedom：Atray[1．．4】ofInteger；／／重新定义接触结点的自由度ContactMode：Integer；／／接触状态，1为粘接，2为滑动，3为张开ExtraLoad：Array[1．．4]ofSingle；／／接触结点的不平衡荷载Alpha：Single；臌触结点局部坐标系与整体坐标系的夹角“Resistance：sin91e；臌触结点的初始抵抗强度FunctionMakeAlpha0：Single；／／形成局部坐标系s轴与整体坐标系r轴的交角ⅡProcedureMakeResistance0；Ⅳ形成接触结点的初始抵抗强度ProcedureContactModeCheckO；，，分析当前接触结点的接触状态ProcedureMakeExtraLoadO；fl形成接触不平衡荷载；End：自由度类的定义如下；TContactFreedom=Class(TFreedom)Inheriated；N形成接触相关单元的各结点自由度ProcedureElementFreedom(ElementNum：Integer)；，，供单元类调用ProdedureNodeFreedom(Status：Integer)；End；在结点类和自由度类中，将所有的方法全部定义为Virtual方法。对于结点类，也需要和单元类～样，增加一个接触结点指示变量ContactNode(Boolean型)，如果为True，则表明该结点为接触结点，否则为一般结点。其它方法和属性与普通的有限元没有什么大的区别，考虑到目前已经有很多文献对结点类和自由度类有国大量介绍，本文就不再此详细展开了，具体可参考文献[52～58]。当接触结点类和接触自由度分别继承了结点类和自由度类之后，还可以重新定义父类中已有的方法，如接触自由度类中的NodeFreedom(ElementNum)方法。这样做有一个好处，就是程序可以根据实际变量的类型，动态的调用对象的方法，实际上也就是面向对象的多态性特征。其具体实现过程如下：vatNode：Array[1．MaxNode】ofTNode；，，声明变茸forj：=1toMaxNodedobegin砒处为调用语句 亘重至望盔堂巫主受塞生兰焦迨室篁箜耍ifNode[i]，ContactNodethenNode[i]：=TContaetNode．CreateelseNode[i]：=TNode．Create；end；在面向对象编程中，子类可以先定义为父类类型，然后在程序中根据需要再动态地创建需要的子类。由于子类继承了父类的方法，上述语句才可以进行，如果是任意两个类，程序则会报错。在for循环中，先对结点进行判断，看是否是接触结点，如果是接触结点类，则根据接触结点类创建Node[i]，以后对Node[i1的所有调用都是采用的是接触结点类中的方法，否则则是对普通结点类的调用。这样，我们就不必重写程序分别进行处理，大大提高了程序的效率和可读性。图3-15接触结点类与相关类的消息传递图3-15给出了接触结点类与其它相关类的消息通讯图。从大的方面来讲，整个有限元程序将作为一个对象类，并定义为Domain。接触结点类通过Domain类来初始化，得到结点坐标、荷载等信息。TContactFreedom类通过接触结点的接触状态，分配给接触结点自由度，一次接触迭代完成后，将本次的迭代后的接触状态再反馈给接触自由度类。接触相关单元计算单刚时需要接触结点信息，在分析完成后，接触相关单元又将接触信息反馈给接触结点。接触结点通过接触相关单元反馈的信息，通过调用方法MakeExtraLoad0计算接触不平衡荷载，然后提供给荷载类。在所有的消息传递过程中，各对象内部的数据可以根据需要进行不同级别的封装。对只有在对象内部使用的，就将其定义为Private型，如果需要和外部交换。则将其放到接口部分，并将其定义为Public型。这样以来，各个对象之间的相互关系也就很明确，程序的可读性和可扩充性大大增强。 西南交通大学硕士研究生学位论文第46页第4章接触问题算法在基桩计算中的应用——算例在前述的理论研究和程序编制工作完成后，本章将给出几个算例。这些计算的主要目的有两个：一是验证本文计算方法的正确性及程序的可靠性。二是试验用方法求解基桩问题的可行性。对第一个目的，作者针对同一个问题，通过对弹性理论解和有限元解的比较来判断程序的正确性，同时利用计算结果指出了本文方法的特点；对第二个目的，则通过与现场测试、观测结果对比，对计算结果的合理性进行定性的判断。4．1有限元法与弹性理论法结果的比较4．1．1基桩计算的弹性理论法””把桩分为若于个单元，则在J单元作用环形分布力研后，在j段引起的变化(“．)，，可以写为(图4一1)(％)，2詈和，(4-1)其中臣为桩周土的弹形摸量，d为桩身直径，五·为T厂L]单元J的摩擦力使单元j处产生竖向位移的系数。1．1I脒罐黼靴6蒙爹d≯龋翟：，r曰协∽m=瓦加a(4．2’fI一口f乃把桩与土接触面上所有影响力累加起来，可得l1％=苦跏≯d“烈川⋯㈣I口相应地，桩底变位可以写为：Il铲鲁跏+扣。㈣一，一墅加上述式中f。、，扩，矿k代表J点和桩底分布力作用下在桩191i|j点和桩底面引起的无量纲变位，可通图4-1Mindlin解用于单桩分析综合(4一1)～(4-4)，桩所有单元的土位移可用矩阵表示为：H：昙，，f(4—5)H2i1，f”～’ 西南交通大学硕士研究生学位论文第47页式中“一n+l维土位移矢量；r～桩侧剪应力和桩端应力矢量，n+l维；，。一土位移系数矩阵(n+l阶)，详见文献[60]；另一方面，桩受力轴向力P作用，若只计桩轴向力的压缩作用，忽略径向力的影响，则桩的平衡微分方程可表示为：一0w：竺—L(4—6)瑟dE。R』式中W为桩身位移，t为桩侧剪应力，露为桩的弹性模量，尼为桩截面面积，对于实心桩，届=l。上式可以写成有限差分的形式，依次应用于计算点i=l～n。写成矩阵的形式为：f=i≥EPJRd，pH，+】，(4—7)4占2”“9式中f—n+l维剪切应力矢量；炒一n+l维桩位移矢量；，。一桩作用矩阵(n+l阶)，详见文献[60]；y—j旦oo⋯ooll舭j由位移协调条件，界面不发生滑移，界面上各点的土位移等于桩位移，即：／．L=w(4—8)将(4—5)和(4—7)代入上式可得：r=[I-』熹叩^纠一y(¨)4(Z／们29“P”式中j—n+l阶的单位矩阵；^一桩的刚度系数，由下式给出：足=≥R』(4-10)L5从式(4—9)求出桩侧阻力和桩端阻力的分布，然后通过(4—7)得到桩身位移。4．1．2有限元饵与弹性解的比较’为了充分反映本方法的特点，作者采用了三种计算方案。其中前两种方案按弹性接触有限元计算，第三种方案按弹性理论计算。各方案的桩、土参数不变，只改变接触面的参数，最大加载量均为10000KN。桩的弹性模量为3×207Kpa，治松比为0．2，桩长L=25m，桩径d=lm，土取均匀土层，其变形模量 西南交通大学硕士研究生学位论文第48页为3X105Kpa，泊松比为0．2。接触面的参数见表4．1。各计算方案的计算桩顶沉降见表4—2，相对应的荷载一沉降曲线见图4．2。表4-1各计算方案的接触面参数计算方案粘聚力(kPa)摩擦系数备注很大l加载始终结点不滑动_l60O0．90．5参数为虚拟的从表4—2和图4．2可以看出：(1)有限元与弹性理论法求解的桩顶沉降吻合较好，相对误差仅为0．81％。(2)而且在荷载小于某一个数值(方案二为5000KN，方案三为3000KN)的时候，用接触有限元计算的结果与用连续体的有限元计算结果相同，这是由于在荷载比较小的时候，桩周土和桩可以共同变形。但当荷载较大时，利用接触有限元计算的结果呈明显的非线形，随着荷载的增大，用接触有限元计算的结果与连续体有限元计算的结果相差越大，这是由于用连续体的有限元无法考虑桩土的相对滑移而造成沉降偏小的缘故。图4．3则描述了用按接触有限元计算的桩身滑动过程。当荷载小于3000KN的时候，桩和土共同变形，没有发生滑动。当荷载逐渐的增大，桩身滑动范围也逐渐向桩底扩展，直到桩身全部滑动，这个滑动的过程也是荷载向下传递的过程。、＼＼方案荷载、、＼＼(KN)＼(接连续体)(弹性接触有限元)(弹性理论法)10000．492000O．980．980．9930001．471．4840001．962．651．98500021453．682，4760002，944．852．9670003．436．173．4680003．927．653．9590004．419．364．45100004．913．054．94 西南交通大学硕士研究生学位论文第49页￡E一鞋蛙2000荷载(KN)4000∞OO800010000图4-2各方案计算所得的荷载-沉降曲线P=4000KNP=6000KNP=8000KNP=9000KNP=10000KN图4-3方案二计算所得的桩身滑动范围发展过程4．2西昌某工地试桩的有限元计算及分析4．2．1概况在西昌至攀枝花的高速公路的几座大桥工程中，根据桥梁的荷载情况和桥位处的地层结构，拟采用拟采用钻孔灌注桩作为墩台基础。大桥基础所在地层分布中均存在着厚度较大的第四系．冰期昔格达组泥岩，该岩层为河湖相静水沉积，呈半成岩状态，胶结不良，岩石强度极低，基本上由胶结物强度所决定，O2468024 西南交通大学硕士研究生学位论文第50页属极软岩类。现行各规范对极软岩中的桩基尚无明确完整的计算方法，为保证桥梁基础设计的科学合理性，也为类似地层中桩基设计与施工提供依据和积累资料，确定通过现场静载试验等方法研究桩的竖向承载特性。试验场地选在德昌县境内蒲坝场的安宁河左岸河漫滩上。根据提供的地质报告，场地处的地质概况可描述如下：①低液限粉土，层深0H1．2m，棕褐～黄灰，结构不均，O．50ra以上为耕植土，局部夹小块石和卵石，软塑状：②漂石质土，层深12～7．0m，黄灰色，石质成分以花岗岩、石英岩为主，部分漂、卵、砾石强风化，次圆～圆状，砂及粉粒充填，结构不均，局部卵、砾石富集，松散～中密，稍湿～饱和，透水性一般；③泥质粉砂岩，层深7．0～26．73m，褐灰色，矿物成分以石英、长石和粘土类矿物为主，粉粒结构，薄厚层状构造，泥质胶结，胶结极差，软硬不均，层理近水平发育，结构不均，局部砂、泥质分别富集，具遇水泥化、软化和脱水开裂特征。根据室内土工试验的结果，岩石的单轴抗压强度小于1MPa，软化系数在0．4～O．6之间。试桩共有八根，其中两根为人工挖孔桩，其余六根为钻孔灌注桩。由于l#(人工挖孔，空底)和3拌(冲击成孔，实底)具有代表性，本论文选取为他们计算研究对象。为了减小覆盖层对桩侧摩阻力的影响。在试桩穿越覆盖层的区域内设置了铁皮隔离，而且对覆盖层土进行了开挖回填。这里的覆盖层是指地质报告中的④和②土层。试桩桩身混凝土采用C30，纵向受力主筋采用II级螺纹钢筋。根据现场施工记录，计算所用试桩的施工情况如表4．3所示。表4．3计算试桩施工情况试地面标高桩底标高桩入土深度基岩面标高覆盖层厚度入岩实际桩成孔深度编方式H1H3H=fHl-H3)H2h=(H1-H∞桩径号(m)(H．h)(m)人工11216．7661195．86620．9001213．7663．00017．9001．25挖孔冲击●31216．8361196．02020．8161212．0004．83615．9800．85成孔4．2．2计算模型及计算参数桩的计算模型如图4．4所示。该问题显然可作为轴对称问题进行计算。此外，桩一土(岩)之问考虑为接触面，当桩表面一点处的竖向剪应力超过其抵抗强度(桩侧极限摩阻力)时，在该点处，桩一土(岩)处之间将产生相对滑 西南交通大学硕士研究生学位论文第51页动a其相应的判断准则也采用Mohr--Coulomb准则。图4，5为相应的计算网格图，图中的L为桩的长度，d为桩的直径。图4-4桩的计算模烈(轴对称)士中的接触相关单元Io阳【早JL“-Jn一己0d～ 西南交通大学硕士研究生学位论文第52页计算中采用的参数列于表4-4和表4．5，其中，昔格达岩的变形模量依现场压板试验的结果确定，强度指标则根据室内试验结果确定。因没有进行相应的试验，桩一土(岩)接触面上的强度指标按经验选取，即：桩一土(岩)沿接触面发生滑移，其相应的强度指标与岩、土的抗剪强度指标接近，对冲孔桩，考虑到施工时的护壁泥浆等因素的影响，其相应的指标应予适当降低。此外，无论是人工挖孔桩还是冲孔桩，发生滑移后，摩擦系数在原有的基础上降低，而粘聚力则为O。表4—4材料参数表、＼指标变形模量泊松比内摩擦角粘聚力名趴(MPa)(o)(kPa)＼覆盖屡(回填)5．00-35强风化昔格达岩150O．33570弱风化昔格达岩300O．34550桩(钢筋混凝土)30000O．2表4．5桩一土(岩)接触面上的强度指标1#桩3#桩粘聚力(kPa)摩擦系数粘聚力(kPa)摩擦系数滑移前滑移后滑移前滑移后滑移前滑移后滑移前滑移后覆盖层00．200．2(回填)强风化100OO．85O．7560O0．70．45昔格达岩弱风化85O1．000．840O0．9O．5营格达岩4．2．31#桩的计算结果及分析1拌桩为人工挖孔桩，且为空底。桩的直径取为1．25m，覆盖层(经铁皮隔离处理后)厚度4．5m。计算时，分8级将荷载加至15000kN，之后，分4次卸至0。表4-6及图4-6为由计算所得的荷载—沉降关系，在荷载小于14000kN，桩的沉降较小，计算结果与现场试验结果(图4．7)接近。当荷载增至15000kN时，沉降急剧增加而无法收敛，表明桩侧摩阻力达到极限状态，整个桩与岩(土)层之间产生相对滑移，而由于空底桩桩端无法提供端阻，故此时桩的沉降已无法稳定。同样，出图4-6还可看出，卸载后桩的变形基本完全恢复，而卸载完毕后的残存沉降量主要来源于桩一岩(土)之间的相对滑移，这部分沉降在卸载后 西南交通大学硕士研究生学位论文第53页几乎全未恢复。由计算结果知，1#桩的极限承载力应在14000“15000kN之间，此桩试验时未达到极限状态，但在可对比荷载段(图4—7)，计算与试验结果是很接近的。图4—8是计算轴向力沿桩的分布图，同样，在覆盖层范围内，侧摩阻力很小，轴向力的分布显得较为平缓，进入岩层后，轴向力沿桩长衰减速度很快，在桩底基本为0(计算时有一定误差)。与实测的轴力分布曲线(图4．9)比较可知，两者在昔格达泥岩中的变化规律相似，其轴力在泥岩中衰减速度都较快。表4-61#桩由计算所得的荷载一沉降关系荷载桩顶端沉降桩底端下沉量桩的变形量(kN)(mm)(ram)(mm)020001．310．580．7340002．631．161．4760003．941．742．280005．262．322．94100006．572，93．67120008．023．694．33140009．664．734．931500065．1260．234．891100065．0760．874．2700062．0359．282．75300059．3858．1l1．27O57．3857．23O．15 西南交通大学硕士研究生学位论文第54页荷载(KN)。o．200040006000800010000120001400016000ulO．一总沉降—、1IE+桩底沉降＼＼￡，；—-．一桩的变形＼＼一，邀’墨55’I龟≮≤芝=一60—RE一⋯一70一图4巧1#桩计算所得的荷载沉降曲线荷载(KN)o20004∞06000800010000126001400016000u’之三+实测桩顶沉降＼5·10。一E邑15：世蜉55。I‘～～～～760一～一70一国4-71#桩计算所得的荷载沉降曲线与实测曲线的比较 西南交通大学硕士研究生学位论文第55页言’1。一冀啪．20一250‘5童m嫠舶桩身轴力(KN)O2000400060008000100001200014000图4-8l#桩由计算所得的桩身轴向力分布桩身轴力(kN)02000400060008000100001200044000图4-9l#桩实测的桩身轴向力分布4．2．43≠}桩的计算结果及分析3号桩为冲孔桩，桩的直径取为0．85m，覆盖层厚度为6m．计算时，分7级将荷载加至14000kN，之后，分3次卸至0。 西南交通大学硕士研究生学位论文第56页表4—7及圈4-10为计算所得的荷载一沉降关系，由图可知，在荷载较小时，P-s曲线呈线形，之后随着荷载的增加，沉降速率加快，曲线里非线形。在荷载较小时，桩身的变形在总沉降中占据较大比例，但随着荷载的增大，桩端的岩土中的塑性区域逐渐增大，致使桩端沉降逐渐在总沉降中占据的比例越来越大，最后超过桩的变形，卸载完毕后的残存沉降量主要来源于桩底岩层的塑性变形。图4．11对计算荷载一沉降曲线与实际测试曲线进行了对比。由图可知，在荷载小于6000KN的时候，计算曲线和实测曲线吻合的较好，但当荷载继续增大时，实测曲线的斜率增加较计算曲线快，原因可能是在桩底部的岩(土)出现了剪切破坏面，而程序计算时只能模拟塑性破坏，而无法模拟因破坏产生的滑移。图4—12是计算所得的轴向力沿桩的分布图，荷载较小时，轴向力沿桩长衰减速度较快，荷载逐渐增大时，桩一岩之间发生相对滑移的范围出上向下逐渐发展，荷载也逐渐向桩底转移，故轴向力的分布显得较为平缓。与桩的实测桩身轴力(图4—13)比较可知，两者规律相似。在覆盖层中，计算桩身轴力变化较小；进入昔格达泥岩后，桩身轴力随深度增加而开始较大幅度的衰减；到桩端附近，由于桩断土的影响，桩断附近的桩身轴力则表现为：越靠近桩端，轴力衰减的越快。。图4．14为计算所得的桩侧摩阻力分布图。从图中可以看出，覆盖层提供的摩阻力很小，在桩荷载向下部传递的过程中基本保持不变。在荷载刚传递到强风化的泥岩时，摩阻力分布特征为中间大、两端小，随着荷载的迸～步增加，强风化的泥岩提供的摩阻力逐步达到极限，其后则基本保持不变。在桩端附近，摩阻力增加较快，尤其是当荷载较大时更为明显。桩身的全部摩阻力随荷载的增加而增大，桩一岩之间发生相对滑移的范围也在逐渐加大，由于6000kN左右时，整个桩与岩层的接触面均已发生相对滑移，所以其侧摩阻力也相应地达到极限状态，此后已基本不随荷载的增加而增大，剩余部分的荷载则由桩底岩层承担。图4．15给出了桩的总侧阻力和端阻力随外荷载的变化曲线。从中可以看出，端阻力在6000KN以前发展缓慢，而6000KN以后发展速度较快，且呈非线形增长的态势。总的侧阻力则相反，6000KN以前，发展迅速，以后基本保持不变。由图4—16可知，加载至6000kN时，滑动范围已扩展到全部桩身，沉降曲线开始出现明显弯曲，此后沉降以较大的速率增加，至12000kN时沉降量达到27．1mm，相当于现场试验8000kN时的沉降量。荷载继续加至14000kN时，沉降增加至34．5mm，此时，桩底岩层己出现较大的塑性区(见图4．19)，计算迭 西南交通大学硕士研究生学位论文第5"7页代时的收敛已经相当困难。全部荷载卸完后，残存的沉降量为8．7ram。图4—17和4—18给出了在14000KN时的桩周土位移等值线。对X方向位移而言，最大值发生在桩底的土层内，而且从等值线的密度也可以看出桩端附近一定范围内的位移都较大。而Y方向的等值线则表现为：桩端附近的桩周土位移和桩底附近土层的位移相等，在接触相关单元取内，位移的等值线呈漏斗型，这是由于桩向下滑移的过程中，带动了土的变形。表4—73#桩由计算所得的荷载一沉降关系荷载桩顶端沉降檄厩端r讥量桩的变形量(kN)(him)(mm)O0O20002．40．81．640004．91．73．2600010．94．66．3800015．87．18．7100002110111200027．113．813．31400034．518．815．71000029．817．812600023．215．37．908．780．7 西南交通大学硕士研究生学位论文籀58页E一遨基E一世娉荷载(KN)0200040006000800010000120001400016000幽4．103#桩由计算所得的荷载沉降曲线荷载(KN)o200040006000800010000120001400016000图4-113#桩计算所得的荷载沉降曲线与实测曲线的比较o∞帖加筠∞拍o伯任{己为如弘 西南交通大学硕士研究生学位论文第59页～0桩身轴力(KN)O20004000600080001000012000440000·5^E．100图4—123#桩由计算所得的桩身轴向力分布一+一12000KN—--Ⅸ——14000KN桩身轴力(kN)400060008000圈4—133#桩实测桩身轴向力分布巧伯帖一S谜避 西南交通大学硕士研究生学位论文第60页ZYR舞010桩侧摩阻力(KPa)0100200300400500600700800圈4．143#桩由计算所得的桩侧摩阻力分布02000400060008000100001200014000荷载(KN)图4．153#桩的端阻力和侧阻力随外荷载的变化一Ⅲ)越囊舢蛐㈣姗撇。 西南交通大学硕士研究生学位论文第61页圈4-163#桩的桩身滑动过程示意图图4一173#桩周土x方向位移等值线幽4—183#桩周土Y方向位移等值线 西南交通大学硕士研究生学位论文第62页P=800咄NP=1000C略Ip=-12000KHP=140131]咧"图4-193#桩由计算所得的桩周岩土塑性区发展过程4．3本章小节由上述计算结果和分析，我们得出以下几点结论：(1)出计算所得的荷载沉降曲线和实测的曲线符合的较好，空底桩在大于14000KN后由于滑动范围已经扩展至全部桩身，桩侧摩阻力全部发挥，桩达到极限承载力；而实底桩在荷载沉降曲线上则表现为：小于6000KN时，两者吻合地很好，大于6000KN后，荷载传至桩底，实测曲线下降速率大于计算曲线，可能是由于桩底出现了剪切破坏面，而程序只能模拟塑性破坏，而无法模拟剪切破坏发生的滑移。(2)实底桩和空底桩的荷载一沉降曲线有明显不同。荷载较小时，实底桩的沉降随荷载基本呈线性增加：继续增大时，则呈非线性增大，由前文的分析知，其分界点对应于桩侧摩阻力达到极限时。对空底桩，桩侧摩阻力达到极限就意味着承载力已达到极限，桩的沉降将急剧增大。(3)对空底桩，在加荷过程中，桩身变形在桩的总沉降中始终占据较大比例：而实底桩则表现为：加载初期，荷载较小时，桩身变形大于桩端沉降；荷载较大时，桩身变形小于桩端沉降，这是由于桩端岩土的塑性区逐渐扩大而造成的。(4)由实底桩侧摩阻力分布及荷载一侧摩阻力(平均值)曲线可以看出，当荷载增加时，桩侧摩阻力由上至下，逐渐达到极限值，直至整个桩的摩阻力均处于极限状态。此后，若荷载继续增大，将“剩余荷载”转由桩底岩层承担，相应地，桩底应力迅速增加。(5)从实底桩的竖向位移等值线可以看出，桩端附近的桩周土位移和桩底附近土层的位移相等，在接触相关单元内，位移的等值线里漏斗型，这是由于 西南交通大学硕士研究生学位论文第63页桩向下滑移的过程中，带动了土的变形。(6)由计算结果知，l#桩的极限承载力应在14000’15000kN之间，此桩试验时未达到极限状态，但在可对比荷载段，计算与试验结果是很接近的；3#桩的极限承载力约为14000kN，高于实测结果的10000kN。(7)本计算的结果与试验结果在定性上与定量上均吻合得较好，说明本文所采用的计算方法是合理、可靠的，可作为研究桩受力及承载特性等问题的一个有效工具。 堕宣窒塑大_芋硕士研究生学位论文第64页结论本论文的目的是建立一种轴对称接触问题的有限元计算方法，并将其用于基桩竖向工作特性的计算分析。为此作者完成了以下工作：(1)利用虚功原理推导了接触问题的平衡方程，为有限元的计算公式奠定了基础。所得到的计算公式比较直观，计算也比较容易实现。(2)建立了轴对称接触问题等参元有限元计算公式。(3)分析了“自由度法”的特点，并将其应用到桩土相互作用的求解当中。(4)研制了轴对称接触有限元的计算程序。该程序综合考虑了材料的弹塑性、不连续面等因素。另外，与接触有关的几个问题，如接触结点的滑动方向，接触面的初始抵抗强度，接触结点位于模型边界时的处理方法等都在程序中得到了体·现。最后，利用面向对象编程的思想，建立起轴对称接触问题有限元的对象类框架，并对程序中可能遇到的问题进行了探讨。(5)利用编制的程序分析了西昌某工地的两根试桩，验证了方法的可行性。通过计算得到了试桩的荷载沉降曲线、轴力分布、桩侧摩阻力分布、桩周土的位移等值线、桩周土的塑性区等结果，并与现场试验结果进行了对比，两者吻合较好，说明本文所采用的计算方法是合理、可靠的。通过上述工作，作者得出以下几点结论：(1)“自由度法”是建立在系统自由度上的平衡方程，每一个自出度与结点有关，但并不完全依赖与结点。针对于接触问题，由于接触状态的改变，系统的自由度也随之改变，“自由度法”通过对增加与回收来实现的，处理接触问题时比传统基于结点的分析方法更合理，也更方便。(2)将桩土之间相互作用作为接触问题来处理，可以更好地反映出桩土相互作用的力学特性。通过对西昌某工地的两根试桩的分析结果也显示，得到的计算结果是合理的。(3)对桩基而言，桩土之间的接触面对桩的工作特性是至关重要的。考虑接触面的弹塑性接触有限元计算结果比弹性理论解跟趋予实际，尤其是在荷载较大的时候。(4)面向对象编程思想在数据结构复杂，更新快的程序中得到了广泛的应用。弹塑性接触有限元程序正具有以上特点，因此很适合在该类程序中引入面向对象 西南交通大学硕士研究生学位论文第65页编程的思想。桩土之涮的不连续面是影响基桩工作特性的重要因素，将该面作为接触问题来研究，是一种有效而合理的途径。桩基工作特性不但依赖于地层情况．还强烈依赖于施工方法，因此没有一种计算方法能适应各类桩基础，本文方法亦是如此。总的看来，作者认为还有以下几个问题需要解决：(1)针对方法本身而言，由于接触状态的变化会对系统的总刚维数发生变化，增大了计算的工作量，因此怎样提高求解效率是亟需解决的一个问题。(2)对打入桩而言，结构的失效可能并不发生在桩土界面上，而是发生整体或局部剪切破坏。这样，破坏面的位置事先并不知道，破坏面是从无到有，逐渐发展的，本文的方法目前还不能对打入过程及桩周岩土体的变化进行模拟。(3)本文方法所需要的计算参数虽然较易获得，但计算参数的精度对结算结果影响较大，尤其是接触面的摩擦系数。因此，怎样准确的确定接触面的各项参数也是需要研究的一个问题。 西南交通大学硕士研究生学位论文第66页致谢本论文是在导师毛坚强的悉心指导和耐心帮助下完成的，从论文的选题直至完成，都凝聚着导师的心血，而每次在程序编写最困难的时候，都是导师的耐心启发和循循善诱的讲解才使我度过了难关。另外，导师那幽默豁达的生活态度也令我终身难忘。值此论文付梓之时，特向我的导师致以深深的敬意和诚挚的感谢!本论文的算例取自作者在西昌某工地参加的一个试验。在试验过程中，教研室吴兴序老师给了我很多的帮助和支持，在此深表感谢。另外，作者还要感谢于志强老师给了我很多参加科研实践的机会，让我锻炼了自己的工作能力。感谢教育过我的各位老师，是您们的教导才使我走入了知识的殿堂。同时，还要感谢论文参考文献的作者，没有你们的辛苦劳动，我也不可能这么快就能完成论文。在论文的编写过程中，还得到了师兄徐骏，许敬以及学友刘志勇，彭世江的关心与帮助，在此一并表示感谢。作者还要感谢中国建筑西南勘察研究院总工程师康景文，技术中心及人力资源部的各位领导对我的关心和帮助。最后，深深感谢所有支持和帮助过我的老师，师兄弟，亲朋以及好友们!王新2003年10月 亘塑至遒查兰塑主壁窒生堂焦鲨塞签!!蔓参考文献1《桩基上科手册》编写委员会．桩基．_J：程手册．北京：中国建筑工业出版社．19952毛坚强．接触问题的一种有限元计算方法及其在岩土工程中的应用．西南交通大学博士研究生学位论文．20023钱家欢，殷宗泽．土一I：数值分析．北京：中国铁道出版社．19914BatheKJ，WilsonEL．NumericalMethodsinFiniteElementAnalysis．Prentice-Bin，Inc．，19765BatheKJ．FiniteElementProceduresinEngineeringAnalysisPrentice—Hall，Ine．，19826CloughRW．TheFiniteElementinPlaneStressAnalysis．AircraftEng．19547wilsonELStructuralAnalysisofAxisymmetricSolids．AIAAJ．1965：V01．3(12)：2269～22748CloughRWandWoodwardRJ．AnalysisofEmbankmentStressesandDeformations．J．SojlMech．Founs．Div．，ASCE，V01．93，SM4，19679王勖成，邵敏．有限单元法基本原理和数值方法(第2版)．北京：清华大学出版社．199710李开泰，黄艾香，黄庆怀．有限元方法及其应用(1)一方法构造和数学基础．西安：西安交通犬学出版社．198411苏煜城，吴启光．偏微分方程数值解法．北京：气象出版社．198912龚晓南．土塑性力学(第二版)．杭州：浙江大学出版社．200113何君毅，林祥都．工程缩构非线形问题的数值解法．北京：国防工业出版社．199414JohnsonKL．ContactMechanics．CambridgeUniversityPress，198515孔样安，江晓禹，金学松．固体接触力学．北京：中国铁道出版社．199916张丙印．师瑞锋，王刚．高面板堆石坝面板脱空问题的接触力学分析．岩土二l程学报．2003：V01．25(3)17卢廷浩，鲍伏波．接触面薄层单元耦合本构模型．水利学报．2002年第二期18PuvaultGandLionsJLInequalitiesinMechanicsandPhysics。Berlin：Springer．1976198augD，SaxceG．FrictionlessContactofElasticBodiesbyFiniteElementMethodandMathematicalProgrammingTechnique．ComputStruct，11：55～67，198020KlarbringA．GeneralContactBoundaryConditionsAndtheAnalysisofFrictionalSystems，IntJSolidsStruct，22：1377～1398，198621LiM．ShaD，TammaKK．LinearComplementaryFormulationsInvolvingFrictionalContactforE1astoplasticDeform8bleBodies．JApplMech，64：80～89，199722ZhongWanxie．SunSuming．AFiniteElementMethodforElastoplasticStructureandContactProblelDbyParametricQuadraticProgramming,IntJNumerMethEngng,26：2723～2738．199823PianT1tHandKubomuraK．FormulationofContactProblemsbyAssumedStressHybridElements．inNonlinearFiniteE1ementAnalysisinStructuralMechanies。ⅣWanderlichetalEds．Springer—Verlag，198124ShyuSC．ChangTYPandSaleedAF，FrictionContactAnalysisUsingaMixedFiniteE1ementMethod．ComplStruct．V01．32(NO．1)：223-242，198925HeyligerPRandReddyJN．AMixedComputationalAlgorithmforPlanElasticContact 西南交通大学硕士研究生学位论文第68页Problems—IFormulationComputers&Structures，26：62i～634，198726PottsDM．MartinsJP．TheShaftResistanceofAxiallyLoadedPileinClay．Geotechnique，V01．32(3)：145～158，198227E1lisonRD．LoadDeformationMechanismforBoredPiles．J．oftheSollMechanicsFoundationDivision，V01．97(4)：66～677，197128TrochanisAM，BielakJ，ChristianoP．SimplifiedModelforAnalysisofOneorTwoPiles．ASCE，1991，GT3：448～44629DesaiCSNumericalDesign-AnalysisforPilesinSands．J．oftheGeotechnicalEngineeringDivision，V01．100(6)：613～635，197430GoodmanRE，TaylorRL，BrekkeTLAModelfortheMechanicsofJointedRock．．J．oftheSoilMech．AndFoundal，ionsDivision(SM3)。ASME．1968．Vik．94：638～55931Zienkiewicz0CTheFiniteElementMethod．3rded，London：McGraw-Hill．197732朱伯芳．有限单元法原理与应用(第二版)．北京：中国水利水电出版社．199833殷有泉．固体力学非线形有限元引论．北京：北京大学出版社，清华大学出版社．198734DesaiCS，ZamanMM，LightnerJGetal．ThinLayerElementForInterfaceandJoints[J]。Int。Jorun。ForNum．＆Analy．Meth．InGeomech．1982；V01．8(1)：19--4335段文峰，金菊顺，佟德生．在竖向荷载作用下桩一土数值模拟与试验验证．吉林建筑工程学院学报．1997；V01．19(1)：1l～1436孙索泰．弹性力学中的变分原理导引．南京：东南大学出版社．198937卓家寿，陈振雷，颜炳乐．非线性介质摩擦接触问题的变分不等式及其应用．第三届全国计算力学会议论文集一计算力学的理论与应用．北东：科学出版社．1992：100～10338张焱，孔祥安，陈平，施太和．接触力学中的变分原理．谣南石油学院学报．1998；Y01．20(1)：99～10439梁昆淼．数学物理方法(第三版)．北京：高等教育出版社．199540朱百里，沈珠江．计算土力学．上海：上海科学技术出版社．19904l沈珠江．理论土力学．北京：中国水利水电出版社．200042卓家寿．弹性力学中的有限单元法．北京：高等教育出版社．198743D．R．J．欧文，E．辛顿．塑性力学有限元一原理与应用．曾国平，刘忠，徐家礼译．北京：兵器工业出版社．198944MarrWA，ChristianJT．FiniteElementAnalysisofElasto—PlasticSoiis．ReportN0．4．FiniteElementAnalysisofElasto～PlasticFrictionalMaterialsforApplicationtoLunarEarthSciences，June197245潘昌实．隧道力学数值方法．北京：中国铁道出版社．199546雷晓燕．岩土工程数值计算．北京：中国铁道出版社．199947JEAKIN．ApplicationAndImplementatiDnofFiniteElementMethods，AcademicPress198248陆严清．塑性变形理论及应用．北京：国防工业出版社．199849R．D．库克．有限元分析的概念和应用．程耿东，何穷，张国荣译．北京：科学出版社．198950李景渭．有限元法．北京：北京邮电大学出版杜．199951谢康和，周健．岩土工程有限元分析理论与应用．北京：科学出版社．200252吴晓涵．面向对象结构分析程序设计．北京：科学出版社．200253孔祥安，瞿己．面向对象有限元程序的数据设计．西南交通大学学报．1996；V01．3l(4)：355～360 亘塑至夔查堂亟主堕塞圭堂焦迨窒篁!!夏54崔俊芝t梁俊．现代有限元软件方法．北京：国防工业出版社．1995年55曹中清．面向对象的有限元程序设计方法研究．西南交通大学硕士研究生学位论文．199556李会平，曹中消，周本宽．弹塑性分析的面向对象有限元方法．话南交通火学学报，1997；V01．32(4)：401～40657朱晓光．面向对象的非线形有限元程序框架设计．火连理工大学硕士研究生学位论文．200258习俊通，梅雪松，吴序堂．面向对象的接触问题的有限元分析．机械强度．1999：V01．2l(1)：39～4159MarcoCanto．MasteringDelphi5．SybexPress．199960HGPoulos．EHDavis．PileFoundationAnalysisAndDesign．JohnWiley&Sons19806l左名麒，圳人礼，毛洪渊．桩基础工程．中国铁道出版社．199662BrajaM．Das．PrinciplesofFoundationEngineering．Brooks／ColeEngineeringDivision．198463成昆铁路技术总结委员会．成昆铁路2一线路工程地质及路基，人民铁道出版社，198064Zhaoshanrui，AnalysiSandCalculationofState—vectorsForAxiallyLoadedPiles．Proc．OfEASE(I)．Thailand．1986：15～1765陈开旭，安关峰，鲁亮．采用有厚度接触单元对桩基沉降的研究．岩土力学．2000：V01．21(1)：92～9666吴兴序．灌注桩的轴向阻抗特性及压力灌浆技术研究．西南交通大学博士研究生学位论文200167吴起星．广西第三系泥岩桩端承载力研究，广西大学硕士研究生学位论文．200268张靖宇．对桩土相互作用及单桩竖向承载力的研究．西安建筑科技大学硕士研究生学位论文．200169葛崇勋，张永胜．关于嵌岩桩竖向承载力计算方法的探讨．江苏建筑．2001；第四期：47～5270ChangCSandDuncanJM．AComputerProgramForAnalysisofStressesandMovementinExcavations，ReportNo．TE77—4，UniversityofCalifornia，Berkeley7lArgyriSJH，FaustG，SzimmatJ．WarnkeEP．andV／i11iamKJ．RecentDevelopmentsinTheFiniteElementandAnalysiSofPCRV．2ndInt．Conf．SMTRT。Berlin，197372XYLei，GSwobodaAndGZens．ApplicationofContactFrictionInterfaceElementtoTunnelExcavationinFaultedRock．ComputersAndGeotechnicy．1995：V01．17：349～37073DruckerDCandPragerW．SollMechanicsandPlasticAnalysiSorLimitDesign．Q，Appl．Math．10，195274DuncanJMandChangCY．NonlinearAnalysiSofStress—StrainForSoilS，Proc．ASCE，V96，SM5．197075PoscoeKH，SchofiedANandThurairajahA．YieldingofClaysinStatesWeterThanCritical．Geotechnique13(3)，196376RoscoeKHandBurlandJB．0ntheGeneraliezedStress-StrainBehaviourof“Wet”Clay．InEngineeringPlasticity．CambridgeUniv．Press，535．196877Yuzhiqiang，Wangxu，Wuxingxu．GroutedPileandItsBearingCapacity．JournalofSouthWestJiaotongUniversity．2000：V01．8(1)：25～3278VesicAS．Expans．ionofCavitiesinInfiniteSoilMass．Proe．ASCE，V01．98，SM3．】972 ．亘壹至塑查兰亟±堑窒生堂鱼迨窭篁!!里79毛坚强．AnImprovedMethodFortheSolutionofLinearViscoelasticProblemsWith删．JournalofSouthwestJiaotongUniversity．1990：280郑宏，葛修润，谷先荣，丰定详．关于岩土工程有限元分析中的若干问题．岩j：力学．1995：V01．16(3)：7～1281康清梁．钢筋混凝土有限元分析．北京：中国水利水电出版社．199682狄生林．钢筋混凝土梁的非线形有限元分析．南京工学院学报．1984(2)83Xianzhaoran，Anhuining，LeeKM．Non一1inearRnalysiSofaSinglePileBehaviour．JournalofYantaiUniversity．2002；V01．15(3)：201～20684罗战友单桩竖向极限承载力的灰色预测．西安建筑科技大学硕士研究生学位论文．200l85刘兴远，康景文。林文修．桩基工作特性分析的神经网络模型．北京：中圆建筑J：业出版社．199986刘利民，舒翔，熊巨平．桩基工程的理论进展与工程实践，北京：中国建材工业出版社．200287甘德福，陈孝培．论桩侧向压力的黄金分割．工程勘察．1997；第二期88黄强．桩基工程若干热点技术问题．北京：中国建材工业出舨社．199689陶全心，李著璨．结构优化设计方法．北京：清华大学出版社．198590中华人民共和国建设部．建筑桩基技术规范(JGJ94—94)．北京：中国建筑工业出版社．199591LeeIK．ApplicationofFiniteElamentMethodinGeot．Engg．。PartI—LinearAnalysiS，Ch．17inFiniteElementTechniques—AShortCourseofFundamentalSandApplication，Univ．ofN．S．W，。Aust，197392RahmanMU，eta1．AnIteratireProcedureForFiniteE1amentStressAnalysisofFrictionalContactProblems[J]．ComputerStructure，18：947～954，198493AyariML，SaoumaVE．StaticandDynamicContact／ImpactProblemsUsingFictitiousForces[J]．IntJNumerMethEngng，32：623～643，1991 西南交通大学硕士研究生学位论文第71页攻读硕士学位期间发表的论文及科研成果1发表的论文[1]王新．钻孔灌注桩桩底压浆后桩侧摩阻力的增强因素分析．四川建筑．2003年第23卷第六期[2]王新．昔格达泥岩竖向承载力试验及有限元分析．隧道建殴．2003年第23卷第三期：16～18[3]刘智勇，姚勇，王新．昔格达泥岩中冲孔灌注桩的竖向承载力试验分析．话南科技人学学报．已录用2参与的科研项目及试验[1]交大靓典碎石桩复合地基检测【2]九寨沟黄龙机场填方阶段检测[3]金菊花园桩基深层静载试验[4]深圳地铁大轴力桩基托换及重叠隧道施工技术的数值计算[5]西攀高速公路提高桩基承载力试验研究[6]九寨沟黄龙机场跑道地基土反应模量测试[7]高海公路自平衡法桩基承载力试验'