7
$\begingroup$

I have the following problem: I am currently doing calculations with some physical background in research and need a short and sweet finite element simulation of a 2D deformation problem. Usually I use Mathematica for fitting and evaluation of scattering data and the development and application of various scattering models. Hence, I am a bit lost when it comes to the use of Mathematica in structural mechanics. In essence, radial forces are applied in a circular hole in an arbitrarily shaped plate (please note, that there are actually multiple holes in the structure).

My problem is the definition of the Neumann conditions on the circular cut out, which apply these forces. By definition these forces act on the wall of the cut-out towards or from the center and should consequently extend or contract the structure. Depending if I apply contractive or expansive forces, the structure deforms into the direction of the force-value specified by the Neumann Value and not radially.

Hence my question: How would one go about programming these forces, as they apparently aren't projected correctly on the normalvector of the boundary element?

I think there is definitely an issue with my approach concerning the definition or the implementation of the boundary conditions. I am actually wondering if markers could do the trick? Still, the set up is more or less as follows:

s = RegionDifference[pol, circ[[1]]];
For[i = 2, i <= Length[all], i++,
s = RegionDifference[s, circ[[i]]];
];


Clear[mesh];
mesh = ToElementMesh[s,"MaxCellMeasure" -> 1.5, "MaxBoundaryCellMeasure" -> 
1.5,"MeshOrder" -> 2, PrecisionGoal -> 3, MeshQualityGoal -> 1, 
"MeshElementBlocks" -> 10, "NodeReordering" -> True, 
"ImproveBoundaryPosition" -> True];`

Now the definition of the input:

nr = ToNumericalRegion[s];
sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[mesh]];
vd = NDSolve`VariableData[{"DependentVariables", 
 "Space"} -> {{u, v}, {x, y}}];
methodData = InitializePDEMethodData[vd, sd]

(*Definition of the operator and elacstic constants, initialization \
of the coefficients and descritization of the PDE, ps assigns the stress 
operator*)

em = 150*10^3;

