土的本構理論與模型研究的困難與探索

土的本構理論與模型研究的困難與探索

楊光華1，2

（1. 廣東省水利水電科學研究院，廣東 廣州 510610；2. 廣東省巖土工程技術研究中心，廣東 廣州 510610）

土的本構模型是土力學未來發展的核心，也是現代土力學發展的瓶頸。自從1963年劍橋大學Roscoe教授等在塑性力學理論基礎建立了土的劍橋彈塑性本構模型以來，50多年過去了，雖然取得了很大的進步，建立了數以百計的本構模型或更多，但其真正在工程設計中能像規范標準這樣作為設計依據的應用還很少，多數還只是作為工程設計的參考?，F代科技的發展，已解決了數學計算的困難，現代土力學理論的發展，關鍵還是要解決好土的本構模型，提高工程計算的準確性，使設計更科學合理。土的本構理論與模型歷經了半個多世紀的發展，目前也是遇到了發展的瓶頸。

那么，接下來該如何發展呢？制約發展的困難何在？突破口在哪里？是值得我們認真思考的一個問題。為此，對多年的研究探索做了一下梳理和思考，希望能拋磚引玉，促進現代土力學的發展：

一．什么是土的本構模型？

二．如何開展研究的？

三．有什么困難？

四．應如何發展？

一．什么是土的本構模型？

狹義的土的本構模型是指土的應力與應變關系方程。對于分量而言就是六個分量應力與應變的關系：

通常土是非線性材料，一般用增量應力和應變關系表示：

土的本構模型的研究就是建立（1）或（2）式的應力應變關系方程。

二．如何進行土的本構模型的研究？

以上式（1）或式（2）是六維的應力分量與應變分量的關系，如果我們通過試驗，獲得了6個應力分量和6個應變分量對應關系的試驗結果，則對試驗結果采用數學擬合或統計就可以建立其數學關系方程了。但實際上，我們不能通過試驗直接得到六個分量應力與應變的關系，試驗最多只能在主空間上實現，所以通常最多也只能得到主空間上3個主應力與3個主應變的關系：

這樣（3）或（4）式總是可以通過試驗所得到的。但當我們進行數值計算時所需的是（1）或（2）式的六個分量應力與應變的方程。如何由主空間上的（3）或（4）式得到所需的一般空間下的（1）或（2）式呢？其實這是一個數學問題。通常所說的本構理論就是在各種假定條件下建立的由主空間（3）、（4）式推廣到一般坐標系下（1）、（2）式的數學表達式。

因此，所謂的本構模型的研究要解決的問題就是兩個：一個就是如何通過試驗合理方便的確定主空間的本構方程（3）、（4），因為即使是三維空間，也有很多變化，有不同的應力路徑、不同的應力或應變狀態的組合，不同的組合得到的關系也可能是不同的，這樣就要通過試驗發現其本構的特性，選擇合適的試驗方案，獲得試驗曲線，然后用數學擬合或統計的方法獲得主空間的應力應變關系方程（3）、（4）式。當研究的材料比較簡單時，則可以用簡單的試驗確定其關系，如對于各向同性線性材料，則應力與應變的關系即為廣義虎克定律：

這種材料只有彈性模量E和泊松比V兩個參數，在主空間上試驗確定這兩個參數就可以了。

二是建立各種由主空間到一般坐標空間的轉換方法，即構成了各種本構理論。如彈塑性的塑性位勢理論，各種非線性彈性理論等。

或也可以采用直接的數學坐標變換方法等。

所以，土的本構模型的研究其實就是兩個問題：一是通過試驗解決主空間上的關系，二是建立由主空間到一般坐標空間的變換關系，這種變換其實就是本構理論（楊光華，1990）。

三．土的本構模型研究的困難和解決方法

1.主空間規律性的研究

土的變形特性的復雜性只能通過試驗的方法來揭示，目前已知的主要特性就有：非線性、彈塑性、壓硬性、剪脹性、應力路徑相關性等。

主空間規律性的研究就是通過試驗揭示其更多的力學特性，如壓硬性、剪脹性，非線性等。這些都是巖土材料特有的特性。

