Point contact - Revision history

E6WikiAdmin: /* Numerical solution with FreeFem++ */

2022-11-23T17:39:10Z

‎Numerical solution with FreeFem++

2022-11-23T17:38:01Z

‎Explicit solution with Meshless Local Strong Form Method (MLSM)

2022-10-22T19:35:16Z

‎Implicit solution]

2022-10-22T19:35:00Z

‎Implicit solution with Meshless Local Strong Form Method (MLSM)

2018-02-07T04:50:20Z

‎Convergence

2017-09-14T16:40:13Z

‎Closed form solution

2017-09-14T16:35:14Z

‎Closed form solution

2017-09-14T01:14:06Z

‎Purely Dirichlet case

2017-09-14T01:13:18Z

‎Closed form solution

2017-09-14T01:10:43Z

‎Closed form solution

← Older revision		Revision as of 17:38, 23 November 2022
Line 89:		Line 89:


−	==Explicit solution ~~with [[Meshless Local Strong Form Method (MLSM)]]==~~	+	==Explicit solution

	<syntaxhighlight lang="c++" line>		<syntaxhighlight lang="c++" line>

@@ Line 126: / Line 126: @@
 Note that an operator grad(div(u)) that requires also mixed derivatives is used. To obtain stable solution a minimal second order monomial basis is required.
-==Implicit solution]==
+==Implicit solution==
 <syntaxhighlight lang="c++" line>
      /// [Field declarations]

@@ Line 126: / Line 126: @@
 Note that an operator grad(div(u)) that requires also mixed derivatives is used. To obtain stable solution a minimal second order monomial basis is required.
-==Implicit solution with [[Meshless Local Strong Form Method (MLSM)]]==
+==Implicit solution]==
 <syntaxhighlight lang="c++" line>
      /// [Field declarations]

← Older revision		Revision as of 04:50, 7 February 2018
Line 182:		Line 182:

	[[File:crossY.png\|400px]][[File:crossX0.png\|400px]][[File:crossX.png\|400px]]		[[File:crossY.png\|400px]][[File:crossX0.png\|400px]][[File:crossX.png\|400px]]
−
−	~~==== Convergence ====~~
−
−
−	~~[[File:Point cont conv.png\|800px]]~~
−
−	~~Note, that this are preliminary results and there might be some minor mistakes. Updates are coming soon :)~~

@@ Line 36: / Line 36: @@
      double x = p[0], y = p[1], r2 = x*x + y*y, factor = -force / (4.0*M_PI*mu);
      double ux = 2*x*y/r2 + 2*mu/(lam + mu)*std::atan2(y, x);
-     double uy = (y*y-x*x)/r2 - (lam + 2*mu)/(lam + mu)*std::log(x*x + y*y);
+     double uy = (y*y-x*x)/r2 - (lam + 2*mu)/(lam + mu)*std::log(r2);
      return Vec2d(factor*ux, factor*uy);
 };
@@ Line 43: / Line 43: @@
 Copyable Matlab code for closed form solution:
 <syntaxhighlight lang="matlab" line>
-function [sxx, syy, sxy, u, v] = cantilever_beam_analytical(x, y, P, L, D, E, nu)
+function [ux, uy] = point_contact_analytical(x, y, lam, mu, P)
-% Evaluates closed form solution to timoshenko cantilever beam.
+% Evaluates closed form solution for point contact.
-% P = loading force
+% lam and mu are Lame parameters, P is the point force magnitude
-% L = length
+r2 = x.^2 + y.^2; factor = -P / (4.0*pi*mu);
-% D = height
+ux = factor*(2*x.*y./r2 + 2*mu/(lam + mu)*atan2(y, x));
-% E = Young's modulus
+uy = factor*((y.^2-x.^2)./r2 - (lam + 2*mu)/(lam + mu)*log(r2));
 end
 </syntaxhighlight>

