修正剑桥模型

1.剑桥模型（Cam-clay Model ）剑桥模型是由英国剑桥大学Roscoe 等于1963年提出的，这个模型基于正常固结土和超固结土试样的排水和不排水三轴实验基础上，提出了土体临界状态的概念，并在实验基础上，再引入加工硬化原理和能量方程，提出剑桥模型。

这个模型从试验和理论上较好的阐明了土体弹塑性变形特征，尤其考虑了土的塑性体积变形，因而一般认为，剑桥模型的问世，标志着土本构理论发展的新阶段的开始。

（1） 剑桥模型。

剑桥模型基于传统塑性位势理论，采用单屈服面和关联流动法则屈服面形式也不是基于大量的实验而提出的假设，而是依据能量理论提出的。

依据能量方程，外力做功dW 一部分转化为弹性能e dW ，另一部分转化为耗散能（或称塑性能）p dW ,因而有dW =e dW +p dW（1-154） e dW =e eV qd d p γε+' （1-155）p pV p qd d p dW γε+'= （1-156）剑桥模型中，由各向等压固结实验中回弹曲线确定弹性体积变形p p d e k d eV ''+=1ε （1-157）式中，k 为膨胀指数，即 p In e '-回弹曲线的斜率。

同时，假设弹性剪切变形为零，即0=e d γ （1-158）则弹性能 p d ekp p d k dW e '+=''=1υ （1-159）剑桥模型中还建立如下的能量方程，即塑性能等于由于摩擦产生的能量耗散，则有 p p p V d p qd d p γνγε'=+'- （1-160）式中第一项改用负号，是因为p V d ε取以压为正。

代入式（1-161）⎪⎭⎪⎬⎫==ij pij p d s d d λεεθθσ （1-161） 并考虑式（1-158），则有γγγνd p M d p M d p dW p p p '='='= （1-162） 式中，M 为q p '-'平面上的破坏线的斜率，即ϕϕ'-'=sin 3sin 6M （1-163） 式中ϕ'为土体有效摩擦角。

1修正剑桥模型介绍土体本构理论是岩土工程学科的重要基础理论。

随着对土体力学特性的不断深入，塑性理论逐渐被应用于土体本构关系的研究中来。

Roscoe 于1963 年提出著名的剑桥粘土模型，是应用塑性理论的代表，被看做现代土力学的开端，在本构理论研究发展过程中, 各种建模思想不断涌现，出现了各种不同形式的土体本构模型，但弹塑性模型中得到公认的还只有剑桥模型。

现在国际岩土本构的一大发展趋势是又回到剑桥模型，在剑桥模型基础上进行改进和修正，修正剑桥模型是由罗斯科（Roscoe）和伯兰特（Burland）于1968年对剑桥模型作了修正后提出的一个土的弹塑性模型。

主要是对剑桥模型的弹头形屈服面形状作了修正，认为屈服面轨迹应为椭圆。

修正后的模型通常称为修正剑桥模型。

随后又修正了剑桥模型认为在完全状态边界面内土体变形是完全弹性的观点。

认为在完全状态边界面内，当剪应力增加时，虽不产生塑性体积变形，但产生塑性剪切变形。

这可认为是对修正剑桥模型的再次修正。

剑桥模型是英国剑桥大学的Roscoe和Burland根据正常固结粘土和弱超固结粘土的三轴试验，采用状态边界面的概念，由塑性理论的流动法则和塑性势理论，采用简单曲线配合法，建立塑性与硬化定律的函数。

它考虑了静水压力屈服特性、压硬性、剪缩性，但破坏面有尖角，该点的塑性应变方向不易确定。

其假定的弹性墙内加载仍会产生塑性变形。

原始的剑桥模型中存在一个缺点，即p'轴上各向同性压缩的屈服点p'的屈服面正交方x向（塑性流动方向）与水平坐标轴方向不一致。