主空間上通過試驗獲得規律，然后用數學擬合或統計的方法表示其規律，建立擬合方程。不同的數學擬合方法與不同的本構理論的結合，則得到不同的本構模型。主空間的數學擬合方法有直接擬合法，通常較多的是用經驗函數去擬合試驗結果，典型的如Duncan-Chang模型，是建立于廣義虎克定律基礎上，采用經驗函數雙曲線函數擬合試驗曲線，然后由擬合曲線方程求取廣義虎克定律的兩個參數Ｅ、ν。彈塑性模型通常是依據彈塑性本構理論的要求，通過主空間上的試驗去確定塑性勢函數、硬化規律等，當假設不同的塑性勢函數或用不同的硬化函數時，則得到不同的彈塑性模型。

2.本構理論的研究

各種本構理論其實就是一個數學上的坐標變換方法，這是本構理論研究的一個重要的新認識。

楊光華（1991）從張量分析的角度得到，當兩個二階張量的主方向一致時，則存在數學上的變換關系：

反之：

為二階張量的三個主值，當應力與應變主方向一致時，則可有

另一方面，也可以從數學矢量擬合的角度建立用勢函數表述的坐標變換關系。假設在一個三維坐標空間上表達一個已知的應變矢量，在應力空間可以采用三個線性無關的勢函數的梯度矢量來擬合這一個已知矢量：

      當與主軸方向相同時，則由（9）式有

      把（10）式代入(11)式，則得到：

      如果矢量在主應力空間上是一個有勢場，則可以用單一的一個勢函數表示：

      將其代入（11）式則由上式得到

      這就是傳統的彈性位勢理論，也稱為格林（Green）超彈性理論。以往多是從能量角度推導得到，Q即是彈性勢能函數。而這里是從數學的方法得到，從數學方法得到我們就可以清楚其作了什么樣的數學假設，其實它是一定數學假設條例下的特殊情況，這樣，應用時就清楚了其適用的條件和范圍。 

      在式（12）中，可以表示為三個應力不變量的函數，展開得到

      進一步展開就是通常的柯西（Cauchy）本構理論，同樣，Green理論也是其一種特例。

      假設對于塑性應變增量也把其看作為應力空間上的一個三維矢量。同樣可以用三個線性無關的勢函數的梯度矢量來擬合

      假設塑性應變增量與應力總量主方向一致，則由以上的分解準則（9）式可以有

      把（16）式代入上式（17）中，則有

      這就是楊光華（1991）提出的廣義（塑性）位勢理論之一，當然主空間的廣義位勢理論也可以在應變空間中表述，也可以用塑性應力增量表述：

      這樣，式（18）到式（21）就構成了一套系統對偶的廣義位勢理論體系。當是一個有勢場矢量時，則可以采用單一的一個勢函數Q的梯度矢量擬合，

      同樣，假設塑性應變增量與應力總量主方向一致，則由以上的分解準則（9）式可以有

      把（22）代入（23），則可得到：

      這就是通常的塑性位勢理論，Q即是塑性勢函數，是廣義位勢理論的特例，這里也是通過數學的方法而得到的，與傳統的塑性位勢理論得到的方法不同。從數學方法上，則我們可以知道，傳統的塑性位勢理論是有兩個數學假設的：一是一個有勢場矢量，二是假設了塑性應變增量方向與應力總量主方向是一致的?？梢?，我們常用的彈塑性本構理論其數學背景是一個特殊的數學坐標變換方法。同樣也可有應變空間上的數學關系，如：

同樣，對于塑性應力增量也有應變空間和應力空間的正交關系：

（24）—（27）就是楊光華（1988）提出的廣義塑性位勢理論。

這就是本構理論的數學原理。本構理論的實質其實就是解決主空間到一般坐標空間的一個數學變換的方法，傳統的各種本構理論都可以這樣用數學的方法統一來得到。數學變換主要是有兩種形式，一種是應力與應變總量的主方向一致的分解準則，或者是增量與增量主方向一致的分解準則，或稱同級分解準則，另一種是增量與總量主方向一致的分解準則。廣義虎克定律是最簡單的同級線性分解準則，數學背景是應力增量和應變增量主方

向相同的最簡單形式。傳統塑性位勢理論是假設了塑性應變增量與應力總量主方向一致的數學條件而得到的，是增量與總量主方向一致的分解準則。

四． 如何構建土的本構模型？

目前應用較多的模型主要是兩種，一是建立于廣義虎克定律基礎上的增量非線性模型，經典的如Duncan-Chang模型；二是建立于經典彈塑性理論基礎上的彈塑性模型，經典的如劍橋模型。

