技术教程 编程教程 修正剑桥模型预测-用python3.4 正文

修正剑桥模型预测-用python3.4

2023-06-02,,

下面是预测结果：

 #!/usr/bin/env python

 # -*- coding:utf-8 -*-

 # __author__ = "blzhu"

 """

 python study

 Date：2017

 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——围压sigma3=常数

 根据ε1求其余的量

 """

 # from numpy import *

 import numpy as np

 import string

 import matplotlib.pyplot as plt

 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称

 from pylab import *

 mpl.rcParams['font.sans-serif'] = ['SimHei']

 # 基本参数赋值

 λ = 0.0853

 κ = 0.0188

 M = 1.36

 ν = 0.3

 e0 = 0.68

 # sigma3 = 196

 p = np.array(np.zeros((18, 1)))

 p[0] = 196.0

 q = np.array(np.zeros((18, 1)))

 q[0] = 0.0

 η = np.array(np.zeros((18, 1)))

 Dpp = np.array(np.zeros((18, 1)))

 Dpq = np.array(np.zeros((18, 1)))

 Dqp = np.array(np.zeros((18, 1)))

 Dqq = np.array(np.zeros((18, 1)))

 dε1 = np.array(np.zeros((18, 1)))

 dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]])

 dε1[4:] = 0.02

 dε3 = np.array(np.zeros((18, 1)))

 dσ1 = np.array(np.zeros((18, 1)))

 dσ3 = np.array(np.zeros((18, 1)))

 dp = np.array(np.zeros((18, 1)))

 dq = np.array(np.zeros((18, 1)))

 dεv = np.array(np.zeros((18, 1)))

 εv = np.array(np.zeros((18, 1)))

 εv[0] = 0.0

 dεd = np.array(np.zeros((18, 1)))

 εd = np.array(np.zeros((18, 1)))

 εd[0] = 0.0

 σ1 = np.array(np.zeros((18, 1)))

 σ1[0] = 196.0

 σ3 = np.array(np.zeros((18, 1)))

 σ3[0] = 196.0

 ##############################

 # 计算

 for i in range(0,17):

     dσ3[i + 1] = 0.0

     σ3[i] = 196

     σ3[i + 1] = 196 + dσ3[i + 1]

     # # 柔度矩阵元素

     # 先求Ck和Cp

     Ck = κ / (1 + e0)

     Cp = (λ - κ) / (1 + e0)

     # 再求柔度矩阵元素D

     η[i] = q[i] / p[i]

     Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2))

     Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2)

     Dqp[i + 1] = Dpq[i+1]

     Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4)

     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])

     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]) / (

     Dpp[i+1] + 3 * Dpq[i+1] + 3 * Dqp[i+1] + 9 * Dqq[i+1])

     σ1[i + 1] = σ1[i] + dσ1[i + 1]

     p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0

     q[i + 1] = σ1[i + 1] - σ3[i + 1]

     η[i + 1] = q[i + 1] / p[i + 1]

     dp[i + 1] = 1 / 3 * (dσ1[i + 1] + 2 * dσ3[i + 1])

     dq[i + 1] = dσ1[i + 1] - dσ3[i + 1]

     dεv[i + 1] = dε1[i + 1] + 2 * dε3[i + 1]

     εv[i + 1] = εv[i] + dεv[i + 1]

     dεd[i + 1] = 2 / 3 * (dε1[i + 1] - dε3[i + 1])

     εd[i + 1] = εd[i] + dεd[i + 1]

 # 数据输出

 print('p:')

 print(p)

 print('q:')

 print(q)

 print('εd:')

 print(εd)

 print('εv:')

 print(εv)

 print('η:')

 print(η)

 print('dp:')

 print(dp)

 print('dεd:')

 print(dεd)

 print('dεv:')

 print(dεv)

 # 绘图

 plt.figure(1)#创建图表1

 ax1=plt.subplot(111)

 # plt.plot(p, q, 'b*')

 plt.xlabel('p(kPa)')

 plt.ylabel('q(kPa)')

 plt.title(U'应力路径')

 plt.plot(εd,η,'r--')

 plt.plot(εd,εv,'r--')

 plt.show()