这会导致各向同性加载（初始固结）所产生的塑性（体积）应变增量方向（它应该与水平坐标p'轴的方向一致）与屈服面的正交方向（塑性流动方向）不一致，如图1所示，图中虚线为原始剑桥模型的屈服面。

这是原始剑桥模型的屈服面与试验结果不一致的地方，也是该屈服面不足的地方。

图1 原始剑桥模型和修正剑桥模型在点处的流动情况纵观剑桥模型40 多年的发展，总结其局限性主要有：（1）受制于经典塑性理论，采用Drucker公设和相关联的流动法则，在很多情况下与岩土工程实际状态不符；破坏面有尖角，该点的塑性应变方向不易确定。

土的弹塑性模型近年来，根据弹塑性理论建立的土的弹塑性模型发展很快，各国学者提出的弹塑性本构模型很多。

下面几节分别介绍剑桥模型，修正剑桥模型，Lade-Duncan 模型，以及清华模型的基本概念。

一．剑桥模型英国剑桥大学Roscoc 和他的同事（1958~1963）在正常固结粘土和超固结粘土试样的排水和不排水三轴试验的基础上，发展了Rendulic （1937）提出的饱和粘土有效应力和孔隙比成唯一关系的概念，提出完全状态边界面的思想。

他们假定土体是加工硬化材料，服从相关联流动规则，根据能量方程，建立剑桥模型。

剑桥模型从理论上阐明了土体弹塑性的变形特性，标志着土的本构理论发展新阶段的开始。

1.临界状态线和Roscoe 面各向等压固结过程中，孔隙比e 或比容()1e υυ=+与有效应力的关系可用下式表示：ln N p υλ'=-（1）式中N ——当 1.0p '=时的比容。

因此exp N p υλ-⎛⎫'= ⎪⎝⎭（2）（a）,p q ''平面（b）,ln p υ'平面图1临界状态线正常固结粘土排水和不排水三轴试验表明：它们有条共同的破坏轨迹，与排水条件无关。

破坏轨迹在,p q ''平面上是一条过原点的直线，在,ln p υ'平面上也是直线，目与正常固结线平行，分别如图（a）和（b〕所示。

破坏轨迹线可用下式表示：cs csq Mp '=（3）ln cs cs p υλ'=Γ-（4）式中CS ——表示临界状态；M——,p q''平面上临界状态线斜率；p'=时土体的比容；Γ—— 1.0csυ'平面上临界状态线斜率。

λ——,ln p一旦土体的应力路径到达这条线，土体就会发生塑性流动。

这时土体被认为处于临界状态，破坏轨迹被称为临界状态线。

临界状态线在,,''空间为一条空间曲线，如下图2所示。

正常固结黏土本构关系研究综述摘要:正常固结黏土是一种普遍存在于土工和岩土工程中的地质材料，在地基工程和地质灾害评估等领域具有重要作用。

本文通过文献综述的方式总结了已有研究中适用于正常固结黏土的本构模型：Mohr-Coulomb模型、Drucker-Prager模型、Cam-Clay模型、Bingham模型，并对其进行分类评述。

关键词：黏土；本构模型；综述在土壤力学领域中，对于黏土的研究一直是一个重要而又复杂的课题。

黏土是黏性土的典型代表，具强塑性、吸水性、膨胀性、收缩性、吸附性、冻胀性、烧结性、耐火性等特殊性质。

其性质和行为对于工程建设、地质灾害预测与调控等方面具有重要意义[1]。

本构关系指的是描述土体力学性质与其应力-应变关系之间联系的数学模型。

对于正常固结黏土的本构关系研究，旨在揭示黏土在径向和剪切方向上的变形特性。

本综述旨在回顾当前关于正常固结黏土本构关系的研究进展，并对其中涉及的主要观点、方法进行综合概述。

