I encountered similar graphs like the one shown in the attached figure in many papers and researches. Does anyone use a certain software to draw this type of graphs?
Presuming you are talking about concrete model yield surface. I1, J2 and J are stress invariants and
sigv=I1/3; rho=sqrt(2*J2); theta=1/3*acos(J);
are principal stress coordinates in High-Westergard coordinate system. Value of yield surface for (sigv, rho, theta) is labeled with f, and positive values correspond to stress inside the surface, zero values on the surface and negative values outside surface. Tension strength is marked with ft and compression strength with fc.
1. Determine is yield surface has start point on hydrostatic axis tension (positive) side. Concrete yield surfaces should have starting point. Procedure would be to:
1. Start from stress point (sigv, rho, theta) =(0,0,0)
2. Determine f for that point (should be positive = inside surface)
3. Increase sigv for small increment e.g. 0.05·fc, leave rho and theta 0. This way you are moving on hydrostatic axis on positive (tension) side. Determine f and keep increasing sigv until f is negative. This means you are outside of yield surface.
4. Best way now is to use bisection method between current and previous increment from 3. to find point close enough to yield surface f>>0. This is starting point, store its coordinates.
2. Now start double for loop where outer iterates moving on hydrostatic axis towards negative side from starting point, and inner for loop iterates angle from 0° to 360° for current hydrostatic point (from outer loop).
1. Start first for loop of sigv from starting point from 1.4 towards i.e. -2.5·fc (should be enough) with the chosen number of steps i.e. 20.
2. Start second for loop from 0° to 360° corresponding to theta angle with chosen number of steps i.e. 24.
3. Increase rho for small increments of i.e. 0.25·fc, starting from 0. Determine f and keep increasing rho until f is negative. This means you are outside of yield surface.
4. Best way now is to use bisection method between current and previous increment from 3. to find point close enough to yield surface f>>0. Store this point, this is one vertex of your surface.
Finally plot or deviatoric sections from vertexes of outer for loop and plot all meridians that correspond to same theta angle.
Actually this type of graph can be plotted using differents software application such as Mathematica, MatLab or MicrocalOrigin among many others that have already recommended to you. For a certain researcher, the best software to solve a task like this can be Mathematica, while for another researcher it can be any other software. This selection depends on many factors that including, of course, the tastes and preferences of the person who does it. I prefer the MicrocalOrigin. I only can recommend you that if are solving a specific problem and one of your tasks is plotting a similar graph to this, i.e., a three-dimensional graph then, try it with the software that you have on hand and if you know how to work with this software much better, do not worry about whether it is better with one or the other software.