ps = {Inactive[Div][({{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 
\(1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] +   
Inactive[Div][({{-(Y/(1 - \[Nu]^2)), 0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \
[Nu]^2)))}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
Inactive[Div][({{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \\[Nu])/(1 
- \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] +Inactive[     
Div][({{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),0}, {0, -(Y/(1 - \
[Nu]^2))}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}]} /. {Y -> em, \[Nu] -> 
0.28};

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 -> em, \[Nu] 
-> 0.28};

initCoeffs = InitializePDECoefficients[vd, sd, {diffusionCoefficients}]
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];

Now the definition of the Neumann boundary conditions: (Note, that there are multiple holes in the actual structure)

Subscript[\[CapitalGamma], Nv1] = 
NeumannValue[ns[[1]], (y - all[[1, 2]])^2 + (x - all[[1, 1]])^2 <= (rs[[1]] 
+ .7)^2];
Subscript[\[CapitalGamma], Nu1] = 
NeumannValue[ns[[1]], (y - all[[1, 2]])^2 + (x - all[[1, 1]])^2 <= (rs[[1]]   
+ .7)^2];

To set the structure into place, several Dirichlet conditions are specified at the corners of the polygon:

Subscript[\[CapitalGamma], Dv] = DirichletCondition[v[x, y] == 0., y <= -8.2 
&& x <= -4.25];
Subscript[\[CapitalGamma], Dv1] =DirichletCondition[v[x, y] == 0., y <= -8.2 
&& x >= 4.25];
Subscript[\[CapitalGamma], Dv2] =DirichletCondition[v[x, y] == 0., y >=  8.2 
&& x >= 4.25];
Subscript[\[CapitalGamma], Dv3] =DirichletCondition[v[x, y] == 0., y >= 8.2 
&& x <= -4.25];
Subscript[\[CapitalGamma], Dv4] =DirichletCondition[v[x, y] == 0., x <= 
-9.2];
Subscript[\[CapitalGamma], Dv5] =DirichletCondition[v[x, y] == 0., x >=  
9.2];
Subscript[\[CapitalGamma], Du] =DirichletCondition[u[x, y] == 0., y <= -8.2 
&& x <= -4.25];
Subscript[\[CapitalGamma], Du1] =DirichletCondition[u[x, y] == 0., y <= -8.2 
&& x >= 4.25];
Subscript[\[CapitalGamma], Du2] =DirichletCondition[u[x, y] == 0., y >=  8.2 
&& x >= 4.25];
Subscript[\[CapitalGamma], Du3] =DirichletCondition[u[x, y] == 0., y >= 8.2 
&& x <= -4.25];
Subscript[\[CapitalGamma], Du4] =DirichletCondition[u[x, y] == 0., x <= 
-9.2];
Subscript[\[CapitalGamma], Du5] =DirichletCondition[u[x, y] == 0., x >=  
9.2];

initBCs = InitializeBoundaryConditions[vd, sd, {{Subscript[\[CapitalGamma], 
Du], Subscript[\[CapitalGamma],Du1], Subscript[\[CapitalGamma], Du2], 
Subscript[\[CapitalGamma],Du3], Subscript[\[CapitalGamma], Du4], Subscript[\
[CapitalGamma],Du5], Subscript[\[CapitalGamma], Nu1]}, {Subscript[\
[CapitalGamma], Dv],Subscript[\[CapitalGamma], Dv1], Subscript[\
[CapitalGamma], Dv2],Subscript[\[CapitalGamma], Dv3], Subscript[\
[CapitalGamma], Dv4],Subscript[\[CapitalGamma], Dv5], Subscript[\
[CapitalGamma], Nv1]}}]

The solution is calculated via LinearSolve, hence the definition of the system matrices:

{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]

Application of the boundary conditions to the mesh via discretization:

discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

DeployBoundaryConditions[{load, stiffness}, discreteBCs];

Now the solution stage of the problem:

 Short[MemoryConstrained[solution = LinearSolve[stiffness, load],8*10^9]]

Depending on the sign of my force ns, the structure deforms accordingly:

"Negative Load"

And finally also:

"Positive Load"

Actually, a contraction/expansion is expected.

Thank you very, very much in advance for any kind of help!

Best regards,

u53r

Edit

I found a solution to the problem, which might be quite trivial for users on this site. For a circle, one simple defines the projection of the force on the surface. The direction of the force vector needs to be known, which in the case of a radial, uniform distribution on a circular cut-out reduces to:

Subscript[\[CapitalGamma], Nu] = 
NeumannValue[ns[[1]]*Dot[{1,0}, {Cos[ArcTan[(center[[1, 1]] - x), 
(center[[1, 2]] - y)]],Sin[ArcTan[(center[[1, 1]] - x), (center[[1, 2]] - 
y)]]}], (y - center[[1,2]])^2 + (x - center[[1, 1]])^2 <=  (radius[[1]] + 
.7)^2];

For the v direction the same follows.

For arbitrary curves on would need to describe it implicitly and calculate the normal vector via:

Then this vector will be projected into the u and v direction with the {1,0} or {0,1} vector. I didn't program this out, but it should work this way.

Best regards,

u53r

$\endgroup$
5
$\begingroup$

Here is a way to do it directly with NDSolve. No need to get your hands dirty unless you need to.

I use a slightly different region; it's no restriction.

r = RegionDifference[Rectangle[{-2, -2}, {2, 2}], Disk[]];

Define the plane stress operator operator:

op = {Inactive[Div][{{0, -((nu*Y)/(1 - nu^2))}, 
             {-((1 - nu)*Y)/(2*(1 - nu^2)), 0}} . 
           Inactive[Grad][v[x, y], {x, y}], {x, y}] + 
       Inactive[Div][{{-(Y/(1 - nu^2)), 0}, 
             {0, -((1 - nu)*Y)/(2*(1 - nu^2))}} . 
           Inactive[Grad][u[x, y], {x, y}], {x, y}], 
     Inactive[Div][{{0, -((1 - nu)*Y)/(2*(1 - nu^2))}, 
             {-((nu*Y)/(1 - nu^2)), 0}} . Inactive[Grad][
             u[x, y], {x, y}], {x, y}] + 
       Inactive[Div][{{-((1 - nu)*Y)/(2*(1 - nu^2)), 0}, 
             {0, -(Y/(1 - nu^2))}} . Inactive[Grad][v[x, y], 
             {x, y}], {x, y}]} /. {Y -> 10^3, nu -> 33/100};

Boundary conditions: Fix the structure that the corner points of the rectangle region.

Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0., 
    v[x, y] == 
     0.}, (x == -2 && y == -2) || (x == -2 && y == 2) || (x == 2 && 
      y == -2) || (x == 2 && y == 2)];

Specify a force that acts in the inner circle:

force = -100;

Add a NeumannValue on the boundary for the x- and the y-direction based on your idea and solve:

{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[force*Cos[ArcTan[x, y]], 
       x^2 + y^2 == 1], 
      NeumannValue[force*Sin[ArcTan[x, y]], x^2 + y^2 == 1]}, 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] r];

Show the deformed structure:

mesh = ufun["ElementMesh"];
Show[{
  mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
  NDSolve`FEM`ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

enter image description here

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ can you let me know how to calculate the elastic strain energy in this cad part using intergral? $\endgroup$ – ABCDEMMM Jul 11 '18 at 7:41
  • 1
    $\begingroup$ @ABCDEMMM, I think you should ask a question for that and make sure you include the code that you have so far. $\endgroup$ – user21 Jul 11 '18 at 7:48
  • 1
    $\begingroup$ @ABCDEMMM, no one can answer that question to your satisfaction because it is missing information, the code you have so far etc. For example add a mathematical and mathematica definition of the energy integral you want to solve. Where do you want to solve this. Unless you write a better question you will not get a good answer. Make it easy for people to help you by not having them look up thinks like the mathematical definition of the energy integral you want to solve. $\endgroup$ – user21 Jul 11 '18 at 8:13

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.