9
$\begingroup$

I found the deflection of an orifice plate (circular plate with a hole) subject to uniform pressure using Mathematica's NDSolve functionality. The plate is fixed at its outer circumference and is free at its inner circumference (hole).

The orifice plate (domain: $\Omega$ for NDSolve) is shown (all dimensions are meter): enter image description here

The deflection of the plate follows closely the assumption of ``small strains''. The deflection can be easily found by solving the biharmonic equation (with only radial terms; assuming no change in $\theta$ component of deflection and neglccting $z$ variation because of small thickness) with "fixed" boundary conditions on the outer circumference and "free" boundary conditions (for shear force and moment in terms of deflection) on the inner hole circumference. In fact my solution matches tabulated values for deflections of.

My solution method is shown below

a = 10*10^-3;
b = 5*10^-3;
ν = 1/3;
p0 = 0.1*10^6;
Ey = 200 *10^9;
h = 1*10^-3;
De = (Ey h^3)/(12 (1 - ν^2));

sol = NDSolve[
  {
   w''''[r] + (2/r) w'''[r] - (1/(r^2)) w''[r] + (1/(r^3)) w'[
       r] == -p0/De,
   w[a] == 0,
   w'[a] == 0,
   -(Derivative[1][w][b]/b^2) + Derivative[2][w][b]/b + 
     Derivative[3][w][b] == 0,
   (ν Derivative[1][w][b])/b + Derivative[2][w][b] == 0
   },
  w, {r, b, a}]

Plotting is done

Plot[Evaluate[w[r] /. sol], {r, b, a}, 
 PlotRange -> {{0, a}, Automatic}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 18}, 
 AxesLabel -> {"Radial location [m]", "Transverse deflection [m]"}, 
 PlotStyle -> {Black, Thick}, ImageSize -> 800]

The radial deflection of the plate is plotted and matches very closely, tabulated data from the "Mechanical engineer's handbook by James Carvill" (data table is available on request if needed)

enter image description here

My question is, how do I assign boundary conditions in the form of DirichletCondition and NeumannValue etc.

I would really like to solve the biharmonic equation in $r$ on an axisymmetric domain of my choosing and plot the deflection. This would need proper application of boundary conditions for the domain boundaries that I have not been able to do. I would like to use a domain as specified by a region difference. the idea is to plot the deflections on the domain to show the utility of FEM to solve such problems. An example would be Stokes flow solved using FEM in Mathematica.

My attempt so far is to do the following:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
Needs["NDSolve`FEM`"]
a = 10*10^-3;
b = 2*10^-3;
ν = 1/3;
p0 = 0.1*10^6;
Ey = 200 *10^9;
h = 1*10^-3;
De = (Ey h^3)/(12 (1 - ν^2));
Ω = RegionDifference[Disk[{0, 0}, a], Disk[{0, 0}, b]];
RegionPlot[Ω, AspectRatio -> Automatic, 
 ImageSize -> 400, LabelStyle -> {24, GrayLevel[0]}]

Boundary conditions (How do I correctly obtain the second and third conditions, viz., $c_3, c_4$)

Subscript[Γ, d1] = DirichletCondition[0, r == a]
Subscript[Γ, n1] = NeumannValue[0, r == a]
Subscript[c, 3] = -(Derivative[1][w][b]/b^2) + (w^′′)[
   b]/b + 
\!\(\*SuperscriptBox[\(w\), 
TagBox[
RowBox[{"(", "3", ")"}],
Derivative],
MultilineFunction->None]\)[b]
Subscript[c, 4] = (ν Derivative[1][w][b])/
  b + (w^′′)[b]

Application of Boundary conditions in NDSolve for my domain $\Omega$ (UNABLE TO ACCOMPLISH THIS)

sol = NDSolve[
  {
   w''''[r] + (2/r) w'''[r] - (1/(r^2)) w''[r] + (1/(r^3)) w'[
       r] == -p0/De,
   Subscript[Γ, d1] == 0,
   Subscript[Γ, n1] == 0,
   -(Derivative[1][w][b]/b^2) + Derivative[2][w][b]/b + 
     Derivative[3][w][b] == 0,
   (ν Derivative[1][w][b])/b + Derivative[2][w][b] == 0
   },
  w, r ∈ Ω,
  Method -> {"FiniteElement", "InterpolationOrder" -> 2,
    "MeshOptions" -> {"MaxCellMeasure" -> 0.5, 
      "MaxBoundaryCellMeasure" -> 0.02}}
  ]