以廣義虎克定律為基礎的各種模型均是通過主空間的試驗求兩個參數，用廣義虎克定律即可得到六個應力與應變分量關系的本構方程。關鍵是采用什么的試驗和用什么樣的函數來擬合試驗曲線。

彈塑性模型過去多是設法去推求塑性勢函數和求取硬化參數。傳統塑性理論中通過Drucker公設，得到屈服函數與塑性應變增量方向具有正交關系，與塑性位勢理論的數學條件一致，這樣就可以用屈服函數代替塑性勢函數。因為對于金屬材料，其具有明顯的屈服點，屈服函數易確定，但對于巖土材料，沒有明顯的屈服點，難以明確確定屈服函數，于是就只能利用數學正交條件，由塑性應變的增量方向反推求塑性勢函數，求得塑性勢函數和和塑性勢函數隨塑性應變而變化的規律，即硬化規律，再用彈塑性本構理論構建彈塑性矩陣。最典型的是劍橋模型和清華模型，但這種求塑性勢函數的方法除劍橋模型相對簡單外，其他各種彈塑性模型都過于復雜，工程實踐中較少得到應用。同時這種建模方法對土的剪脹反映不太夠理想。廣義虎克定律也不能反映土的剪脹。

廣義虎克定律采用擬合試驗曲線的方法求模型參數，簡單直觀，彈塑性模型表述能力強但復雜，硬化土模型（Harding soil model）是結合了兩者的優點的模型，用雙曲線表述試驗曲線，用彈塑性理論構建剛度矩陣。

1. 用廣義位勢理論構建統一的本構模型

有限元等數值方法最終需要的是六維應力與應變分量關系的計算矩陣。按廣義位勢理論有：

勢函數為應力不變量的函數，忽略的影響，則為，代入（28）式得

這樣，并不需要通常的塑性勢函數或數學勢函數。

設主空間上試驗擬合所得方程為：

同時通過回彈試驗得到彈性矩陣為，則由（29）~ (31)建立六維的本構關系，可以直接利用四個參數建立統一的本構模型計算矩陣：

這樣，只要有四個參數即可以建立彈塑性的計算矩陣，不需要像通常的彈塑性模型那樣去推求塑性勢函數和硬化參數等。

顯然，上式在一定的條件下又可退化為傳統的彈塑性本構模型。例如：

當，則退化為傳統的關聯流動彈塑性模型。劍橋模型、修正劍橋模型滿足這個關系。

當，則退化為傳統的非關聯流動彈塑性模型。

數學上的意義就是p-q空間上塑性應變的二階矩陣是非滿秩的矩陣，物理力學意義就是塑性應變增量方向與應力增量方向無關。

若，則表示塑性應變增量方向與應力增量方向有關，說明模型可以反映塑性應變增量方向與應力增量方向相關的特性。

可見，廣義位勢理論用于建立彈塑性模型更具一般性，且無需傳統的彈塑性理論那樣需要去推求塑性勢函數的麻煩，可以統一建立計算彈塑性矩陣，當滿足一定的數學條件時，則退化為傳統彈塑性理論的簡單模型，可見其能表述更廣泛的本構特性。

還可以由（21）式構建應變空間的本構模型的統一計算矩陣（楊光華，1996），本構模型的參數也是。

2.主空間上廣義位勢模型參數的確定

利用廣義位勢理論建立本構模型，最后都是回到主空間上確定模型參數。以p-q空間為例，最后表達式都是相同的形式：

考慮彈性變形，則變成：

      為線彈性體積模量和剪切模量。

或寫成

廣義虎克定律時：可由或表達。這樣，只要通過主空間試驗確定

這四個參數即可以。關鍵問題是：

用什么樣的應力路徑試驗或試驗方案來確定這幾個系數呢？

土存在應力路徑相關性，不同的應力路徑試驗獲得的參數有差異，而實際工程中各點的應力路徑可能都不同，但是不可能對所有點的應力路徑進行各種試驗。采用的試驗方案應是：一是要反映主要應力路徑；二是試驗要簡單易實現；三是取決于模型參數的多少，當模型的參數多時，則要采用更多的不同應力路徑的試驗才能確定?？梢杂幸韵乱恍┓椒ǎ?/span>

1） 擬合試驗方法

由p-q空間的試驗，可以擬合得到以下方程：

對其求導，可得：

這樣，對比（36）-（39）即可以得到這4個參數。

2). K-G模型法

通常的K-G模型可以得到以下的一般關系：

