When solving some (vector) BVP for 2d elasticity, I’d like to specify a Dirichlet boundary condition for the normal component of a vector field [u,v], as, say, …+on(1, u*N.x+v*N.y=0).
I understand that the syntax of on command does not allow this. Is there a good workaround? Thanks!
1 Like
Hi @michaellevitin , I have run up against the same problem as yours and I got around it using a penalty-based approach.
If you are integrating your PDE on a 2D domain Th, and you are interested in imposing a homogeneous Neumann boundary condition on a portion of the boundary labelled as 1, you can try
int1d(Th,1)(tgv*dn(u)*v)
where the normal derivative is defined via the following macro
macro dn(u)(dx(u)*N.x+dy(u)*N.y)//EOM.
I hope this helps.
2 Likes
Thanks, Paolo! It is not exactly what I need - I am solving a vector eigenvalue problem for elasticity, see FreeFem manual and want to impose mixed boundary conditions: tangential component of traction to vanish (which is automatic due to the weak formulation), and the normal component of the vector field [u,v], that is u*N.x+v*N.y to vanish as well. It’s the second one I’m struggling with.
M.
@michaellevitin Hi, do you work out this problem? I encountered the same problem with you. I also use the code on(1, unx+vny=0), but it does not work. I solve the Navier slip boundary condition at curved wall. the tangential velocity is given by 1/bvt, and the normal velocity vanishes, which means that unx+vny=0 at curved wall.
Could you tell me how to address this. I have found that the FreeFEM can handle this problem. A Ph.D. thesis entitled Effective Slip Lengths for Stokes Flow over Rough, Mixed-Slip Surfaces use the FreeFem++ to solve the Navier slip BC and no-penetration BC on curved wall. I contact the author but recieve no reply.
Many thanks.
@michaellevitin
Dear Michael
Did you find an answer for your question?
I have the same trouble as I have the BC on 1 I have a jet normal to the boundary with the amplitude of Aj.
I don’t know how to define it.
Best regards,
1 Like
Dear all,
I have the same problem as @michaellevitin : I have a 3D elasticity problem with symmetry conditions on a equilateral triangle. Displacement vector is [ux,uy,uz]. On the base “1” of the triangle (y=0, x=[-a/2 a/2]) I can simply add …+on(1, uy=0). However on the other two sides of the triangle, being not aligned to x or y axis, I tried …+on(2, uxN.x+uyN.y) but of course is not working. In fact I should put into relations ux with uy and therefore it is not possible to simply putting the diagonal elements to a tgv value.
Is there a way to impose symmetry conditions (e.g., Dirichlet boundary condition for the normal component of a vector field) on boundaries that are NOT aligned to x, y and z axis?
Thanks very much in advance
The problem the put fixed (Dirichlet) boundary condition on normal componant and Neuman on tangent for example.
At math level the correct (not to bad) weak formulation is take basic fonction in space Vh0
of function with a zero on u.N =0 on boundary node. If N (the normal ) is not a parallel of axis then
you have a coupling of all componante.
You can use a kind of penalization to enforce the u.N =0 by add term like in 2d
- int1d(Th)( penal* (N’[u1,u2])( (N’*[v1,v2]))
where penal est a big constant but not to big …
a small example for Stokes problem
and the associed example
Stokes-BC.edp (1.9 KB)
2 Likes
Hi Paolo, what is tgv? I’ve seen it mentioned in parts of the documentation related to Dirichlet BCs, but I’m not sure of what it is. Is it a very large value?
It is a technique to in force simply Dirichlet Boundary Condition. (call BC)
See the section 7.6 (sorry in french) of https://www.ljll.math.upmc.fr/hecht/ftp/NM406/NotesdeCoursMN406.pdf
otherwise if the make real penalization not exact (only on the diagonal term of the matrix )
you need to have enough digit in you floating number to take 2 things ons for the BC and one for the equation. (see Penalty method - Wikipedia )
Hello, Prof. Hecht,
I’m facing a similar situation, but something seems to go wrong.
Here, do you mean the following?
macro u [u1,u2] // displacement vector
macro v [v1,v2] // test function for displacement
macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))] // strain tensor
// define governing equation
problem gov(u,v)
=
-int2d(Sh)((D*e(u))'*e(v))
+int1d(Sh,2)(penaltyParameter * (u1*N.x + u2*N.y) * (v1*N.x + v2*N.y)) // u.N=0
+on(1,u1=0.1*N.x,u2=0.1*N.y) // other boundary conditions
;
I tried the above, but did not get the expected results.
The following is the result plotting [u1,u2].
In this case, the results seem reasonable, but if we only rotate the boundaries and domains, we see a different result as follows.
Below is the code I ran.
2Delastic-slider-Single-240319a.edp (2.1 KB)
Do you have any ideas how to solve this strange behavior?
Thank you for your kind cooperation.
There are two mistakes in your code:
– the coef 4*mu should be 2.*mu, this is why your code is not frame invariant
– the sign of the first term in your variational formulation should be “+” instead of “-”
Also probably you miss a constant sqrt(2.):
macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))/sqrt(2.)]
1 Like
You can use the Lagrange multiplier method if you don’t like the penalization.
For example, let \Omega be a bounded domain of \mathbb{R}^d with boundary \partial\Omega = \overline{\Gamma_N \cup \Gamma_S} and consider the following linear elasticity problem:
-\mathrm{div}(C:\varepsilon(u)) = 0 \text{ in }\Omega
\\
(C:\varepsilon(u))n = t \text{ on }\Gamma_N
\\
u\cdot n = 0 \text{ on }\Gamma_S
with the standard fourth-order elasticity tensor C, given traction t\in H^{-1/2}(\Gamma_N)^d and second-order strain tensor \varepsilon(u). This problem can be rewritten as the variational problem on H^1(\Omega)^d:
\min_{u\in H^1(\Omega)^d} \frac{1}{2}\int_\Omega \varepsilon(u):C:\varepsilon(u) dx - \int_{\Gamma_N} t\cdot u\mathrm{d}s
\\
\text{subject to }u\cdot n = 0 \text{ on }\Gamma_S.
To solve this, we define the Lagrange functional \mathcal{L}:H^1(\Omega)^d\times H^{-1/2}(\Gamma_S)\to\mathbb{R} by
\mathcal{L}(u,\varphi) = \frac{1}{2}\int_\Omega \varepsilon(u):C:\varepsilon(u) dx - \int_{\Gamma_N} t\cdot u\mathrm{d}s + \langle \varphi,u\cdot n\rangle_{H^{1/2}(\Gamma_S)},
where \langle\cdot,\cdot\rangle_{H^{1/2}(\Gamma_S)} is the usual duality pairing. A solution of the variational problem is sought as a stationary point of \mathcal{L}, i.e.,
0 = \langle d_u\mathcal{L}(u,\varphi),v\rangle_{H^1(\Omega)^d} = \int_\Omega \varepsilon(u):C:\varepsilon(v) dx - \int_{\Gamma_N} t\cdot v\mathrm{d}s + \langle \varphi,v\cdot n\rangle_{H^{1/2}(\Gamma_S)},
\\
0 = \langle d_\varphi\mathcal{L}(u,\varphi),\psi\rangle_{H^{-1/2}(\Gamma_S)} = \langle \psi,u\cdot n\rangle_{H^{1/2}(\Gamma_S)},
for all (v,\psi)\in H^1(\Omega)^d\times H^{-1/2}(\Gamma_S). In summary, our problem is to find (u,\varphi)\in H^1(\Omega)^d\times H^{-1/2}(\Gamma_S) such that
\int_\Omega \varepsilon(u):C:\varepsilon(v) dx + \langle \varphi,v\cdot n\rangle_{H^{1/2}(\Gamma_S)} + \langle \psi,u\cdot n\rangle_{H^{1/2}(\Gamma_S)} = \int_{\Gamma_N} t\cdot v\mathrm{d}s
for all (v,\psi)\in H^1(\Omega)^d\times H^{-1/2}(\Gamma_S).
To solve this using FreeFEM, we need an fespace corresponding to the product space H^1(\Omega)^d\times H^{-1/2}(\Gamma_S). Thanks to the new feature in the latest release, this can be easily implemented.
lagrange.edp (2.4 KB)
penalty.edp (2.1 KB)
1 Like