ABAQUS混凝土塑性损伤模型参数生成

admin

import tkinter
import openpyxl
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import matplotlib
matplotlib.use('agg')
stress_c = []       # 压应力数组
strain_c = []       # 压应变数组
dc = []             # 受压损伤数组
stress_t = []       # 拉应力数组
strain_t = []       # 拉应变数组
dt = []             # 受拉损伤数组
# 定义各非弹性应变步的数值
scp_list = [0.000, 0.0004, 0.0008, 0.0012, 0.0016, 0.002, 0.0024, 0.0028, 0.0036, 0.005, 0.01, 0.02]
stp_list = [0.000, 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0008, 0.0010, 0.0020, 0.0030, 0.0040, 0.0050]


def cal_con(f, er, var0, var1, var2, var3, var4, text1, text2, text3, text4, frame):
        try:
                fcuk = float(f.get())  # 获取输入的混凝土抗压强度
                epc = float(er.get())  # 获取输入的抗压弹性强度与极限强度的比值
        except:
                var0.set('数据输入有误！请重新输入。')
        fck = 0.88*0.76*1*fcuk
        cv = 0.67/(fcuk**0.45)      # 变异系数根据混凝土抗拉强度和立方体抗压强度的关系拟合得到,无理论依据...
        ac2 = min(1.0, 1-(fcuk-40)*0.13/40)     # C40以上的混凝土考虑脆性折减系数ac2
        ftk = 0.88*0.395*fcuk**0.55*(1-1.645*cv)**0.45*ac2
        e0 = 10**5/(2.2+34.7/fcuk)
        str_fck = str(fck).split('.')[0] + '.' + str(fck).split('.')[1][0:2]
        str_ftk = str(ftk).split('.')[0] + '.' + str(ftk).split('.')[1][0:2]
        str_e0 = str(e0).split('.')[0] + '.' + str(e0).split('.')[1][0:2]
        var1.set(str_fck)
        var2.set(str_ftk)
        var3.set(str_e0)
        s_cr = (700 + 172 * fcuk ** 0.5) / 10 ** 6
        ac = 0.157 * fcuk ** 0.785 - 0.905
        n = e0 * s_cr / (e0 * s_cr - fck)
        s0 = Cal_con_c0(s_cr, ac, n, fck, e0, epc).s0     # 调用求取混凝土受压曲线上，非弹性应变为零位置处应力值的类
        ec = fck * epc / s0  # 计算割线模量
        str_ec = str(ec).split('.')[0] + '.' + str(ec).split('.')[1][0:2]
        var4.set(str_ec)  # 将割线模量传递给窗口标签变量
        s_tr = ftk ** 0.54 * 65 / 10 ** 6
        at = 0.312 * ftk ** 2
        global stress_c
        global strain_c
        global dc
        global stress_t
        global strain_t
        global dt
        stress_c = [0]  # 压应力数组
        strain_c = [0]  # 压应变数组
        dc = [0]  # 受压损伤数组
        stress_t = [0]  # 拉应力数组
        strain_t = [0]  # 拉应变数组
        dt = [0]  # 受拉损伤数组
        for scp in scp_list:
                tmp = cal_con_c(s_cr, ac, n, fck, e0, ec, epc, scp)
                stress_c.append(tmp[0])
                strain_c.append(tmp[1])
                dc.append(tmp[2])
        for stp in stp_list:
                tmp1 = cal_con_t(ftk, s_tr, at, e0, ec, stp)
                stress_t.append(tmp1[0])
                strain_t.append(tmp1[1])
                dt.append(tmp1[2])
        for i in range(len(strain_c)-1):      # 将混凝土受压的应变、应力、损伤数据加入到图表数组中
                text1.insert('end', str(stress_c[i+1])[:5] + '\t' + str(scp_list[i])+'\n')     # 应力-非弹性应变输入文本框
                text2.insert('end', str(dc[i+1])[:5] + '\t' + str(scp_list[i]) + '\n')         # 损伤-非弹性应变输入文本框
        for i in range(len(strain_t)-1):      # 将混凝土受拉的应变、应力、损伤数据加入到图表数组中
                text3.insert('end', str(stress_t[i+1])[:5] + '\t' + str(stp_list[i]) + '\n')  # 应力-非弹性应变输入文本框
                text4.insert('end', str(dt[i+1])[:5] + '\t' + str(stp_list[i]) + '\n')  # 损伤-非弹性应变输入文本框
        plt.clf()
        figure = plt.figure(figsize=(4.9, 5.5))
        ax = []
        for i in range(4):
                ax.append(figure.add_subplot(2, 2, i + 1))
        figure.subplots_adjust(left=0.13, right=0.95, bottom=0.15, top=0.95, wspace=0.4, hspace=0.35)                # 调整绘图布局
        ax[0].plot(strain_c, stress_c)
        ax[0].set_xlabel('compression strain')
        ax[0].set_ylabel('compression stress')
        ax[1].set_xlabel('tension strain')
        ax[1].set_ylabel('tension stress')
        ax[2].set_xlabel('compression strain')
        ax[2].set_ylabel('compression damage')
        ax[3].set_xlabel('tension strain')
        ax[3].set_ylabel('tension damage')
        ax[2].plot(strain_c, dc)
        ax[1].plot(strain_t, stress_t)
        ax[3].plot(strain_t, dt)
        cvs1 = FigureCanvasTkAgg(figure, master=frame)
        cvs1.get_tk_widget().place(width=400, height=400)
        var0.set('计算完成!')