1.Mohr-Coulomb模型土塑性力学理论始于1773年法国科学家库伦提出的Coulomb屈服准则1776年库伦总结土的破坏现象和影响因素，提出土的破坏公式为：根据砂土实验结果得到：对于黏性土，可给出更为普遍的表达式：式中c为内聚力，为内摩擦角。

Mohr-Coulomb模型[2]是基于莫尔-库仑准则的一种用于描述土体和岩石的强度和变形行为的本构模型。

该模型假设土体或岩石在发生破坏时遵循弹塑性行为，其中包括两个主要参数：内聚力和内摩擦角。

内聚力越大，材料越难破坏；内摩擦角越大材料越抗剪切破坏。

应力状态：屈服准则：Mohr-Coulomb模型通过在应力空间中绘制Mohr圆来描述材料的破坏准则。

根据Mohr圆上某一点的位置，可以判断材料处于弹性区还是塑性区，以及是否达到破坏条件。

2．Drucker-Prager模型Drucker等提出在Mohr-Coulomb锥形屈服面上再加上一族强化帽形屈服面。

常见地基模型总结常见地基模型总结地基模型是描述地基土在受力状态下应力和应变之间关系的数学表达式。

广义的讲，是描述土体在受力状态下的应力、应变、应变率、应力水平、应力历史、加载率、加载途径以及时间、温度等之间的函数关系。

通常模型有线弹性地基模型、非线弹性地基模型和弹塑性地基模型等。

一、线弹性地基模型地基土在荷载作用下，应力应变关系为直线关系，用广义胡克定律表示。

常用的有三种，温克勒地基模型、弹性半空间地基模型、分层地基模型。

1、温克勒地基模型假定地基由许多独立且互不影响的弹簧组成，即地基任一点所受力只与该点的地基变形成正比，而且该点所受的力不影响该点以外的变形。

表达式为p=k·s（式中k为地基基床系数，根据不同地基分别采用现场载荷班试验或室内三轴、固结试验获得）。

该方法计算简便，只要k值选择得当，可获得较为满意的结果，但在理论上不够严格，未考虑土介质的连续性，忽略了地基中的切应1力，按这一模型，地基变形只发生在基底范围内，而在基底范围外没有地基变形，这与实际不符使用不当会造成不良后果。

该法在地基梁和板以及桩的分析中广泛采用，如台北101大楼采用了广义温克勒地基模型。

由于该模型未考虑剪力作用，故主要使用于土层薄、结构大、土层下为基岩（剪切模量小、可压缩层薄）的地基，而上硬下软的地基不适用。

2、弹性半空间地基模型假定地基为均匀、各向同性的弹性半空间体。

采用Boussinesq公式求解。

对于均布荷载下矩形中点的竖向变形以及对于荷载面积以外的任一点的变形可以通过积分求得。

该法考虑了压力的扩散作用，比温克勒模型更合理，但未反应地基土的分层特性，且认为压力可以扩散到无限远处，造成计算的沉降量和地表沉降范围都较实测结果为大。

3、分层地基模型分层地基模型即是我国地基基础规范中用以计算地基最终沉降量的分层总和法。

该模型能较好的反应地基土扩散应力和变形的能力，能较容易的考虑土层非均匀性沿深度的变化和土的分层，通过计算表明，分层地2基模型的计算结果比较符合实际情况。

万方数据　万方数据　万方数据第３期关云飞等：ＡＮＳＹＳ软件中修正剑桥模型的二次开发（１２）更新塑性应变增量｛ｄ￡），判别本增量步是否收敛，若不收敛，将时间二分，重新在该高斯点求解。

（１３）程序结束。

３．３利用ＵＰＦｓ对修正剑桥模型进行二次开发本文利用二次开发工具ＵＰＦｓ将修正剑桥模型添加到标准ＡＮＳＹＳ程序中，具体做法如下：（１）安装ＡＮＳＹＳ—ＵＰＦｓ二次开发工具以及Ｆｏｒｔｒａｎ编译器，不同ＡＮＳＹＳ版本对于编译器的版本有相应的要求。

