# How to solve the wave velocity in steel

In the simulation of stress wave propagation, I have the following two problems.

First question:

This question comes from page 69 of this book.

The Lame equation for a linear elastic body without volume force is as follows:

Grad[Div[u[x, t], {x, t}], {x, t}] + μ*
Laplacian[u[x, t], {x, t}] == ρ*D[u[x, t], {t, 2}]


The relationship between the parameters λ and μ, Young's modulus E and Poisson's ratio ν in the above figure are shown as follows:

If the material is steel, then $$E = 2.10*10^5MPa$$, $$ν = 0.28$$, $$ρ = 0.01 g/mm^3$$. Assuming that the stress wave is an irrotational wave, how to calculate the velocity of the wave propagating in the steel.

Second question:

I saw the propagation animation of Rayleigh wave here, but the post only provided CDF file at the end of the post. I want to make this animation with MMA. How can I get the same effect as in the animation.

• There should be nothing that prevents you from tearing a CDF apart to use the underlying code, at least from my understanding/limited experience. Jul 24, 2020 at 7:01
• @CATrevillian Thank you very much for your reply, but I don't know how to extract the relevant code from the CDF file. Could you give me a post about how to extract the code? Jul 24, 2020 at 7:08
• I’ll see what I can muster :) it may be several hours from now, but I’ll try to show you what I mean Jul 24, 2020 at 7:16
• @CATrevillian Thank you very much for your help. Jul 24, 2020 at 7:17
• @Ordinaryusers68 Are looking for some code generating numerical solution of the elastic wave propagation like this one physics.stackexchange.com/questions/524928/… ? Jul 24, 2020 at 10:52

To visualize 3D elastic P-,R-,S-wave we use standard 2D FEM solver (see tutorial) and ListPointPlot3D[]. This code generates S-wave:

 Needs["NDSolveFEM"]; \[CapitalOmega] = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{0, 5}, {0, 1}}, "MaxCellMeasure" -> 0.03];
mesh["Wireframe"];
diffusionCoefficients = "DiffusionCoefficients" -> {{{{-(Y/(1 - \[Nu]^2)), 0}, {0, -((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)))}}, {{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))), 0}}}, {{{0, -((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)))}, {-((Y*\[Nu])/(1 - \[Nu]^2)), 0}},
{{-((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))), 0}, {0, -(Y/(1 - \[Nu]^2))}}}} /. {Y -> 10^2, \[Nu] -> 33/100}; massCoefficients = "MassCoefficients" -> {{1, 0}, {0, 1}};
vd = NDSolveVariableData[{"Time", "DependentVariables", "Space"} -> {t, {u, v}, {x, y}}];
sd = NDSolveSolutionData[{"Time", "Space"} -> {0., ToNumericalRegion[mesh]}];
methodData = InitializePDEMethodData[vd, sd];
Subscript[\[CapitalGamma], Nv] = NeumannValue[0., y == 1]; Subscript[\[CapitalGamma], Du] = DirichletCondition[u[x, y] == 0, y == 0]; Subscript[\[CapitalGamma], Nu] = NeumannValue[0, x == 5];
Subscript[\[CapitalGamma], Dv] = DirichletCondition[v[x, y] == 0.05*Sin[Pi*(x + t)], y == 0 || y == 1];

initCoeffs = InitializePDECoefficients[vd, sd, {diffusionCoefficients, massCoefficients, loadCoefficients}];
initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[\[CapitalGamma], Du], Subscript[\[CapitalGamma], Nu]}, {Subscript[\[CapitalGamma], Dv], Subscript[\[CapitalGamma], Nv]}}];
sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Stationary"];
rhs[t_?NumericQ, uv_, duv_] := Module[{l, s, d, m, tdpde, tbcs, rayleighDamping},
NDSolveSetSolutionDataComponent[sd, "Time", t];
{l, s, d, m} = sdpde["SystemMatrices"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, "Transient"];
{l, s, d, m} += tdpde["SystemMatrices"];
rayleighDamping = 0.1*m + 0.04*s;
DeployBoundaryConditions[{l, s, rayleighDamping, m}, tbcs];
DeployBoundaryConditions[{l, s, rayleighDamping, m}, sbcs];
l - s.uv - rayleighDamping.duv
]
dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];
mass = sdpde["MassMatrix"];
stiff = sdpde["StiffnessMatrix"];
rd = 0.1*mass + 0.04*stiff;
sparsity = ArrayFlatten[{{mass["PatternArray"], mass["PatternArray"]}, {rd["PatternArray"], rd["PatternArray"]}}];
Dynamic["time: " <> ToString[CForm[currentTime]]];

tfun = NDSolveValue[{
mass.uv''[ t] == rhs[t, uv[t], uv'[t]]
, uv[ 0] == init, uv'[ 0] == dinit}, uv, {t, 0, 8}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {Automatic, Sparse -> sparsity}
, EvaluationMonitor :> (currentTime = t;)
];
split = Span @@@ Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];


Visualization

Do[mesht = Function[t,
dmesh =
ElementMeshDeformation[mesh, Part[tfun[t], #] & /@ split,
"ScalingFactor" -> 2]] /@ {k};
c3D = Table[{dmesh["Coordinates"][[i, 1]], j/6,
dmesh["Coordinates"][[i, 2]]}, {i,
Length[dmesh["Coordinates"]]}, {j, 0, 6}];
frame3D[k] =
ListPointPlot3D[Flatten[c3D, 1], ColorFunction -> "Rainbow",
BoxRatios -> Automatic, Boxed -> False, Axes -> False];, {k, 4, 8,
1/10}];


To save picture as a gif file we use

Export["C:\\...\\waveS3D.gif",
Table[frame3D[k], {k, 4, 8, 1/10}], AnimationRepetitions -> Infinity]


• Thank you very much, but when I ran the second piece of code(Do[mesht = Function[...1/10}];), the program output some error messages (version 12.1). Jul 25, 2020 at 23:51
• @Ordinaryusers68 Sorry, it should be  "ScalingFactor" -> 2`. Code been updated. Jul 26, 2020 at 11:16