如果進一步分解出塑性應變增量

則可以得到：

這個方法可以通過反映土的剪脹性。

3). 模型法

對于彈塑性模型，彈性參數或可以用卸載試驗確定，對于這4個參數，相當于要有4個獨立方程才能確定。但當采用一定的假設后，就可以簡化，如假設符合關聯流動法則，則可獲得兩個已知方程：，這樣4個參數只需再補充兩個方程，這樣就由4個方程可以確定4個獨立參數了，這時可以像Duncan-Chang模型一樣，采用. 的常規三軸試驗并擬合試驗結果，可以得到兩個參數，補充了兩個方程，其參數確定與Duncan-Chang模型相同，此時與Duncan-Chang模型不同的是，這里可以允許泊松比，從而可以表述土的剪脹特性；當由常規三軸試驗得到兩個參數時，由常規三軸試驗的應力條件，由（46）、（48）可以得到（47）、（49）的兩個補充方程：

聯合這兩個方程求解，即可以得到這四個參數（楊光華，1993）。

當只有時，則要增加一個試驗。同樣，如對這4個參數不做任何假設限定，則需要有4個獨立參數的試驗。      

4).類劍橋模型方法

當假設這4個參數滿足關聯流動法則時，則有：

由這兩個已知方程可以減少兩個未知參數，引入參數β，

這樣（34）、（35）式可以改寫為：

這樣只有A和β兩個參數。參考劍橋模型的方法，這兩個參數可以用三向壓縮試驗和破壞線上的塑性變形特點來確定。

對軸上，，由壓縮試驗得： ,臨界線CSL上，要求：A=0，可以構造。

對：P軸上，，A已知，則必須β=∞，臨界線CSL上，要求β=0這樣構造可以滿足。幾何意義如圖所示。

這樣就可以得到A、β兩個參數的一種形式：

代回（52）式即得到：

從而可以得到這4個參數（楊光華等，2013）。

修正劍橋模型相當于取

，其滿足了數學條件：p軸上M=η時，A=0，β=0。其相應p-q空間上的關系式：

而原劍橋模型則相當于?。?img src="/Upload/images/201907/5d3963d4ad57c.png" title="圖片未命名" alt="圖片未命名" /> 。顯然，原劍橋模型不滿足p軸上η=0時β=∞的條件，即在p軸上有，這是不符合試驗結果的。其相應p-q空間上的關系式：

這樣，通過劍橋模型，我們也可以得到這4個參數。同時，在這里，我們無需像傳統劍橋模型那樣通過假設不同的能量函數來推導建立了，而可以從數學上更直觀方便的建立，并且可以更清楚劍橋模型的數學背景。

這樣，按廣義位勢理論，只要通過p-q空間的試驗，然后擬合試驗曲線，確定彈性矩陣和A、B、C、D這四個參數，然后直接代入（35）式即可得到土的本構模型方程的計算矩陣，直觀方便、理論明確，不需要像傳統的彈塑性理論那樣去推求塑性勢函數和硬化參數等。這樣不僅簡單方便，并且可以建立更具普遍性的本構模型。

這樣，按廣義位勢理論建立模型，土的本構模型的研究就剩下如何選擇主空間的試驗方案和采用什么樣的數學方程擬合主空間的試驗曲線來確定A、B、C、D這四個參數的問題了。

而依據傳統彈塑性本構理論所建立的模型中，由于不同的人采用了不同的塑性勢函數和不同的硬化參數，從而產生了眾多的不同模型。同時，由于傳統的彈塑性本構理論是對A、B、C、D這四個參數作了一些假設限定的，如關聯流動模型假設了：，非關聯流動模型假設了：，但實際土的試驗結果可能是不符合這些假設的，如果還是采用這些理論去表述，顯然是不合理的，但又希望能使表述結果與試驗結果更接近，這就勢必會構造更復雜的模型，試圖去符合試驗結果，就如人們想用廣義虎克定律構建的模型去表述土的剪脹性那樣，勢必會構造復雜而別扭的模型，但結果并不理想，因為土的變形特性已經超出了理論的表述能力了，而科學正確的做法就是創造具有更好表述能力的新理論，才能更好的解決問題。

3.土的本構模型研究的困難：

1） 目前剪脹性表達沒有很好解決