本文基于ＡＮｓＹＳ７．０，在ＷｉｎｄｏｗＸＰ下，使用ＣｏｍｐａｑＶｉｓｕａｌＦｏｒｔｒａｎ６．６Ｂ编译器。

（２）根据前述的方法编写修正剑桥模型弹塑性本构关系的用户子程序ＵＳＥＲＭＡＴ．Ｆ。

（３）将本构模型用户子程序ＵＳＥＲＭＡＴ．Ｆ连接到ＡＮＳＹＳ中。

利用ＵＰＦｓ和数据接口生成含修正剑桥模型的ＡＮＳＹＳ程序。

（４）运行自定义版本的ＡＮＳＹＳ程序。

有两种方式运行二次开发后的ＡＮＳＹＳ程序。

一是用第（３）步生成的可执行文件覆盖安装目录下的ＡＮＳＹＳ．ｅｘｅ，然后按正常方式启动；二是利用ａｎｓｙｓ７０ｃｕｓｔ命令调用，运行”ａｎｓｙｓ７０ｃｕｓｔ．ｃｕｓｔｏｍ／ｐａｔｈｎａｍｅ／ａｎｓｙｓ．ｅｘｅ”，进入ＡＮＳＹＳ界面。

在激活ＵＰＦｓ运行自定义版本的ＡＮＳＹＳ程序后，使用ＴＢ，ＵＳＥＲ命令通知程序将使用自定义本构关系，并通过命令ＴＢＤＡＴＡ将本构模型的相应参数输入。

由于本文在修正剑桥模型的子程序中使用了状态变量，包括各方向的塑性应变，总的等效塑性应变及等效应力，因此，在定义材料属性时必须使用ＴＢ，ＳＴＡＴＥ命令定义状态变量的个数和大小，这里设置状态变量个数为６，状态变量的初始值为０。

进行有限元计算时，在每次Ｎｅｗｔｏｎ．Ｐａｐｈｓｏｎ迭代过程中，修正剑桥模型的子程序在每个单元积分点被调用。

子程序通过在每次时间增量开始时该积分点的应力、应变、状态变量及当前的应变增量进行计算，得出本次时间增量结束时的应力、应变和各状态变量的值，并输出相应的弹塑性刚度矩阵。

（1）计算理论（2）计算所使用的数据注意：在运行代码时，应将数据复制为txt文件，保存在代码所在目录，并命名为”number.txt”。

