5
$\begingroup$

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:

enter image description here

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:

enter image description here

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.

enter image description here

$\endgroup$
6
  • 1
    $\begingroup$ There should be nothing that prevents you from tearing a CDF apart to use the underlying code, at least from my understanding/limited experience. $\endgroup$ Jul 24, 2020 at 7:01
  • 1
    $\begingroup$ @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? $\endgroup$ Jul 24, 2020 at 7:08
  • $\begingroup$ I’ll see what I can muster :) it may be several hours from now, but I’ll try to show you what I mean $\endgroup$ Jul 24, 2020 at 7:16
  • $\begingroup$ @CATrevillian Thank you very much for your help. $\endgroup$ Jul 24, 2020 at 7:17
  • 1
    $\begingroup$ @Ordinaryusers68 Are looking for some code generating numerical solution of the elastic wave propagation like this one physics.stackexchange.com/questions/524928/… ? $\endgroup$ Jul 24, 2020 at 10:52

1 Answer 1

10
$\begingroup$

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

 Needs["NDSolve`FEM`"]; \[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 = NDSolve`VariableData[{"Time", "DependentVariables", "Space"} -> {t, {u, v}, {x, y}}]; 
sd = NDSolve`SolutionData[{"Time", "Space"} -> {0., ToNumericalRegion[mesh]}];
methodData = InitializePDEMethodData[vd, sd];
loadCoefficients = "LoadCoefficients" -> {{0}, {0}};
    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},
  NDSolve`SetSolutionDataComponent[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]

Figure 1

$\endgroup$
2
  • $\begingroup$ 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). $\endgroup$ Jul 25, 2020 at 23:51
  • $\begingroup$ @Ordinaryusers68 Sorry, it should be ` "ScalingFactor" -> 2`. Code been updated. $\endgroup$ Jul 26, 2020 at 11:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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