8

I work on geological problems and I use the Finite Element Method. But this question can be applied on other classical mechanical problems.

I work on implicit 3D surfaces (which represent the limits between two geological layers aka two media).

I have to impose displacement on these surfaces (Dirichlet Condition). The condition that I impose on surfaces are transfered on the surrounding nodes. I make the assumption that I can use the " classical" way to apply Dirichlet Condition (penalty method in my case) on these nodes, that are NOT on the the boundaries.

My question is: do you think this assumption is valid? Have you got some references about the subject? As I understand the problem I try to apply Boundary Conditions not on the boundaries...

Thank you in advance.

Best.

EDIT

We try to solve this problem: $$ \left\{ \begin{aligned} \sigma_{ij,j} + F_{i} &= 0 & &\text{in the domain $\Omega$}&\\ u_i&= q_{i}& &\text{on the boundary $\Gamma_{q}$}&\\ \sigma_{ij}n_{j} &= h_{i}& &\text{on the boundary $\Gamma_{h}$}&\\ u_i &= b_{i}& &\text{punctually in the domain $\Omega$}&\\ \end{aligned} \right. $$

We run the FEA on a mesh. We basically try to impose Dirichlet Condition on the nodes which are inside the domain $\Omega$ (fourth line). We constrain the displacement value of $u_i$ to $b_i$ for these nodes.

AntoineMazuyer
  • 95
  • 5

2 Answers2

9

The solution of the equation you are looking for is in the space $H^1$ of functions that have one weak derivative, but in 2d and 3d, this does not imply that the solution is in fact continuous. As a consequence, it is not possible to define the value of a solution at individual points, and equally consequentially, it is mathematically not possible to prescribe (boundary) values at individual points.

But I suspect that you're not actually interested in prescribing values for the displacement at individual points, but for all points along a line (in 2d) or surface (in 3d), and this is perfectly valid. (Mathematically speaking, this is so because the trace operator is well defined on $H^1$). At the discrete level, this is then equivalent to prescribing the values on all nodes that lie on this line/surface.

Conceptually, you can think of this as a crack in your domain that is occupied by an infinitely thin device whose displacement you can manipulate and that attached to the two sides of the crack.

Wolfgang Bangerth
  • 55,373
  • 59
  • 119
  • OK thank you very much for these informations ! – AntoineMazuyer Oct 05 '15 at 12:15
  • Revisiting this post, I have seen many papers (e.g. https://www.mdpi.com/2073-4360/14/8/1559/htm) applying only one Dirichlet boundary condition in 2D elasticity problem (plates). I am trying to do the same but I ended up with a stiffness matrix that is not a positive definite matrix. Does anyone know how they (also commercial packages) were able to solve the problem this way? – Riobaldo Tatarana Oct 10 '22 at 20:11
  • @ProfessorP.CosmoKlunk Does the matrix have zero eigenvalues, or in fact eigenvalues of both signs? – Wolfgang Bangerth Oct 10 '22 at 21:16
  • @WolfgangBangerth after applying the Dirichlet boundary conditions on the matrix (as you show on the deal.II's videos) the negative one is -4.12e-17 (can I consider it as zero?). The other ones are close to zero (8.08e-17) and positive. – Riobaldo Tatarana Oct 10 '22 at 23:51
  • -4e-17 is small only if it is small compared to the other eigenvalues. If the other eigenvalues are of the order +4e-15, then it is not small, for example. If the other eigenvalues are of the order +4e3, then it is small. – Wolfgang Bangerth Oct 11 '22 at 03:31
  • @WolfgangBangerth, I increased the number of nodes and now I have the lowest eigenvalue equal to 8.9e-17 and the biggest one is 1. So the matrix is semi-positive definite. However, my results are not correct. The solution looks like I have three surfaces shifted. – Riobaldo Tatarana Oct 11 '22 at 09:02
  • What does the visual inspection tell you about what the nullspace of the solution might be? – Wolfgang Bangerth Oct 11 '22 at 18:25
  • I cannot find any references expanding on your last statement. E.g. applying $u_i=b_i$ on a boundary $\Gamma_b$, that exists within the domain $\Omega$, such as a crack, how to you formulate this? From the weak form of a boundary value problem we can apply integration by parts that yields a term that integrates on the domain $\Omega$ and a term that integrates on its boundary $\partial\Omega$. The weighing functions are chosen to satisfy essential boundary conditions. But if $\Gamma_b$ is not within $\partial\Omega$, what do you do? – Breno Dec 02 '22 at 12:35
  • I guess you have to explicitly define two separate regions $\Omega_1$ and $\Omega_2$ around the crack, such that $\Omega=\Omega_1\cup \Omega_2$ and then handle the kinematics of the surface between them. Is there no way around this? – Breno Dec 02 '22 at 12:37
  • Correction: the weighing functions are chosen to be homogeneous at the application of the boundary conditions $\Gamma_q$ – Breno Dec 02 '22 at 12:42
0

Well you can model slip between geologic faults (a discontinuity) using constraint equations which can be thought of as loading at the interface instead of the boundary. You simply specify the relative slip and solve the problem. E.g., see the top-left figure in https://bitbucket.org/stali/defmod/wiki/Gallery

If that is indeed what you want then I would suggest that you use Lagrange multipliers instead of the penalty method. Off course you can use any constraint equation you want but it should make physical sense.

stali
  • 1,759
  • 1
  • 10
  • 13
  • Thank your for suggestions but my initial question is more like "Is FEA, which is based on variational formulation of boundary values problem, still a valid tool to solve a constrained equation with constraints on values inside the domain (which is no longer a boundary value problem)?" EDIT: I take a look to tour code it's pretty interesting – AntoineMazuyer Oct 01 '15 at 14:37
  • Yes it is. In theory you can think of the fault interface as a boundary (there are coincident nodes). – stali Oct 01 '15 at 14:54
  • Ok but the surfaces on which I impose the boundary conditions are horizons. (interface between two layers) – AntoineMazuyer Oct 01 '15 at 15:10