e0landa ka M miu p q pc00.50.140.05 1.250.2100.00330.011000.50.140.05 1.250.2101.666751000.50.140.05 1.250.2103.3333101000.50.140.05 1.250.2105151000.50.140.05 1.250.2106.6667201000.50.140.05 1.250.2108.3333251000.50.140.05 1.250.2110301000.50.140.05 1.250.2111.6667351000.50.140.05 1.250.2113.3333401000.50.140.05 1.250.2115451000.50.140.05 1.250.2116.6667501000.50.140.05 1.250.2118.3333551000.50.140.05 1.250.2120601000.50.140.05 1.250.2121.6667651000.50.140.05 1.250.2123.3333701000.50.140.05 1.250.2125751000.50.140.05 1.250.2126.6667801000.50.140.05 1.250.2128.3333851000.50.140.05 1.250.2130901000.50.140.05 1.250.2131.6667951000.50.140.05 1.250.2133.33331001000.50.140.05 1.250.21351051000.50.140.05 1.250.2136.66671101000.50.140.05 1.250.2138.33331151000.50.140.05 1.250.21401201000.50.140.05 1.250.2141.66671251000.50.140.05 1.250.2143.33331301000.50.140.05 1.250.21451351000.50.140.05 1.250.2146.66671401000.50.140.05 1.250.2148.33331451000.50.140.05 1.250.21501501000.50.140.05 1.250.2151.66671551000.50.140.05 1.250.2153.33331601000.50.140.05 1.250.21551651000.50.140.05 1.250.2156.66671701000.50.140.05 1.250.2158.33331751000.50.140.05 1.250.21601801000.50.140.05 1.250.2161.66671851000.50.140.05 1.250.2163.33331901000.50.140.05 1.250.21651951000.50.140.05 1.250.2166.66672001000.50.140.05 1.250.2168.33332051000.50.140.05 1.250.21702101000.50.140.05 1.250.2171.6667215100（3）python 相关代码import numpy as npfrom math import *import matplotlib.pyplot as pltfrom matplotlib import animationdef calculate_pq(a1,a2,a3):#计算pq(三维)p = (a1+a2+a3)/3q1 = sqrt(((a1-a2)**2)+((a2-a3)**2)+((a3-a1)**2))q = sqrt(0.5)*q1return p,qdef calculate_matrix(e0,a,b,m,u,p,q,p0):#计算模型柔度矩阵v = 1+e0l = q/pDp = b/(p+p*e0)+((a-b)*((m**2)-(l**2)))/((p+p*e0)*((m**2)-(l**2)))Dpq = ((a-b)*2*l)/((p+p*e0)*(m*m+l*l))Dq = (2*b+2*b*u)/(9*v*p-18*u*v*p)+((a-b)*4*l*l)/((p+p*e0)*((m**4)-(l**4)))D = np.array([[Dp,Dpq],[Dpq,Dq]])return Ddef calculate_drained_pc_sigma(e0,a,b,m,u,p,q,p0):pc = p+q*q/(m*m*p)sigmav = ((a-b)*log(p/p0))/(1+e0)+((a-b)/(1+e0))*log(1+(q*q/(m*m*p*p)))return pc,sigmavdef plot_data(fr):xdata.append(pc_x[fr])ydata.append(sigmav_y[fr])line.set_data(xdata,ydata)return line,def plot_data2(fr):xdata2.append(sigmas[fr])ydata2.append(q_y[fr])line2.set_data(xdata2,ydata2)return line2,def plot_data3(fr):xdata3.append(sigmas[fr])ydata3.append(sigmav[fr])line3.set_data(xdata3,ydata3)return line3,#计算硬化规律相关数据number = np.loadtxt('number.txt',skiprows=1)pc_x = []sigmav_y=[]for i in range(len(number)):arr = number[[i]]arr1 = list(arr.flatten())e0,a,b,m,u,p,q,p0 = arr1pc,sigmav = calculate_drained_pc_sigma(e0,a,b,m,u,p,q,p0) pc_x.append(pc)sigmav_y.append(sigmav)#ev,ep,q数据计算number1 = number[:,[5,6]]sigmav = []sigmas = []q_y = []for i in range(len(number1)-1):x = i-1arr = number1[[i]]pi,qi = list(arr.flatten())q_y.append(qi)if i > 0:px,qx = list(number1[[x]].flatten())dp = pi-pxdq = qi-qxelse:dp = pi-100dq = qidpq = np.array([dp,dq])arr1 = number[[i]]e0,a,b,m,u,p,q,p0 =list(arr1.flatten())D = calculate_matrix(e0,a,b,m,u,p,q,p0)sigma = np.dot(D,dpq)dvi,dsi = list(sigma.flatten())if i > 0:vi = sigmav[-1]+dvisi = sigmas[-1]+dsielse:vi = dvisi = dsisigmav.append(vi)sigmas.append(si)result_1 = np.array([q_y,sigmav,sigmas])result1 = np.round(result_1.T,6)np.savetxt('result1.txt',result1,fmt='%10.8f',delimiter='\t\t')#画图fig = plt.figure(figsize=(15,8),dpi=100)plt.rcParams['font.sans-serif'] = ['SimHei']xdata,ydata = [],[]xdata2,ydata2 = [],[]xdata3,ydata3 = [],[]plt.subplot(1,3,1)e('grayscale')plt.axis([0,400,0,0.08])plt.title('硬化规律(三轴压缩)',fontsize=20)plt.xlabel('屈服压力',fontsize=10)plt.ylabel('塑性体积应变',fontsize=10)line, = plt.plot(xdata,ydata,linewidth=1,color='red',marker='o',linestyle='--',markersize=3,label='kexiv-pc' )ani = animation.FuncAnimation(fig=fig,func=plot_data,frames=range(len(sigmav_y)),interval=100,repe at=False,blit=True)plt.legend()plt.subplot(1,3,2)e('grayscale')plt.axis([0,0.4,0,220])plt.xlabel('剪应变',fontsize=10)plt.ylabel('偏应力',fontsize=10)line2, = plt.plot(xdata2,ydata2,linewidth=1,color='red',marker='<',linestyle='--',markersize=3,label='q-kexi _s')ani2 = animation.FuncAnimation(fig=fig,func=plot_data2,frames=range(len(sigmas)),interval=100,repea t=False,blit=True)plt.legend()plt.subplot(1,3,3)e('grayscale')plt.axis([0,0.4,0,0.11])plt.xlabel('剪应变',fontsize=10)plt.ylabel('体应变',fontsize=10)line3, = plt.plot(xdata3,ydata3,linewidth=1,color='red',marker='<',linestyle='--',markersize=3,label='kexi_v -kexi_s')ani3 = animation.FuncAnimation(fig=fig,func=plot_data3,frames=range(len(sigmas)),interval=100,repea t=False,blit=True)plt.legend()plt.show()。