以廣義虎克定律為基礎的模型限定了泊松比，是不能反映土的剪脹性的，土的剪脹是指泊松比的情況。傳統的劍橋彈塑性模型主要是表述剪縮特性的，也是不能反映土的剪脹性的，后人做了一定的擴展，在一定程度上可以反映剪脹性。也有在Zienkiewicz廣義塑性模型的基礎上建立一些考慮剪脹性的新模型，直接用塑性流動方向建模，避免了求塑性勢函數，但其本質還是傳統的彈塑性本構理論。其實在廣義位勢理論基礎上可以很方便的建立剪脹模型，可以直接數學擬合p-q空間上具有剪脹性的試驗曲線，獲得A、B、C、D四個參數即可以建立剪脹模型，當在廣義位勢理論基礎上引入傳統理論的關聯流動假設時，即假設AD-BC=0，B=C，這時只需要像Duncan-Chang模型那樣確定切線模量Et和泊松比即可建立模型，當出現剪脹時則泊松比可以，從而可以簡便有效的建立反映剪脹性的模型。

2） 原狀性沒有解決

除了水利工程的筑壩材料是人工填筑的材料外，大部分土木工程遇到的是天然形成的巖土材料。巖土本構模型通常是通過鉆孔取樣進行室內試驗確定模型的參數，這就存在著取樣擾動的影響，使室內土樣與現場土樣不同，由此獲得的模型參數與現場土參數存在差異，從而影響模型的準確性。同樣，即使是大壩的人工填土，現場碾壓的結果與室內試驗模擬的也可能存在差異，對于通常的堆石壩，則室內試驗的顆粒尺寸要遠小于現場碾壓的石料尺寸，室內也很難模擬真實尺寸和碾壓的情況，這些都影響模型參數與實際情況的差異。有必要發展適合于現場原位試驗的本構模型，或通過原位試驗和工程實測結果率定本構模型的參數，可能是提高模型計算準確性的途徑。

3） 選擇合適的試驗確定A、B、C、D四個參數

由于巖土材料的復雜性，影響因素多，如何選擇合適的試驗方案來確定這四個參數，可能要根據工程土體的受力狀態，應力路徑等綜合設計試驗方案，這可能要根據具體工程情況來考慮會更合適，也是需要進一步研究的問題。另一方面，土的復雜性是否還存在人們沒有認識到的特性？如果還有人們沒有認識到的特性，一旦工程中遇到會出現超出人們認識的特性時，也會存在風險，因此，也還需要深入試驗，更多的揭示和認識土的變形特性，如主軸旋轉的影響、次生各向異性、小應變特性等。

4) 缺乏邊值問題的驗證

由以上的研究可知，土的本構模型研究分為主空間的研究和本構理論的研究兩部分，目前通常的模型驗證主要是在主空間上來驗證比較模型預測與試驗結果，這實際上只是驗證了主空間的研究，對本構理論的假設不能驗證。應該采用邊值問題來比較，才可以驗證兩部分。但對于實際工程，由于現場地質分布的不均勻性和現場土體的結構性，用取樣室內試驗確定模型參數用于現場原位邊值問題，也還是不一定合適?？赡懿捎檬覂冗呏的Ｐ驮囼瀬眚炞C會比較合適，均勻性也容易實現，土樣與模型土的一致性也易控制一些。因此，如能組織邊值模型試驗用于驗證可能比現在通常用主空間的土樣試驗驗證更有意義一些。

五、結論

1. 土的本構模型的研究分為兩個內容：一是主空間的試驗和試驗結果的數學擬合，二是把主空間的試驗結果轉換為數值計算所需的六個應力分量與應變分量的關系，第二個問題即是本構理論問題。

2. 土的本構理論與模型研究的困難在于我們對土的變形特性和規律的認識還不夠全面，同時以往采用傳統的本構理論只是一種特殊理論，用于表述土的變形特性的能力不足，需要發展和采用新的更具普遍性的理論。

3. 從數學原理角度建立的廣義位勢理論是更具普遍性的本構理論，解決了由主空間到一般空間轉換的本構理論問題，傳統的本構理論只是特例。從數學原理出發建立的本構理論可以更清楚所采用的本構理論所作的數學物理假設，明確其適應條件，如傳統的以塑性位勢理論為基礎的彈塑性本構理論，其作了兩個假設：一是假設了塑性應變增量矢量場是一個有勢場，二是假設增量塑性應變主方向與應力全量主方向相同，是廣義位勢理論的特例。

