2
$\begingroup$

Background, here are the equations that I am trying to solve: Full Equations

Where R, E1, E2, V1, V2, P are all user inputs. X/A goes from -2 to 2 and Z/A goes from 0 to -2. Below is the code that I have so far. I created a list of inputs. Then created two arrays for the x and z inputs. The last is where I am having trouble. I'm trying to create a code such that it will hold a value for X constant in SX, SZ, and TXZ and plug in all the values for Z. Then move to the next value for X and plug all the values in for all the Z. The end goal is to create a density plot that for SX, SZ, and TXZ. Thank you!

R = .1;
E1 = 200*10^9;
E2 = 550*10^9;
P=1000;
V1 = 0.3;
V2 = 0.3;
E = 1/(((1-(V1^2))/E1)+((1-(V2^2))/E2));
A = ((.75*P*R)/(1.61172*10^11))^(1/3);
X = Range[-2 A, 2 A, 0.01*3*A];
Z = Range[0,-2 A, 0.005*3*A];

ZZ = ConstantArray[Z[[Range[Length[Z]]]], Length[X]];
XX = ConstantArray[X[[Range[Length[X]]]],Length[Z]];

For[i=1,i=Length[XX],
For[j=1,j = Length[ZZ],
M = Sqrt(.5*(((A^2-i^2+j^2)^2+4*i^2*j^2)^(.5)+(A^2-i^2+j^2)))
N = Sqrt(.5*(((A^2-i^2+j^2)^2+4*i^2*j^2)^(.5)-(A^2-i^2+j^2)))
SX = (-P/A)*M*((1-((j^2+N^2)/(M^2+N^2)))-2*N)
SZ = (-P/A)*M*((1-((j^2+N^2)/(M^2+N^2))))
SY = V1*(SX+SZ)
TXZ = (-P/A)*N*((M^2-j^2)/(M^2+N^2)),
DensityPlot[SX/P,XX/A,ZZ/A] 
]
]
$\endgroup$
3
  • $\begingroup$ You don't need to discretize the formula yourself, check the document of DensityPlot carefully. Also, notice e.g. E and N are built-in symbol in Mathematica (See their color? They're black, rather than blue ), you can't use them as variable names. $\endgroup$
    – xzczd
    Oct 30 '18 at 12:37
  • $\begingroup$ What do you want to determine from these equations? $\endgroup$ Oct 30 '18 at 13:15
  • $\begingroup$ So I just want to graph the equations. Im not looking for any numerical results. The formulas will tell me the stress and shear below the surface of a material, and I'd like to see the counter plots. $\endgroup$
    – Alc
    Oct 30 '18 at 13:19
5
$\begingroup$

After correcting all errors and normalization to A, we have

R = .1;
    E1 = 200*10^9;
    E2 = 550*10^9;
    P = 1000;
    V1 = 0.3;
    V2 = 0.3;
    E3 = 1/(((1 - (V1^2))/E1) + ((1 - (V2^2))/E2));
    A = ((.75*P*R)/(1.61172*10^11))^(1/3);
    X = Range[-2 A, 2 A, 0.01*3*A];
    Z = Range[0, -2 A, 0.005*3*A];




m = Sqrt[.5*(((1 - i^2 + j^2)^2 + 4*i^2*j^2)^(.5) + (1 - i^2 + j^2))];
 n = Sqrt[.5*(((1 - i^2 + j^2)^2 + 4*i^2*j^2)^(.5) - (1 - i^2 + 
        j^2))]; 
SX = (-P)*m*((1 - ((j^2 + n^2)/(m^2 + n^2))) - 2*j); 
SZ = (-P)*m*((1 - ((j^2 + n^2)/(m^2 + n^2)))); SY = V1*(SX + SZ);
 TXZ = (-P)*n*((m^2 - j^2)/(m^2 + n^2)); {DensityPlot[
  SX/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(x\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 DensityPlot[SY/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(y\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 DensityPlot[SZ/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 DensityPlot[TXZ/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Tau]\), \(xz\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}]}
{ContourPlot[SX/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(x\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 ContourPlot[SY/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(y\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 ContourPlot[SZ/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}], 
 ContourPlot[TXZ/P, {i, -2, 2}, {j, -2, 0}, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Tau]\), \(xz\)]\)/P", 
  PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}]}

fig1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.