← Older revision		Revision as of 16:35, 14 September 2017
Line 30:		Line 30:

	The proof of correctness and plots of the solution are available in this notebook: [[:File:point_contact.nb]]		The proof of correctness and plots of the solution are available in this notebook: [[:File:point_contact.nb]]
		+
		+	Copyable C++ code for closed form solution:
		+	<syntaxhighlight lang="c++" line>
		+	std::function<Vec2d(Vec2d)> analytical = [=](const Vec2d& p) {
		+	double x = p[0], y = p[1], r2 = xx + yy, factor = -force / (4.0M_PImu);
		+	double ux = 2xy/r2 + 2mu/(lam + mu)std::atan2(y, x);
		+	double uy = (yy-xx)/r2 - (lam + 2mu)/(lam + mu)std::log(xx + yy);
		+	return Vec2d(factorux, factoruy);
		+	};
		+	</syntaxhighlight>
		+
		+	Copyable Matlab code for closed form solution:
		+	<syntaxhighlight lang="matlab" line>
		+	function [sxx, syy, sxy, u, v] = cantilever_beam_analytical(x, y, P, L, D, E, nu)
		+	% Evaluates closed form solution to timoshenko cantilever beam.
		+	% P = loading force
		+	% L = length
		+	% D = height
		+	% E = Young's modulus
		+	% nu = Poisson ratio
		+	I = D^3 / 12;
		+	sxx = I.^(-1).P.x.*y;
		+	syy = zeros(size(x));
		+	sxy = (1/2).I.^(-1).P.((1/4).D.^2+(-1).*y.^2);
		+
		+	u = (1/24).E.^(-1).I.^(-1).P.y.(3.D.^2.*(1+nu)-...
		+	4.(3.L.^2+(-3).x.^2+(2+nu).y.^2));
		+	v = (-1/24).E.^(-1).I.^(-1).P.(3.D.^2.(1+nu).(L+(-1).x)+...
		+	4.(2.L.^3+(-3).L.^2.x+x.^3+3.nu.x.*y.^2));
		+	end
		+	</syntaxhighlight>

	==Numerical solution with [http://www.freefem.org/ FreeFem++]==		==Numerical solution with [http://www.freefem.org/ FreeFem++]==

@@ Line 138: / Line 138: @@
 == Results ==
 === Purely Dirichlet case ===
-For starters we check the solution of the problem with Dirichlet BC. The conditions are obtained from closed form solution [[#Point contact on a 2D half-plane]].
+For starters we check the solution of the problem with Dirichlet BC. The conditions are obtained from closed form solution [[#Closed form solution]].
 [[File:image_1b19ssjsn1e4t7qvq0ffnb1pft9.png|600px]]

@@ Line 26: / Line 26: @@
 u_y &= -\frac{P}{4\pi\mu}\left(\frac{y^2-x^2}{r^2}-\frac{\lambda+2\mu}{\lambda + \mu} \log(r^2)\right), \label{eq:dispy}
 \end{align}
-where again $r^2 = x^2+y^2$ and $\lambda$ and $\mu$ are Lame parameters. For plane stress, do the substitution $\lambda \gets \frac{2\mu\lambda}{\lambda+2\mu}$ in the solution and in the Navier equation.
+where again $r^2 = x^2+y^2$ and $\lambda$ and $\mu$ are Lame parameters. For plane stress, do the substitution
 ==Numerical solution with [http://www.freefem.org/ FreeFem++]==

@@ Line 11: / Line 11: @@
 == Closed form solution ==
-The analytical relations for the stresses can be found from the [https://en.wikipedia.org/wiki/Flamant_solution Flamant solution] (stress distributions in a linear elastic wedge loaded by point forces a the tip. When the "wedge" is flat we get a half-plane) and are given as:
+The analytical relations for the stresses can be found from the [https://en.wikipedia.org/wiki/Flamant_solution Flamant solution] (stress distributions in a linear elastic wedge loaded by point forces a the tip. When the "wedge" is flat we get a half-plane). For a point $(x, y)$ denote its square distance to the origin by $r^2 = x^2 + y^2$, and the stresses in point $(x, y)$ are given by:
 \begin{equation}
 \sigma_{xx} = -\frac{2P}{\pi} \frac{x^2y}{r^4},
@@ Line 26: / Line 26: @@
 u_y &= -\frac{P}{4\pi\mu}\left(\frac{y^2-x^2}{r^2}-\frac{\lambda+2\mu}{\lambda + \mu} \log(r^2)\right), \label{eq:dispy}
 \end{align}
-where $r^2 = x^2+y^2$ and $\lambda$ and $\mu$ are Lame parameters. For plane stress, do the substitution $\lambda \gets \frac{2\mu\lambda}{\lambda+2\mu}$ in the solution and in the Navier equation.
+where again $r^2 = x^2+y^2$ and $\lambda$ and $\mu$ are Lame parameters. For plane stress, do the substitution $\lambda \gets \frac{2\mu\lambda}{\lambda+2\mu}$ in the solution and in the Navier equation.
 ==Numerical solution with [http://www.freefem.org/ FreeFem++]==