4. 可以采用廣義位勢理論建立一個更普遍適用的本構方程，這樣本構模型的研究就剩下如何由主空間確定四個模型參數A、B、C、D的問題了，從而使本構模型的構建既簡化又更具普遍性。

5. 要使模型能用于工程設計，需要提高模型計算的準確性，要解決模型參數確定的可靠性的問題。由于通常的模型參數是由室內試驗確定的，存在取樣擾動和室內試驗與現場樣品的差異等的影響，會使室內土樣試驗得到的參數與現場土存在差異，要提高模型參數的可靠性，需要發展原位試驗確定或率定模型參數的方法。同時，要使模型能在工程中得到推廣應用，還需要構建參數少易確定精度好的實用模型。  

6. 本構模型的驗證，僅只是在主空間進行驗證還是不夠的，全面的驗證應該要采用邊值問題來驗證，即要驗證主空間上的規律和本構理論兩個方面。

7. 還需要深入試驗研究，全面認識和揭示土的變形特性。

主要參考文獻：

[1] 楊光華,李廣信,介玉新著《土的本構模型的廣義位勢理論及其應用》[M]，北京：中國水利水電出版社，2007年

[2] 李廣信 《高等土力學》[M]，北京：清華大學出版社，2002.

[3] 楊光華. 建立彈塑性本構關系的廣義塑性位勢理論[A]. 第三屆全國巖石力學數值計算與解析方法討論會論文集[C]. 1988

[4] 楊光華. 巖土類工程材料應力～應變本構理論的基本數學問題[A]. 巖土力學數值方法的工程應用——第二屆全國巖石力學數值計算與模型實驗學術研討會論文集[C]. 1990

[5] 楊光華，巖土類材料的多重勢面彈塑性本構模型理論，巖土工程學報，1991年第5期。

[6] 楊光華. 巖土類工程材料本構關系的勢函數模型理論[A]. 第四屆全國巖石力學數值計算與解析方法討論會論文集[C]. 1991

[7] Yang Guanghua, A New Elasto-plastic Constitutive Model for Soils, lnt. Conf. On Soils Eng., Science Press, Nov., 1993. Guang zhou, China.

[8] Yang Guang-hua. A New Strain Space Elastoplastic Constitutive Model for Soils. Proceedings of the 2nd Conf. on Soft Soil Eng, May,1996. Nanjing.

[9] 楊光華，一個無需塑性勢函數的新的土體彈塑性本構模型，《水工結構工程理論與應用》，大連海運學院出版社，1993年6月。

[10] 楊光華，巖土類工程材料本構方程的一個張量普遍形式定律，《水工結構工程理論與應用》，大連海運學院出版社，1993年6月。

[11] 楊光華，土的數學本構理論的研究(綜述報告)，第二屆全國青年巖土力學與工程會議論文集，大連理工大學出版社，1995年9月。

[12] 楊光華. 21 世紀應建立巖土材料的本構理論[J].巖土工程學報. 1997(03)

[13] 楊光華，土體材料本構特性的數學分析，第六屆全國巖土力學數值分析與解析方法討論會論文集，1998年8月。

[14] 楊光華，介玉新，李廣信，黃文峰，土的多重勢面模型及驗證，巖土工程學報，1999年第5期。

[15] 楊光華，李廣信，巖土本構模型的數學基礎與廣義位勢理論，巖土力學，2002年第5期。

[16] 楊光華，李廣信，從廣義位勢理論的角度看土的本構理論的研究，巖土工程學報，2007年第4期。

[17] 楊光華，巖土材料不符合Drucker 公設的一個證明，巖土工程學報[J]，2010年第1期。

[18] 楊光華,溫勇,鐘志輝. 基于廣義位勢理論的類劍橋模型[J].巖土力學. 2013(06)

[19] 楊光華,姚捷,溫勇. 考慮擬彈性塑性變形的土體彈塑性本構模型[J].巖土工程學報. 2013(08)

[20] Guang-hua Yang,Yu-xin Jie and Guang-xin Li,A Mathematical  Approach to Establishing Constitutive Models for Geometerials[J], Journal of Applied Mathematics, Volume,2013,Article ID 739068,10 pages. http//dx.doi.org/10.1155/2013/739068.

[21] 楊光華，土的現代本構理論的發展與回顧展望，巖土工程學報，2018年第8期。

[22] 郭萬里，朱俊高，彭文明，粗粒土的剪脹方程及廣義塑性本構模型研究，巖土工程學報，2018年第6期。