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$
6
  • 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, 2015 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, 2015 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, 2015 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, 2015 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, 2015 at 10:42

3 Answers 3

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$
12
  • $\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, 2015 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, 2015 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, 2015 at 13:52
  • $\begingroup$ Thank you for confirming my deduction of XY plane. $\endgroup$
    – dearN
    Aug 12, 2015 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, 2015 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$
3
  • $\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, 2015 at 17:28
  • 1
    $\begingroup$ Its because you have not initialized your own parameters. $\endgroup$
    – chris
    Aug 8, 2015 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, 2015 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$
1
  • $\begingroup$ Thank you! I never tried DSolve for this problem. $\endgroup$
    – dearN
    Sep 24, 2015 at 13:38

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.