class Cal_con_c0:    # 求解混凝土受压曲线上，非弹性应变为零位置的应力值
        def __init__(self, s_cr, ac, n, fck, e0, epc):
                self.s_cr = s_cr
                self.ac = ac
                self.n = n
                self.fck = fck
                self.e0 = e0
                self.epc = epc
                pc = fck / (e0 * s_cr)
                s0 = s_cr / 2
                res = 1
                m = 0
                s_min = 0
                s_max = s_cr
                while (abs(res) > 10 ** (-8)) and m < 100:  # 求取极限弹性应力所对应的应变，即非弹性应变等于零的位置。
                        res = pc * n * e0 * s0 / (n - 1 + (s0 / s_cr) ** n) - fck * epc
                        if res > 0:
                                s_max = s0
                                s0 = (s0 + s_min) / 2
                        else:
                                s_min = s0
                                s0 = (s0 + s_max) / 2
                        m += 1
                self.s0 = s0


def cal_con_c(s_cr, ac, n, fck, e0, ec, epc, scp):    # 混凝土受压时应力应变曲线定义/求解各非弹性应变对应的应力、应变、损伤
        sin_cr = s_cr - fck / ec  # 计算抗压强度峰值位置的非弹性应变，用于后续迭代求解各非弹性应变值位置时，判别迭代区间。
        pc = fck / (e0 * s_cr)
        s0 = Cal_con_c0(s_cr, ac, n, fck, e0, epc).s0
        if scp < 10**(-8):
                dc = 0.0
                return fck*epc, s0, dc
        elif scp <= sin_cr:     # 求取上升段的应力、应变及损伤参数。
                res = 1
                m = 0
                s_min = s0
                s_max = s_cr
                s = s_cr / 2
                while (abs(res) > 10 ** (-8)) and m < 100:
                        res = s - pc * n * e0 * s / (n - 1 + (s / s_cr) ** n)/ec - scp
                        if res > 0:
                                s_max = s
                                s = (s + s_min) / 2
                        else:
                                s_min = s
                                s = (s + s_max) / 2
                        m += 1
                stress_c = pc * n * e0 * s / (n - 1 + (s / s_cr) ** n)
                dc = 1 - (stress_c/s/ec)**0.5
                return stress_c, s, dc
        else:       # 求取下降段的应力、应变及损伤参数。
                res = 1
                m = 0
                s_min = s_cr
                s_max = 5
                s = s_cr*2
                while (abs(res) > 10 ** (-8)) and m < 100:
                        res = s - pc * e0 * s / (ac*(s/s_cr - 1)**2 + s/s_cr)/ec - scp
                        if res > 0:
                                s_max = s
                                s = (s + s_min) / 2
                        else:
                                s_min = s
                                s = (s + s_max) / 2
                        m += 1
                stress_c = pc * e0 * s / (ac*(s/s_cr - 1)**2 + s/s_cr)
                dc = 1 - (stress_c/s/ec)**0.5
                return stress_c, s, dc