The FEM parameters chosen in NDSolve are not very well thought of currently and I am open to suggestions.

$\endgroup$
  • 1
    $\begingroup$ It seems (but I am not sure I understand what you are trying to achieve) your use of the DirichletCondition is not quite correct; it should be something like DirichletCondition[ w[r] == 0, r \[Element] boundary]? $\endgroup$ – chris Aug 8 '15 at 16:01
  • $\begingroup$ @chris Likely you are correct. I probably haven't used this correctly. However, I am also concerned about the other double and triple derivative conditions. Is there some way I can appoint these boundary conditions on the domain $\Omega$ without needing to resort to the (seemingly cryptic ;)') DirichletCondition` and NeumannValue types. $\endgroup$ – dearN Aug 8 '15 at 16:02
  • $\begingroup$ As an alternative you could use the plane stress Structural Mechanics example in 2D as a basis. $\endgroup$ – user21 Aug 10 '15 at 9:03
  • $\begingroup$ @user21 I am thinking this is a plane strain problem. And even then the boundary conditions may need to be beyond second order. Maybe in wrong about the latter and I will check. $\endgroup$ – dearN Aug 10 '15 at 10:32
  • $\begingroup$ There is also a plane strain example in that same section. If you have an application that needs higher order derivatives I'd be curious to know what it is. $\endgroup$ – user21 Aug 10 '15 at 10:42
9
$\begingroup$

Here is something to get you started: The idea is to use a 2D cross section of the orifice plate in the xz-direction (not the xy-direction, as then you could not apply the surface force).

a = 10*10^-3;
b = 5*10^-3;
ν = 1/3;
p0 = 0.1*10^6;
Ey = 200*10^9;
h = 1*10^-3;
De = (Ey h^3)/(12 (1 - ν^2));

planeStrain = {Inactive[Div][{{0, -((Ey*ν)/((1 - 2*ν)*(1 + ν)))}, 
           {-Ey/(2*(1 + ν)), 0}} . Inactive[Grad][v[x, y], 
           {x, y}], {x, y}] + Inactive[Div][
       {{-((Ey*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, 
           {0, -Ey/(2*(1 + ν))}} . Inactive[Grad][u[x, y], 
           {x, y}], {x, y}], 
   Inactive[Div][{{0, -Ey/(2*(1 + ν))}, 
           {-((Ey*ν)/((1 - 2*ν)*(1 + ν))), 0}} . 
         Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
     Inactive[Div][{{-Ey/(2*(1 + ν)), 0}, 
           {0, -((Ey*(1 - ν))/((1 - 2*ν)*(1 + ν)))}} . 
         Inactive[Grad][v[x, y], {x, y}], {x, y}]}/. {Y -> Ey};

r = Rectangle[{b, 0}, {a, h}];
Subscript[Γ, 
  D] = {DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == a]};
{uif, vif} = 
  NDSolveValue[{planeStrain == {0, NeumannValue[-p0/De, y == h]}, 
    Subscript[Γ, D]}, {u, v}, {x, y} ∈ r];
mesh = uif["ElementMesh"];
Show[{
  mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
  ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10^5][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

enter image description here

And the deformation at h/2 in the y-direction:

Plot[vif[x, h/2], {x, b, a}]

enter image description here

There is a scaling issue, I'll leave that to you to figure out. perhaps some factor is missing or the plane strain model is not quite right in this scenario. But in principal that's how you'd do it. Or, if you want to go and get the hammer you could do it in 3D with a stress operator (and that should give the correct answer with more thinking)

$\endgroup$
  • $\begingroup$ This is interesting. However, the XZ plane in your method is the plane of the screen (deducing from deflection direction)? $\endgroup$ – dearN Aug 12 '15 at 13:36
  • $\begingroup$ @drN, yes, the xy-plane is the plane of the orifice plate. The z direction is then out of the screen - the pressure is applied from that z direction onto the disk in xy. $\endgroup$ – user21 Aug 12 '15 at 13:44
  • 1
    $\begingroup$ I doubt that a 3D model would provide a different solution given the axisymmetric nature of the plate and the load applied. But this is definitely food for thought. It is interesting (and perhaps) to note that this FEM model in Mathematica surpasses in computation speed what is obtained through Hyperworks. Hyperworks doesn't allow the application of boundary conditions as derivatives and one is required to choose loose terms such as "free boundary" or "fixed boundary". My hope is to get mathematica significant recognition in my university's mechanical engineering dept. $\endgroup$ – dearN Aug 12 '15 at 13:52
  • $\begingroup$ Thank you for confirming my deduction of XY plane. $\endgroup$ – dearN Aug 12 '15 at 13:53
  • $\begingroup$ @drN, I very much appreciate that note about Hyperworks. Keeps me motivated to continue implementing the FEM in Mathematica ;-) $\endgroup$ – user21 Aug 12 '15 at 14:06
4
$\begingroup$

This works.

a = 10*10^-3;
b = 5*10^-3;
ν = 1/3;
p0 = 0.1*10^6;
Ey = 200 *10^9;
h = 1*10^-3;
De = (Ey h^3)/(12 (1 - ν^2));

eqn = w''''[r] + (2/r) w'''[r] - (1/(r^2)) w''[r] + (1/(r^3)) w'[
     r] == -p0/De

w1 = NDSolveValue[{eqn, DirichletCondition[w[r] == 0, r == a], 
   DirichletCondition[w'[r] == 0, r == a], DirichletCondition[
    -(Derivative[1][w][r]/r^2) + Derivative[2][w][r]/r + 
      Derivative[3][w][r] == 0, r == b], DirichletCondition[
    (ν Derivative[1][w][r])/r + Derivative[2][w][r] == 0, 
    r == b]}, w, {r, b, a}]


  Plot[w1[r], {r, b, a}]

Mathematica graphics

PS: this might be of relevance to you though..

$\endgroup$
  • $\begingroup$ No it doesn't: Encountered non-numerical value for a derivative at r == 0.005., but I suppose I'll amend my question: I would like to use a domain as specified by a region difference. the idea is to plot the deflections on the domain to show the utility of FEM to solve such problems. I hope that makes sense. The link you provide is relevant. $\endgroup$ – dearN Aug 8 '15 at 17:28
  • 1
    $\begingroup$ Its because you have not initialized your own parameters. $\endgroup$ – chris Aug 8 '15 at 18:11
  • $\begingroup$ Yes, you are correct. Initializing my parameters allowed it to run correctly. However, I am still curious as to how I could use Regions and FEM to solve such a problem. $\endgroup$ – dearN Aug 8 '15 at 18:12
2
$\begingroup$

In the course of considering question 94639, I noticed that the steady-state problem above can be solved analytically by using DSolve instead of NDSolve.

sol = DSolve[{w''''[r] + (2/r) w'''[r] - (1/(r^2)) w''[r] + (1/(r^3)) w'[r] == -p0/De, 
    w[a] == 0, w'[a] == 0, 
    -(Derivative[1][w][b]/b^2) + Derivative[2][w][b]/b + Derivative[3][w][b] == 0, 
    (ν Derivative[1][w][b])/b + Derivative[2][w][b] == 0}, w[r], r]
%[[1, 1]] // FullSimplify
(* w[r] -> (p0 (-16 a^2 b^4 (1 + ν) Log[a]^2 - (a - r) (a + r) (-(a^2 + b^2) 
(a^2 - 6 b^2 - r^2) + (a - b) (a + b) (a^2 - 2 b^2 - r^2) ν + 8 b^4 (1 + ν) Log[b])
- 4 b^2 (-a^2 (b^2 + 2 r^2) (-1 + ν) + a^4 (1 + ν) + 2 b^2 r^2 (1 + ν) + 4 a^2 b^2 (1 + ν)
Log[b]) Log[r] + 4 a^2 b^2 Log[a] (-2 r^2 (-1 + ν) + a^2 (1 + ν) + b^2 (3 + ν) + 
4 b^2 (1 + ν) (Log[b] + Log[r]))))/(64 De (a^2 (-1 + ν) - b^2 (1 + ν))) *)

Inserting the constants provided in the question further simplifies the expression to

(* w[r] -> 2.44965*10^-6 + 0.082687 r^2 - 83.3333 r^4 + 
(4.79836*10^-7 + 0.0166667 r^2) Log[r] *)

which when plotted gives

enter image description here

as expected.

$\endgroup$
  • $\begingroup$ Thank you! I never tried DSolve for this problem. $\endgroup$ – dearN Sep 24 '15 at 13:38

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.