第二个：

 #!/usr/bin/env python

 # -*- coding:utf-8 -*-

 # __author__ = "blzhu"

 """

 python study

 Date：2017

 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——p=常数

 根据ε1求其余的量

 """

 # from numpy import *

 import numpy as np

 import string

 import matplotlib.pyplot as plt

 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称

 from pylab import *

 mpl.rcParams['font.sans-serif'] = ['SimHei']

 # 基本参数赋值

 λ = 0.0853

 κ = 0.0188

 M = 1.36

 ν = 0.3

 e0 = 0.68

 p = np.array(np.zeros((18, 1)))

 p[0] = 196.0

 q = np.array(np.zeros((18, 1)))

 q[0] = 0.0

 η = np.array(np.zeros((18, 1)))

 Dpp = np.array(np.zeros((18, 1)))

 Dpq = np.array(np.zeros((18, 1)))

 Dqp = np.array(np.zeros((18, 1)))

 Dqq = np.array(np.zeros((18, 1)))

 dε1 = np.array(np.zeros((18, 1)))

 dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]])

 dε1[4:] = 0.02

 dε3 = np.array(np.zeros((18, 1)))

 dσ1 = np.array(np.zeros((18, 1)))

 dσ3 = np.array(np.zeros((18, 1)))

 dp = np.array(np.zeros((18, 1)))

 dq = np.array(np.zeros((18, 1)))

 dεv = np.array(np.zeros((18, 1)))

 εv = np.array(np.zeros((18, 1)))

 εv[0] = 0.0

 dεd = np.array(np.zeros((18, 1)))

 εd = np.array(np.zeros((18, 1)))

 εd[0] = 0.0

 σ1 = np.array(np.zeros((18, 1)))

 σ1[0] = 196.0

 σ3 = np.array(np.zeros((18, 1)))

 σ3[0] = 196.0

 ##############################

 # 计算

 for i in range(0, 17):

     # # 柔度矩阵元素

     # 先求Ck和Cp

     Ck = κ / (1 + e0)

     Cp = (λ - κ) / (1 + e0)

     # 再求柔度矩阵元素D

     η[i] = q[i] / p[i]

     Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2))

     Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2)

     Dqp[i + 1] = Dpq[i + 1]

     Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4)

     dσ3[i + 1] = -p[i] * dε1[i + 1] / (Dpq[i + 1] + 3 * Dqq[i + 1])

     σ3[i + 1] = σ3[i] + dσ3[i + 1]

     dσ1[i + 1] = 2 * p[i] * dε1[i + 1] / (Dpq[i + 1] + 3 * Dqq[i + 1])

     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])

     σ1[i + 1] = σ1[i] + dσ1[i + 1]

     p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0

     q[i + 1] = σ1[i + 1] - σ3[i + 1]

     η[i + 1] = q[i + 1] / p[i + 1]

     dp[i + 1] = 1 / 3 * (dσ1[i + 1] + 2 * dσ3[i + 1])

     dq[i + 1] = dσ1[i + 1] - dσ3[i + 1]

     dεv[i + 1] = dε1[i + 1] + 2 * dε3[i + 1]

     εv[i + 1] = εv[i] + dεv[i + 1]

     dεd[i + 1] = 2 / 3 * (dε1[i + 1] - dε3[i + 1])

     εd[i + 1] = εd[i] + dεd[i + 1]

 # 数据输出

 print('p:')

 print(p)

 print('q:')

 print(q)

 print('εd:')

 print(εd)

 print('εv:')

 print(εv)

 print('η:')

 print(η)

 print('dp:')

 print(dp)

 print('dεd:')

 print(dεd)

 print('dεv:')

 print(dεv)

 # 绘图

 plt.figure(1)  # 创建图表1

 ax1 = plt.subplot(111)

 # plt.plot(p, q, 'b*')

 plt.xlabel('p(kPa)')

 plt.ylabel('q(kPa)')

 plt.title(U'应力路径')

 plt.plot(εd, η, 'r*')

 plt.plot(εd, εv, 'r*')

 plt.show()

