下面是预测结果：
1 #!/usr/bin/env python 2 # -*- coding:utf-8 -*- 3 # __author__ = "blzhu" 4 """ 5 python study 6 Date：2017 7 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——围压sigma3=常数 8 根据ε1求其余的量 9 """ 10 # from numpy import * 11 import numpy as np 12 import string 13 import matplotlib.pyplot as plt 14 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称 15 from pylab import * 16 mpl.rcParams['font.sans-serif'] = ['SimHei'] 17 18 # 基本参数赋值 19 λ = 0.0853 20 κ = 0.0188 21 M = 1.36 22 ν = 0.3 23 e0 = 0.68 24 # sigma3 = 196 25 p = np.array(np.zeros((18, 1))) 26 p[0] = 196.0 27 q = np.array(np.zeros((18, 1))) 28 q[0] = 0.0 29 η = np.array(np.zeros((18, 1))) 30 Dpp = np.array(np.zeros((18, 1))) 31 Dpq = np.array(np.zeros((18, 1))) 32 Dqp = np.array(np.zeros((18, 1))) 33 Dqq = np.array(np.zeros((18, 1))) 34 dε1 = np.array(np.zeros((18, 1))) 35 dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]]) 36 dε1[4:] = 0.02 37 dε3 = np.array(np.zeros((18, 1))) 38 dσ1 = np.array(np.zeros((18, 1))) 39 dσ3 = np.array(np.zeros((18, 1))) 40 dp = np.array(np.zeros((18, 1))) 41 dq = np.array(np.zeros((18, 1))) 42 dεv = np.array(np.zeros((18, 1))) 43 εv = np.array(np.zeros((18, 1))) 44 εv[0] = 0.0 45 dεd = np.array(np.zeros((18, 1))) 46 εd = np.array(np.zeros((18, 1))) 47 εd[0] = 0.0 48 σ1 = np.array(np.zeros((18, 1))) 49 σ1[0] = 196.0 50 σ3 = np.array(np.zeros((18, 1))) 51 σ3[0] = 196.0 52 ############################## 53 # 计算 54 for i in range(0,17): 55 dσ3[i + 1] = 0.0 56 σ3[i] = 196 57 σ3[i + 1] = 196 + dσ3[i + 1] 58 59 # # 柔度矩阵元素 60 # 先求Ck和Cp 61 Ck = κ / (1 + e0) 62 Cp = (λ - κ) / (1 + e0) 63 # 再求柔度矩阵元素D 64 η[i] = q[i] / p[i] 65 Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2)) 66 Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2) 67 Dqp[i + 1] = Dpq[i+1] 68 Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4) 69 dσ1[i + 1] = 9 * p[i] * dε1[i+1] / (Dpp[i+1] + 3 * Dpq[i+1] + 3 * Dqp[i+1] + 9 * Dqq[i+1]) 70 dε3[i + 1] = ((1 / 2.0) * dε1[i+1])* (2 * Dpp[i+1] + 6 * Dpq[i+1] - 3 * Dqp[i+1] - 9 * Dqq[i+1]) / ( 71 Dpp[i+1] + 3 * Dpq[i+1] + 3 * Dqp[i+1] + 9 * Dqq[i+1]) 72 σ1[i + 1] = σ1[i] + dσ1[i + 1] 73 p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0 74 q[i + 1] = σ1[i + 1] - σ3[i + 1] 75 η[i + 1] = q[i + 1] / p[i + 1] 76 dp[i + 1] = 1 / 3 * (dσ1[i + 1] + 2 * dσ3[i + 1]) 77 dq[i + 1] = dσ1[i + 1] - dσ3[i + 1] 78 dεv[i + 1] = dε1[i + 1] + 2 * dε3[i + 1] 79 εv[i + 1] = εv[i] + dεv[i + 1] 80 dεd[i + 1] = 2 / 3 * (dε1[i + 1] - dε3[i + 1]) 81 εd[i + 1] = εd[i] + dεd[i + 1] 82 # 数据输出 83 print('p:') 84 print(p) 85 print('q:') 86 print(q) 87 print('εd:') 88 print(εd) 89 print('εv:') 90 print(εv) 91 print('η:') 92 print(η) 93 print('dp:') 94 print(dp) 95 print('dεd:') 96 print(dεd) 97 print('dεv:') 98 print(dεv) 99 # 绘图 100 plt.figure(1)#创建图表1 101 ax1=plt.subplot(111) 102 # plt.plot(p, q, 'b*') 103 plt.xlabel('p(kPa)') 104 plt.ylabel('q(kPa)') 105 plt.title(U'应力路径') 106 plt.plot(εd,η,'r--') 107 plt.plot(εd,εv,'r--') 108 plt.show()
第二个：
1 #!/usr/bin/env python 2 # -*- coding:utf-8 -*- 3 # __author__ = "blzhu" 4 """ 5 python study 6 Date：2017 7 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——p=常数 8 根据ε1求其余的量 9 """ 10 # from numpy import * 11 import numpy as np 12 import string 13 import matplotlib.pyplot as plt 14 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称 15 from pylab import * 16 17 mpl.rcParams['font.sans-serif'] = ['SimHei'] 18 19 # 基本参数赋值 20 λ = 0.0853 21 κ = 0.0188 22 M = 1.36 23 ν = 0.3 24 e0 = 0.68 25 p = np.array(np.zeros((18, 1))) 26 p[0] = 196.0 27 q = np.array(np.zeros((18, 1))) 28 q[0] = 0.0 29 η = np.array(np.zeros((18, 1))) 30 Dpp = np.array(np.zeros((18, 1))) 31 Dpq = np.array(np.zeros((18, 1))) 32 Dqp = np.array(np.zeros((18, 1))) 33 Dqq = np.array(np.zeros((18, 1))) 34 dε1 = np.array(np.zeros((18, 1))) 35 dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]]) 36 dε1[4:] = 0.02 37 dε3 = np.array(np.zeros((18, 1))) 38 dσ1 = np.array(np.zeros((18, 1))) 39 dσ3 = np.array(np.zeros((18, 1))) 40 dp = np.array(np.zeros((18, 1))) 41 dq = np.array(np.zeros((18, 1))) 42 dεv = np.array(np.zeros((18, 1))) 43 εv = np.array(np.zeros((18, 1))) 44 εv[0] = 0.0 45 dεd = np.array(np.zeros((18, 1))) 46 εd = np.array(np.zeros((18, 1))) 47 εd[0] = 0.0 48 σ1 = np.array(np.zeros((18, 1))) 49 σ1[0] = 196.0 50 σ3 = np.array(np.zeros((18, 1))) 51 σ3[0] = 196.0 52 ############################## 53 # 计算 54 for i in range(0, 17): 55 # # 柔度矩阵元素 56 # 先求Ck和Cp 57 Ck = κ / (1 + e0) 58 Cp = (λ - κ) / (1 + e0) 59 # 再求柔度矩阵元素D 60 η[i] = q[i] / p[i] 61 Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2)) 62 Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2) 63 Dqp[i + 1] = Dpq[i + 1] 64 Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4) 65 dσ3[i + 1] = -p[i] * dε1[i + 1] / (Dpq[i + 1] + 3 * Dqq[i + 1]) 66 σ3[i + 1] = σ3[i] + dσ3[i + 1] 67 dσ1[i + 1] = 2 * p[i] * dε1[i + 1] / (Dpq[i + 1] + 3 * Dqq[i + 1]) 68 dε3[i + 1] = ((1 / 2.0) * dε1[i + 1]) * (2 * Dpq[i + 1] - 3 * Dqq[i + 1]) / (Dpq[i + 1] + 3 * Dqq[i + 1]) 69 σ1[i + 1] = σ1[i] + dσ1[i + 1] 70 p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0 71 q[i + 1] = σ1[i + 1] - σ3[i + 1] 72 η[i + 1] = q[i + 1] / p[i + 1] 73 dp[i + 1] = 1 / 3 * (dσ1[i + 1] + 2 * dσ3[i + 1]) 74 dq[i + 1] = dσ1[i + 1] - dσ3[i + 1] 75 dεv[i + 1] = dε1[i + 1] + 2 * dε3[i + 1] 76 εv[i + 1] = εv[i] + dεv[i + 1] 77 dεd[i + 1] = 2 / 3 * (dε1[i + 1] - dε3[i + 1]) 78 εd[i + 1] = εd[i] + dεd[i + 1] 79 # 数据输出 80 print('p:') 81 print(p) 82 print('q:') 83 print(q) 84 print('εd:') 85 print(εd) 86 print('εv:') 87 print(εv) 88 print('η:') 89 print(η) 90 print('dp:') 91 print(dp) 92 print('dεd:') 93 print(dεd) 94 print('dεv:') 95 print(dεv) 96 # 绘图 97 plt.figure(1) # 创建图表1 98 ax1 = plt.subplot(111) 99 # plt.plot(p, q, 'b*') 100 plt.xlabel('p(kPa)') 101 plt.ylabel('q(kPa)') 102 plt.title(U'应力路径') 103 plt.plot(εd, η, 'r*') 104 plt.plot(εd, εv, 'r*') 105 plt.show()
第三个：注意增量步一定要合适，不能太大，因为应力路径太贴近屈服面了。
1 #!/usr/bin/env python 2 # -*- coding:utf-8 -*- 3 # __author__ = "blzhu" 4 """ 5 python study 6 Date：2017 7 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——3-不排水剪切试验 8 根据ε1求其余的量 9 """ 10 # from numpy import * 11 import numpy as np 12 import string 13 import matplotlib.pyplot as plt 14 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称 15 from pylab import * 16 17 mpl.rcParams['font.sans-serif'] = ['SimHei'] 18 19 # 基本参数赋值 20 λ = 0.0853 21 κ = 0.0188 22 M = 1.36 23 ν = 0.3 24 e0 = 0.68 25 p = np.array(np.zeros((18, 1))) 26 p[0] = 196.0 27 q = np.array(np.zeros((18, 1))) 28 q[0] = 0.0 29 η = np.array(np.zeros((18, 1))) 30 Dpp = np.array(np.zeros((18, 1))) 31 Dpq = np.array(np.zeros((18, 1))) 32 Dqp = np.array(np.zeros((18, 1))) 33 Dqq = np.array(np.zeros((18, 1))) 34 dε1 = np.array(np.zeros((18, 1))) 35 # dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]]) 36 # dε1[4:] = 0.02 37 dε1[:] = 0.003 38 dε3 = np.array(np.zeros((18, 1))) 39 dσ1 = np.array(np.zeros((18, 1))) 40 dσ3 = np.array(np.zeros((18, 1))) 41 dp = np.array(np.zeros((18, 1))) 42 dq = np.array(np.zeros((18, 1))) 43 dεv = np.array(np.zeros((18, 1))) 44 εv = np.array(np.zeros((18, 1))) 45 εv[0] = 0.0 46 dεd = np.array(np.zeros((18, 1))) 47 εd = np.array(np.zeros((18, 1))) 48 εd[0] = 0.0 49 σ1 = np.array(np.zeros((18, 1))) 50 σ1[0] = 196.0 51 σ3 = np.array(np.zeros((18, 1))) 52 σ3[0] = 196.0 53 ############################## 54 # 计算 55 for i in range(0, 17): 56 # # 柔度矩阵元素 57 # 先求Ck和Cp 58 Ck = κ / (1 + e0) 59 Cp = (λ - κ) / (1 + e0) 60 # 再求柔度矩阵元素D 61 η[i] = q[i] / p[i] 62 Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2)) 63 Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2) 64 Dqp[i + 1] = Dpq[i + 1] 65 Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4) 66 dσ3[i + 1] = ((p[i]/3) * dε1[i + 1]) *(Dpp[i + 1]+3*Dpq[i + 1])/(Dqp[i + 1]*Dpq[i + 1]-Dqq[i + 1]*Dpp[i + 1]) 67 σ3[i + 1] = σ3[i] + dσ3[i + 1] 68 dσ1[i + 1] = (-p[i]/3.0) * dε1[i + 1] *(2*Dpp[i + 1]-3*Dpq[i + 1])/(Dqp[i + 1]*Dpq[i + 1]-Dqq[i + 1]*Dpp[i + 1]) 69 dε3[i + 1] = (-1 / 2.0) * dε1[i + 1] 70 σ1[i + 1] = σ1[i] + dσ1[i + 1] 71 p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0 72 q[i + 1] = σ1[i + 1] - σ3[i + 1] 73 η[i + 1] = q[i + 1] / p[i + 1] 74 dp[i + 1] = (1 / 3.0) * (dσ1[i + 1] + 2 * dσ3[i + 1]) 75 dq[i + 1] = dσ1[i + 1] - dσ3[i + 1] 76 # 不排水剪切路径 77 dεv[i + 1] = dε1[i + 1]+2*dε3[i + 1] 78 εv[i + 1] = εv[i] + dεv[i + 1] 79 dεd[i + 1] = (2.0 / 3.0) * (dε1[i + 1] - dε3[i + 1]) 80 εd[i + 1] = εd[i] + dεd[i + 1] 81 # 数据输出 82 print('p:') 83 print(p) 84 print('q:') 85 print(q) 86 print('εd:') 87 print(εd) 88 print('εv:') 89 print(εv) 90 print('η:') 91 print(η) 92 print('dp:') 93 print(dp) 94 print('dεd:') 95 print(dεd) 96 print('dεv:') 97 print(dεv) 98 # 绘图 99 plt.figure(1) # 创建图表1 100 ax1 = plt.subplot(111) 101 # plt.plot(p, q, 'b*') 102 plt.xlabel('p(kPa)') 103 plt.ylabel('q(kPa)') 104 plt.title(U'应力应变关系') 105 # plt.plot(p,q) 106 plt.plot(εd, η, 'r*') 107 plt.plot(εd, εv, 'b*') 108 plt.show()
用的python3.4+pycharm编译器，这个编译器可以按列选择，上面的代码可以输出数组，按列选择可以方便的放入excel中，之后处理。
下面是输出的数据在excel中：
1 围压不变 p不变 不排水 2 p q εd q/p εv p q εd q/p εv p q εd q/p εv 3 196 0 0 0 0 196 0 0 0 0 196 0 0 0 0 4 219.8033735 71.41012048 0.00294458 0.32488182 0.00616627 196 121.2569558 0.005 0.61865794 0 196 72.75417349 0.003 0.37119476 0 5 247.0206226 153.0618677 0.00939637 0.61963194 0.01681089 196 179.0710227 0.01284281 0.91362767 0.00647158 176.1387111 133.5818822 0.006 0.75839026 0 6 274.2116436 234.6349307 0.02061517 0.85567092 0.0281545 196 223.903844 0.02578835 1.14236655 0.01263495 153.9005543 162.3173732 0.009 1.05468998 0 7 300.7422139 314.2266416 0.03716315 1.04483716 0.03851054 196 252.3282794 0.04440249 1.28738918 0.01679252 141.7669587 171.3747924 0.012 1.20884862 0 8 319.5248863 370.5746587 0.05496147 1.15976775 0.04511559 196 261.9590164 0.0639265 1.33652559 0.0182205 136.0278412 174.5184164 0.015 1.28296101 0 9 332.7304557 410.191367 0.07353339 1.23280379 0.04939982 196 265.1033676 0.08377088 1.3525682 0.01868736 133.2091415 175.827666 0.018 1.31993694 0 10 341.7615351 437.2846054 0.0926109 1.27950211 0.0521673 196 266.1025194 0.10372143 1.35766592 0.01883572 131.7746429 176.4377929 0.021 1.33893585 0 11 347.7755554 455.3266661 0.11201989 1.30925437 0.05394033 196 266.416703 0.12370587 1.35926889 0.01888238 131.0292971 176.7402229 0.024 1.34886035 0 12 351.6977353 467.0932058 0.13164416 1.32810979 0.05506753 196 266.5151532 0.143701 1.35977119 0.018897 130.6376305 176.8951974 0.027 1.35409068 0 13 354.217724 474.6531719 0.1514067 1.34000401 0.0557799 196 266.5459683 0.16369947 1.35992841 0.01890158 130.4305684 176.976036 0.03 1.35686012 0 14 355.8204297 479.4612891 0.17125726 1.34748106 0.05622822 196 266.5556101 0.183699 1.3599776 0.01890301 130.3207472 177.0186055 0.033 1.35833019 0 15 356.8329394 482.4988181 0.19116348 1.35217006 0.05650956 196 266.5586267 0.20369885 1.35999299 0.01890346 130.2624003 177.0411363 0.036 1.35911158 0 16 357.4698304 484.4094911 0.21110474 1.35510594 0.05668578 196 266.5595704 0.2236988 1.35999781 0.0189036 130.2313729 177.0530933 0.039 1.3595272 0 17 357.8693447 485.608034 0.23106799 1.35694225 0.05679603 196 266.5598656 0.24369879 1.35999931 0.01890364 130.2148652 177.059448 0.042 1.35974835 0 18 358.119518 486.3585541 0.25104501 1.35809005 0.05686496 196 266.559958 0.26369878 1.35999979 0.01890365 130.2060802 177.0628279 0.045 1.35986605 0 19 358.2760029 486.8280086 0.27103066 1.35880719 0.05690803 196 266.5599869 0.28369878 1.35999993 0.01890366 130.2014044 177.0646263 0.048 1.3599287 0 20 358.3738175 487.1214525 0.29102169 1.35925514 0.05693493 196 266.5599959 0.30369878 1.35999998 0.01890366 130.1989156 177.0655834 0.051 1.35996204 0