# Tag Info

1

It has been a long time since I've thought about this problem. It is probably hard to solve with NDSolve, but it works with finite differences. Start with conditions on the surface and work your way down the rod string to the bottom. General wave equation with viscous damping and gravity. pde = D[u[x, t], t, t] == a^2*D[u[x, t], x, x] - c*D[u[x, t], t] - g ...

6

This is not an answer but a comment on Tim's answer. Tim's answer is just fine as it is. However, I'd like to take the opportunity to show how to create a mesh that is an even more accurate representation of the geometry; the additional accuracy is most likely not needed in this case but it makes a nice example to show the functionality. Create a boundary ...

7

Your model appears to have quarter symmetry. If one can take advantage of symmetry, it will be a smaller model and may even be easier to setup. A good place to start to find a good setup is PDEModels Overview. Clicking on the Plane Stress will bring you to a verified operator. It could be useful to use FEMAddOns to difference two boundary meshes so that ...

10

A more suitable solution can be found in the Mathematica documentation for solving plane stress in the section structural mechanics for NDEigensystem. Specify a plane stress PDE: {vals, funs} = NDEigensystem[{ps, DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]}, {u[x, y], v[x, y]}, {x, y} ∈ Ω, 9]; vals {351.293, 369.64, 495.516, 1479.33, 2021....

1

You could precompile the mesh refine function with additional information. For example, cmrf = Compile[{{vertices, _Real, 2}, {area, _Real}}, Module[{fest, ftrue, error, tolerance, test}, fest = Mean[fDirect /@ vertices]; ftrue = fDirect[Mean[vertices]]; error = fest - ftrue; tolerance = 0.43; test = Abs[error] > tolerance; test], {{...

4

As recommended by PaulCommentary and xzczd we put bc = DirichletCondition[z[x, t] == 0, x == 1]; to the end x=1 and apply force to the end x=0. Also we use physical normalization and change phase of force in accordance with zero initial condition, then we have Y = 199*^9;(*Pa*)ρ = 7860;(*kg/m^3*)dia = 1/39.37;(*1" dia converted to meters*)c = Sqrt[Y/ρ]...

2

You assume that NeumannValue and Derivative are interchangeable; they are generally not. This is explained in more detail in the section The Relation between NeumannValue and Boundary Derivatives in the Finite Element Method User Guide. Also the NeumannValue ref page has detailed information what it models. Let's look at the message: The message states that ...

2

OK, let me extend my comment to an answer. You're solving Laplace equation, for which the new-in-v10 FiniteElement is the only available method in NDSolve at the moment. For more info about when FiniteElement will be automatically chosen, check the following post: PDEs : automatic method choice : TensorProductGrid or FiniteElement? When FiniteElement is used ...

4

To use the FEM method it is good to cast your equations into coefficient form as shown FEM Tutorial. $$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$ When you use the coefficient form, there generally is a nice one-to-one ...

3

If you want Function[{x1, x2}, x1 + x2] to be the solution of a PDE, you probably want the following PDE: ClearAll[y, x1, x2]; a = x1 + x2; pde = -D[a D[y[x1, x2], x1], x1] - D[a D[y[x1, x2], x2], x2] + 2; bc = {y[x1, 2] == 2 + x1, y[x1, 3] == 3 + x1}; Y = NDSolveValue[{pde == NeumannValue[-1*a, x1 == 2] + NeumannValue[1*a, x1 == 3], bc}, y, {x1, 2, 3}, {x2, ...

2

Here I am sketching out a solution path. My feeling is that the equations are not quite right, or something with the setup. More on this later. This setup solves the equation given. Make sure it is correct as the examples shown in the question differ slightly. g = Rationalize[0.15]; a = Rationalize[0.2]; b = Rationalize[0.2]; NN = 1500; region = Triangle[]; ...

3

By reversing the sign of the derivative on the left side from that given in NeumannValue, this can be solved by Mathematica analytically as well. ClearAll[y, x1, x2]; pde = Laplacian[y[x1, x2], {x1, x2}] == 0; bc = {y[x1, 2] == 2 + x1, y[x1, 3] == 3 + x1, Derivative[1, 0][y][2, x2] == 1, Derivative[1, 0][y][3, x2] == 1}; solA = DSolve[{...

7

Relatively recently, Wolfram has created a nice Heat Transfer Tutorial and a Heat Transfer Verification Manual. I model with many codes and I usually start the Verification and Validation manual and build complexity from there. It is always embarrassing to build a complex model and find that your setup does not pass verification. The Laplace equation is ...

13

Update: New version of FEMAddOns v1.4.2 includes ExtrudeMesh @user21 updated FEMAddOns and ExtrudeMesh is included. Now, you only need to call FEMAddOns and don't need to worry about namespace collisions with MeshTools. For completeness, here is the updated workflow to get the latest version of FEMAddOns and extend Alex's answer by two lines: ...

11

We can play with numbers to get an ideal mesh. As a stating point you could use this one: Needs["MeshTools"](*Needs["FEMAddOns"]*) Lx = 100; Ly = 50; r[t_] := Piecewise[{{Lx/Cos[t], 0 <= t <= ArcTan[Ly/Lx]}, {Ly/Sin[t], ArcTan[Ly/Lx] < t <= Pi - ArcTan[Ly/Lx]}, {-Lx/Cos[t], Pi - ArcTan[Ly/Lx] < t <= Pi + ...

2

No, unfortunately there is not. The domain decomposition solver DecompositionNDSolveValue you mention from the FEMAddOns works by decomposing the domain, solving the PDE on the subdomain and make use of this subdomain solution as a boundary condition for the solution on another subdomain in clever way. This approach does not work for NDEigensystem or ...

9

Update The real observation is that MeshRefinementFunction does not work for 1D with ToElementMesh. Yes, that's unfortunately the case but you can easily use Needs["NDSolveFEM"]; f = Function[{vertices, area}, If[Mean[vertices] > 1, area > 0.1, area > 0.01]]; mr = DiscretizeRegion[Interval[{0, 2}], MeshRefinementFunction -> f]; ...

6

There is a tutorial called NDSolve Options for Finite Elements that is the main tutorial for Finite Element related options. It names all FEM options (but one) and either gives examples and explanations for them or links to the respective location in the documentation where you can find more information. Concerning NDEigensystem specifically all options are ...

1

It would be nice to update the FEAST solver to the latest (4.0 as of 2020) version to allow for non-Hermitian problems and to benefit from the performance enhancements.

8

If you have some innovative approach in mind… Hmm… I have a solution based on finite difference method (FDM), which is said to be the earliest numeric method for PDE solving. Anyway, it might be viewed as innovative in the sence that it shows FEM isn't the only possible numeric approach, so let me post it here. FDM Based Solution I've used pdetoae for the ...

3

You could start with something like this: eqn = {Laplacian[u[x, y], {x, y}] - v[x, y], Laplacian[v[x, y], {x, y}]}; bcs = DirichletCondition[u[x, y] == 0, True]; ufun = NDEigenvalues[{eqn, bcs}, {u, v}, {x, -Pi/2, Pi/2}, {y, -Pi/2, Pi/2}, 5] This will give a warning message about a possible non positive definite matrix. Which I think is OK, but I ...

4

When you use ToElementMesh the result will always be a mesh of the full region. What you want to do in this case is to generate a boundary mesh were the embedding dimension is not the same as the region dimension: Needs["NDSolveFEM"] ToBoundaryMesh[ "Coordinates" -> Table[{100 Cos[phi], 100 Sin[phi]}, {phi, 0, Pi/4, 1/8 Pi/4}], ...

3

When a symbolic region is meshed there can be discrepancies between what the symbolic region represents and what the numerical approximation represent; that is especially the case for curved regions. "OpenCascade" typically does a good job at preserving the boundary. So if you can, it's advisable to represent the region with graphics primitives ...

1

In solid mechanics, body forces like gravity or inertial forces (like those due to acceleration) are not boundary conditions. Instead, they appear in the actual PDE itself. Surface forces (also called traction forces) are boundary conditions, as are physical constraints such as zero displacement or slope.

3

The units are defined in $N/m³$. 1.Verification Let's look at the left side of the equation. 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[...

10

I have performed the simulations with an AMD Ryzen 3900X CPU, Windows 64 bit, AceFEM version 7.103. Maybe someone can provide some comment to the obtained results and how they compare the the results on Intel i7. Direct solver simulation report: Iterative solver simulation report: EDIT: As per KratosMath request, this is a snapshot of some of my ...

0

This will give a solution but the problem is convection dominated: Needs["NDSolveFEM"]; s[x_] := UnitBox[0.3 x - 1.5] (*initial dispalacement*) ics = {u[x, 0] == s[x]}; nr = ToNumericalRegion[FullRegion, {{0, 10}}]; mesh = ToElementMesh[nr, MaxCellMeasure -> 0.25]; pde = D[u[x, t], t] + D[u[x, t], x] == 0.01*u[x, t]*D[u[x, t], x, x]; bcs = ...

5

Here is an approach based on creating exact regions: a = Rationalize[0.857597, 10^-16]; b = Rationalize[1.653926, 10^-16]; hexagon = Polygon[{{0, (b - a)/2, 1/2}, {(b - a)/2, 0, 1/2}, {1/2, 0, (b - 1)/(2 a)}, {1/2, (b - 1)/2, 0}, {(b - 1)/2, 1/2, 0}, {0, 1/2, (b - 1)/(2 a)}}] // Simplify; octahedron = ImplicitRegion[Abs[x] + Abs[y] + a ...

2

We can use mesh of first order for 3D visualization and short time for visibility. We also change boundary conditions: Needs["NDSolveFEM"]; a = ImplicitRegion[True, {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}]; b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2]; c = Cylinder[{{1, 1, -1/5}, {1, 1, 0}}, 650/1000]; d = Cylinder[{{-1, 1, -1/5}, {-1, 1, 0}}, 650/...

6

This is not ideal, but it gives an approximate resulting region. I first generate random points on the hexagon and add a random vector on the unit sphere. I take the convex hull of the points which is acceptable because the blob must be convex. Finally I discretize the octahedron and intersect with crudehexagonblob: crudehexagonblob = ConvexHullMesh[# + ...

7

Too long for a comment. An easy way to generate a high quality mesh is to replace the Implicitegion with Cubuid and make use of the OpenCascade boundary mesh generator: Needs["NDSolveFEM"] (*a=ImplicitRegion[True,{{x,-1,1},{y,-1,1},{z,0,1}}];*) a = Cuboid[{-1, -1, 0}, {1, 1, 1}]; b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2]; c = Cylinder[{{1, 1, ...

14

Update (Steady-State Solution) I think the fundamental issue is that you are over constraining your system. Whether you are solving the "heat equation" or not, your operator has the same form of the heat equation as shown below: $$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} + \nabla \cdot {\mathbf{q}} = 0$$ If the flux, $\mathbf{q}$, needs ...

5

Thinking for a moment suggests χ == 0 is a numerically problematic boundary condition for two reasons: First, the χ[x]^(3/2) term, which means that χ[x] should not go negative or the integration will run into branch cut difficulties. Second, the solution is concave up so that if it has a turning point (a positive minimum), it will increase; in fact, ...

6

There is also wavelet method for BVP. It is an example with Haar wavelets. It takes 0.36 s to solve this problem with 64 colocation points: ClearAll["Global*"] L = 10; A = 0; B = 1; J = 5; M = 2^J; dx = (B - A)/(2 M); h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}]; p1[x_, n_] := (1/n!)*(x - A)^n; h[x_, k_, m_] := Piecewise[{{1, ...

5

Although I anxiously await Andres' complete write up, I thought that I would post some observations that may help in the investigation of the PeriodicBoundaryCondition. In this case, my initial findings are that a combination @Rodion Stepanov's symmetrized PBC and a triangle mesh lead to more robust results without needing a "Ghost Vicinity". Default ...

6

Answer under construction. Beginning of explanations are coming later (2 days ?). The code below is complete, so one can already evaluate it and enjoy. Short and quick explanations are already possible in in this chatroom, but the subject is really hudge. If you see a problem or some possible simplification anywhere, don't hesitate to comment. ...

Top 50 recent answers are included