基于修正剑桥模型模拟理想三轴不排水试验——两种积分算法的对比分析（CZQ-SpringGod ）1、修正剑桥模型在塑性功中考虑体积塑性应变的影响，根据屈服面一致性原则，假定屈服函数对硬化参数的偏导为0，就获得了以理想三轴不排水试验为基础的修正剑桥模型屈服函数：22(,)()0c q f p q p p p M =+-= （1） 其中3kkp σ=，ij ij ij s p σδ=-，212ij ij J s s =，q =M 为临界线斜率，c p 为前期固结压力。

硬化/软化法则：p c v c dp v d p ελκ=- （2） 式中p v ε为体积塑性应变，v 为比体积，λ为正常固结线斜率，κ为回弹线斜率。

由于不排水屈服面推导过程是基于硬化参数c p 偏导为0，也就是说不排水试验中硬化参数同体积塑性应变无关，屈服面不变化，而若引入硬化法则就同屈服面推导过程中的假定矛盾，因此计算时将模型处理为理想塑性模型。

2、显式和隐式两种积分格式考虑应变增量ε∆驱动下，第n 增量步到第n+1增量步之间的应力积分格式。

显式积分格式的推导参考文献[1]，其中弹塑性矩阵中的塑性硬化模量H=0。

隐式积分格式推导如下：11()n n n p v v p p K εε++=+∆-∆ （3） 111(2)n p n n v c p p ε+++∆=Λ⋅- （4） 12()n n p ij ij ij ij s s G e e +=+∆-∆ （5） 1123n ij p n ij s e M ++∆=Λ （6） 111112(,)()0n n n n n c qf q p p p p M +++++=+-= （7）在这一组方程中没有硬化规律方程表明为理想塑性，并将式（3）-（7）合并化简得到：1112112122(2)06()(1)0n n n n v c n n n trial c p p K K p p G q p p p M Mε++++++⎧--∆+⋅Λ⋅-=⎪⎨+-+Λ=⎪⎩ （8） 式中3(2)(2)2n n trial ij ij ij ij q s G e s G e =+∆+∆ 求解（8）式方程组即可得到n+1增量步的各个增量。