def cal_con_t(ftk, s_tr, at, e0, ec, stp):    # 混凝土受拉时应力应变曲线定义/求解各非弹性应变对应的应力、应变、损伤
        if stp < 10**(-10):
                dt = 0.0
                return ftk, s_tr, dt
        else:
                res = 1
                m = 0
                s_min = s_tr
                s_max = 10
                s = s_tr * 2
                pt = ftk/(e0*s_tr)
                while (abs(res) > 10 ** (-8)) and m < 100:
                        res = s - pt * e0 * s / (at * (s / s_tr - 1) ** 1.7 + s / s_tr) / ec - stp
                        if res > 0:
                                s_max = s
                                s = (s + s_min) / 2
                        else:
                                s_min = s
                                s = (s + s_max) / 2
                        m += 1
                stress_t = pt * e0 * s / (at * (s / s_tr - 1) ** 1.7 + s / s_tr)
                dt = 1 - (stress_t / s / ec)**0.5
                return stress_t, s, dt


def write_out():
        wb = openpyxl.Workbook()
        ws = wb.active
        heard = ['混凝土塑性损伤本构(混规)']
        tt = ['压应力(Mpa)', '非弹性压应变', '受压损伤', '非弹性压应变', '拉应力(Mpa)', '非弹性拉应变', '受拉损伤', '非弹性拉应变']
        ws.append(heard)
        ws.append(tt)
        for n in range(len(stress_c) - 1):
                ws.cell(n + 3, 1).value = stress_c[n+1]
        for n in range(len(stress_c) - 1):
                ws.cell(n + 3, 2).value = scp_list[n]
                ws.cell(n + 3, 4).value = scp_list[n]
        for n in range(len(stress_c) - 1):
                ws.cell(n + 3, 3).value = dc[n+1]
        for n in range(len(stress_t) - 1):
                ws.cell(n + 3, 5).value = stress_t[n+1]
        for n in range(len(stress_t) - 1):
                ws.cell(n + 3, 6).value = stp_list[n]
                ws.cell(n + 3, 8).value = stp_list[n]
        for n in range(len(stress_t) - 1):
                ws.cell(n + 3, 7).value = dt[n+1]
        ws.cell(1, 1).font = openpyxl.styles.Font(name='黑体', size=16)
        ws.cell(1, 1).alignment = openpyxl.styles.Alignment(horizontal='center', vertical='center')
        ws.row_dimensions[1].height = 25
        ws.row_dimensions[2].height = 25
        for i in range(20):
                for m in range(20):
                        ws.cell(i + 1, m + 1).alignment = openpyxl.styles.Alignment(horizontal='center', vertical='center', wrapText=True)
        ws.merge_cells('A1:H1')
        wb.save('混凝土塑性损伤本构.xlsx')


