Трехмерные поверхности выхода и разрушения в matplotlib: поверхность разрушения Мацуока-Накай - PullRequest
Новая
       6

Трехмерные поверхности выхода и разрушения в matplotlib: поверхность разрушения Мацуока-Накай

0 голосов
/ 09 сентября 2018

Я хочу построить трехмерную поверхность выхода и отказа в matplotlib. Пример типа поверхностей, которые я хотел бы построить, можно найти в https://en.wikipedia.org/wiki/Von_Mises_yield_criterion и https://en.wikipedia.org/wiki/Yield_surface.. В Maple я использовал implicitplot3d для генерации таких фигур. В среде Python можно использовать Mayavi.

В качестве примера я смог сгенерировать трехмерную поверхность сбоев для критерия Мацуока-Накаи из уравнений, представленных здесь http://rspa.royalsocietypublishing.org/content/472/2185/20150713, и преобразовав представленный здесь код MATLAB Антона Попова https://www.researchgate.net/post/How_can_I_plot_Mohr-Coulomb_yield_criterion_in_3D_using_matlab в эквивалентный код Python , Я только что искал эквивалентные функции в Python для функций MATLAB в коде Антона, поэтому я уверен, что моя версия может быть улучшена. Я не слишком знаком с созданием 3D-графиков, и мне было интересно, есть ли лучший способ создать тот же график. Также можно найти прикрепленное здесь изображение 3D-график: Поверхность разрушения Мацуока-Накаи . 

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt


def main():
    # plotting parameters
    num_polartheta = 12
    num_p = 6
    num_polarphi = 6 * num_polartheta - 5

    # range of hydrostatic stress (p) to plot p = (x + y + z)/3
    pmin = 0
    pmax = 100
    p = np.linspace(pmin, pmax, num=num_p)

    # 360 degrees in the shperical phi
    polar_phi = np.linspace(0, 2 * np.pi, num_polarphi)

    # along theta angle
    th1 = np.linspace(-np.pi / 6, np.pi / 6, num=num_polartheta)
    th2 = np.linspace(np.pi / 6, -np.pi / 6, num=num_polartheta)
    polar_theta = np.r_[th1, th2[1:], th1[1:], th2[1:], th1[1:], th2[1:]]

    # create grid
    grd_phi = np.tile(polar_phi, (num_p, 1))
    grd_th = np.tile(polar_theta, (num_p, 1))
    grd_p = np.tile(p.reshape(num_p, 1), (1, num_polarphi))

    # Material parameter
    # effective friction angle,
    mat_phi = 30
    mat_phi = np.radians(mat_phi)

    # Matsuoka-Nakai yield surface
    kmn = (9 - np.sin(mat_phi) ** 2) / (1 - np.sin(mat_phi) ** 2)
    A1 = (kmn - 3) / (kmn - 9)
    A2 = kmn / (kmn - 9)
    sqrtJ2 = grd_p / (
        2
        / np.sqrt(3)
        * np.sqrt(A1)
        * np.cos(1 / 3 * np.arccos(A2 / A1 ** (3 / 2) * np.sin(3 * grd_th)))
       )
    X = grd_p + np.sqrt(2 / 3) * sqrtJ2 * np.cos(grd_phi + 2 * np.pi / 3)
    Y = grd_p + np.sqrt(2 / 3) * sqrtJ2 * np.cos(grd_phi)
    Z = grd_p + np.sqrt(2 / 3) * sqrtJ2 * np.cos(grd_phi - 2 * np.pi / 3)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
    ax.view_init(elev=40, azim=45)
    ax.set_xlabel(r"$\sigma_1$")
    ax.set_ylabel(r"$\sigma_2$")
    ax.set_zlabel(r"$\sigma_3$")
    plt.show()


if __name__ == "__main__":
    main()
...