第三个：注意增量步一定要合适，不能太大，因为应力路径太贴近屈服面了。

 #!/usr/bin/env python

 # -*- coding:utf-8 -*-

 # __author__ = "blzhu"

 """

 python study

 Date：2017

 《土的本构关系-罗汀》5.5.2节修正剑桥模型预测——3-不排水剪切试验

 根据ε1求其余的量

 """

 # from numpy import *

 import numpy as np

 import string

 import matplotlib.pyplot as plt

 # 字体的默认设置中并没有中文字体，所以我们只要手动添加中文字体的名称

 from pylab import *

 mpl.rcParams['font.sans-serif'] = ['SimHei']

 # 基本参数赋值

 λ = 0.0853

 κ = 0.0188

 M = 1.36

 ν = 0.3

 e0 = 0.68

 p = np.array(np.zeros((18, 1)))

 p[0] = 196.0

 q = np.array(np.zeros((18, 1)))

 q[0] = 0.0

 η = np.array(np.zeros((18, 1)))

 Dpp = np.array(np.zeros((18, 1)))

 Dpq = np.array(np.zeros((18, 1)))

 Dqp = np.array(np.zeros((18, 1)))

 Dqq = np.array(np.zeros((18, 1)))

 dε1 = np.array(np.zeros((18, 1)))

 # dε1[1:4] = np.array([[0.005, ], [0.010, ], [0.015]])

 # dε1[4:] = 0.02

 dε1[:] = 0.003

 dε3 = np.array(np.zeros((18, 1)))

 dσ1 = np.array(np.zeros((18, 1)))

 dσ3 = np.array(np.zeros((18, 1)))

 dp = np.array(np.zeros((18, 1)))

 dq = np.array(np.zeros((18, 1)))

 dεv = np.array(np.zeros((18, 1)))

 εv = np.array(np.zeros((18, 1)))

 εv[0] = 0.0

 dεd = np.array(np.zeros((18, 1)))

 εd = np.array(np.zeros((18, 1)))

 εd[0] = 0.0

 σ1 = np.array(np.zeros((18, 1)))

 σ1[0] = 196.0

 σ3 = np.array(np.zeros((18, 1)))

 σ3[0] = 196.0

 ##############################

 # 计算

 for i in range(0, 17):

     # # 柔度矩阵元素

     # 先求Ck和Cp

     Ck = κ / (1 + e0)

     Cp = (λ - κ) / (1 + e0)

     # 再求柔度矩阵元素D

     η[i] = q[i] / p[i]

     Dpp[i + 1] = Ck + Cp * ((M ** 2 - η[i] ** 2) / (M ** 2 + η[i] ** 2))

     Dpq[i + 1] = Cp * (2 * η[i]) / (M ** 2 + η[i] ** 2)

     Dqp[i + 1] = Dpq[i + 1]

     Dqq[i + 1] = (2 / 9) * Ck * (1 + ν) / (1 - 2 * ν) + Cp * (4 * η[i] ** 2) / (M ** 4 - η[i] ** 4)

     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])

     σ3[i + 1] = σ3[i] + dσ3[i + 1]

     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])

     dε3[i + 1] = (-1 / 2.0) * dε1[i + 1]

     σ1[i + 1] = σ1[i] + dσ1[i + 1]

     p[i + 1] = (σ1[i + 1] + 2 * σ3[i + 1]) / 3.0

     q[i + 1] = σ1[i + 1] - σ3[i + 1]

     η[i + 1] = q[i + 1] / p[i + 1]

     dp[i + 1] = (1 / 3.0) * (dσ1[i + 1] + 2 * dσ3[i + 1])

     dq[i + 1] = dσ1[i + 1] - dσ3[i + 1]

     # 不排水剪切路径

     dεv[i + 1] = dε1[i + 1]+2*dε3[i + 1]

     εv[i + 1] = εv[i] + dεv[i + 1]

     dεd[i + 1] = (2.0 / 3.0) * (dε1[i + 1] - dε3[i + 1])

     εd[i + 1] = εd[i] + dεd[i + 1]

 # 数据输出

 print('p:')

 print(p)

 print('q:')

 print(q)

 print('εd:')

 print(εd)

 print('εv:')

 print(εv)

 print('η:')

 print(η)

 print('dp:')

 print(dp)

 print('dεd:')

 print(dεd)

 print('dεv:')

 print(dεv)

 # 绘图

 plt.figure(1)  # 创建图表1

 ax1 = plt.subplot(111)

 # plt.plot(p, q, 'b*')

 plt.xlabel('p(kPa)')

 plt.ylabel('q(kPa)')

 plt.title(U'应力应变关系')

 # plt.plot(p,q)

 plt.plot(εd, η, 'r*')

 plt.plot(εd, εv, 'b*')

 plt.show()

用的python3.4+pycharm编译器，这个编译器可以按列选择，上面的代码可以输出数组，按列选择可以方便的放入excel中，之后处理。

下面是输出的数据在excel中：

 围压不变                    p不变                    不排水

 p    q    εd    q/p    εv    p    q    εd    q/p    εv    p    q    εd    q/p    εv

 196    0    0    0    0    196    0    0    0    0    196    0    0    0    0

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

 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

修正剑桥模型预测-用python3.4的相关教程结束。

《修正剑桥模型预测-用python3.4.doc》

下载本文的Word格式文档，以方便收藏与打印。

Copyright 2023. 昆居客 - IT编程,互联网信息技术文档大全!