class Application(tkinter.Frame):
        def __init__(self, master=None):
                super().__init__(master)
                self.master = master
                self.place(x=0, y=0, width=900, height=600)
                self.create_widgets()

        def create_widgets(self):
                tkinter.Label(self, text=' Fcu.k ').place(x=5, y=0, width=30, height=25)  # 输入混凝土立方体抗压强度标签
                f = tkinter.Entry(self, fg='red')
                f.place(x=1, y=25, width=38, height=24)
                tkinter.Label(self, text='Fc0/Fck').place(x=45, y=0, width=45, height=25)  # 输入弹性极限强度与峰值强度比值标签
                var5 = tkinter.StringVar()
                var5.set('0.40')  # 设定混凝土弹性极限强度与峰值强度的比值的默认值为0.4。
                er = tkinter.Entry(self, textvariable=var5, fg='green')
                er.place(x=48, y=25, width=38, height=24)
                tkinter.Label(self, text='Fck').place(x=95, y=0, width=50, height=25)  # 输出混凝土轴心抗压强度标准值标签
                var1 = tkinter.StringVar()  # 轴压强度显示变量
                var1.set('未计算')
                tkinter.Label(self, textvariable=var1, width=9, height=1, bg='LightGray').place(x=95, y=25, width=50, height=24)  # 输出轴压强度
                tkinter.Label(self, text='Ftk').place(x=150, y=0, width=50, height=25)  # 输出混凝土轴心抗拉强度标准值标签
                var2 = tkinter.StringVar()  # 轴拉强度显示变量
                var2.set('未计算')
                tkinter.Label(self, textvariable=var2, width=9, height=1, bg='LightGray').place(x=150, y=25, width=50, height=24)  # 输出轴拉强度
                tkinter.Label(self, text='E0').place(x=205, y=0, width=50, height=25)
                var3 = tkinter.StringVar()  # 弹性模量显示变量
                var3.set('未计算')
                tkinter.Label(self, textvariable=var3, width=9, height=1, bg='LightGray').place(x=205, y=25, width=55, height=24)
                tkinter.Label(self, text='Ec').place(x=265, y=0, width=50, height=25)
                var4 = tkinter.StringVar()  # 弹性模量(割线模量)显示变量
                var4.set('未计算')
                tkinter.Label(self, textvariable=var4, width=9, height=1, bg='LightGray').place(x=265, y=25, width=55, height=24)
                frame = tkinter.Frame(self, width=400, height=392)                # 可通过加relief='ridge', borderwidth=1参数显示窗口轮廓
                frame.place(x=5, y=80)
                canv1 = tkinter.Canvas(frame)
                canv1.place(x=5, y=80)
                # 执行计算混凝土本构程序
                tkinter.Button(self, text='计算', command=lambda: cal_con(f, er, var0, var1, var2, var3, var4, text1, text2,        text3, text4, frame), height=1, width=4).place(x=5, y=50, width=35, height=25)
                tkinter.Button(self, text='导出excel', command=write_out).place(x=45, y=50, width=60, height=25)  # 数据导出excel
                var0 = tkinter.StringVar()
                var0.set('请输入混凝土立方体抗压强度值')
                text0 = tkinter.Label(self, textvariable=var0, font=('宋体', 10), bg='white', fg='blue', anchor='center')
                text0.place(x=108, y=50, width=212, height=25)
                tkinter.Label(self, text='压应力-非弹性应变', font=('宋体', 9)).place(x=407, y=80)
                tkinter.Label(self, text='受压损伤-非弹性应变', font=('宋体', 9)).place(x=524, y=80)
                tkinter.Label(self, text='拉应力-非弹性应变', font=('宋体', 9)).place(x=407, y=280)
                tkinter.Label(self, text='受拉损伤-非弹性应变', font=('宋体', 9)).place(x=524, y=280)
                text1 = tkinter.Text(self, width=15, height=12, bg='white')  # 定义受压应力-非弹性应变输出文本框位置
                text1.place(x=410, y=99)
                text2 = tkinter.Text(self, width=15, height=12, bg='white')  # 定义受压损伤-非弹性应变输出文本框位置
                text2.place(x=530, y=99)
                text3 = tkinter.Text(self, width=15, height=13, bg='white')  # 定义受拉应力-非弹性应变输出文本框位置
                text3.place(x=410, y=300)
                text4 = tkinter.Text(self, width=15, height=13, bg='white')  # 定义受拉损伤-非弹性应变输出文本框位置
                text4.place(x=530, y=300)
                tkinter.Label(self, text='工具说明', font=('', 9)).place(x=470, y=2)
                text5 = tkinter.Text(self, width=44, height=3, bg='yellow', spacing1=2)  # 定义工具说明文本框位置
                text5.place(x=326, y=24)
                text5.insert(1.0, '工具的应力-应变曲线参照《混规》附录C，输入框处Fc0为受压曲线直线段极限应力，\
                Ec为直线段弹性模量(用于输入ABAQUS)，应力输入单位:Mpa。')
                text5.configure(state='disabled')


root = tkinter.Tk()
root.geometry('643x475')
root.title('混凝土塑性损伤本构计算（混规）_V1.0')
root.resizable(width=False, height=False)
app = Application(root)
app.mainloop()