Fix/Pin nodes as boundary condition

General Discussion
1

Hello,

I would like to know if it is possible to pin nodes (instead of surfaces) as boundary conditions to solve thermoelasticity problems.
I want to get the geometric distorsion of the part from the thermal shrinkage. However, I don’t want to pin a surface but just 3 nodes to avoid rigid body translation/rotation.

I looked and try different solution proposed here but it doesn’t work properly.

There is two options :

  • Pin different nodes
  • Extract an element surface and apply the boundary condition only in this surface.

I can illustrate the exemple if necessary.
Here some part of my code : (no shrinkage, it is just to try with a pressure)

mesh3 MTPt = cube(1, 1, 1, [x0.3-0.15, y0.3-0.15, z*0.1]);
fespace TP(MTPt, [P1, P1,P1]);
TP [ut, vt, wt], [uut, vvt, wwt];

real E1 = 5.5e10;
real E2 = 30e8;
real E3 = 10e8;
real nu12 = 0.33;
real nu13 = 0.33;
real nu23 = 0.33;
real nu21 = 0.33;
real nu31 = 0.33;
real nu32 = 0.33;
real G23 = 5.2e9;
real G13 = 5.2e9;
real G12 = 5.2e9;

real Aa = 1 - nu12nu21 - nu23nu32 - nu31nu13 - 2nu21nu32nu13;
Aa = Aa / ( E1E2E3);
real C11 = (1- nu23nu32) / (E2E3Aa);
real C12 = (nu21 + nu31nu23) / (E2E3Aa);
real C13 = (nu31 + nu21nu32) / (E2E3Aa);
real C22 = (1 - nu13nu31) / (E1E3Aa);
real C23 = (nu32 - nu12nu31) / (E1E3Aa);
real C33 = (1- nu12nu21) / (E1E2Aa);
real C44 = 2G23;
real C55 = 2G13;
real C66 = 2*G12;

macro Strain(u,v,w) [dx(u), dy(v), dz(w), sqrt(2)(dy(w)+dz(v)), sqrt(2)(dx(w)+dz(u)), sqrt(2)*(dy(u)+dx(v))]//
func Cc = [ [C11, C12, C13, 0, 0, 0 ] ,
[C12, C22, C23, 0, 0, 0 ] ,
[C13, C23, C33, 0, 0, 0 ] ,
[0, 0, 0, C44, 0, 0 ] ,
[0, 0, 0, 0, C55, 0 ] ,
[0, 0, 0, 0, 0, C66] ];
real Pressure = -1;

solve mechanics([ut,vt,wt],[uut,vvt,wwt])=
int3d(MeshTP)(Strain(uut,vvt,wwt)’(CcStrain(ut,vt,wt)))

  • int2d(MeshTP, 6)(Press*wwt)
  • on(5, wt=0)
  • on(1, vt=0)
  • on(2, ut=0);
    savevtk(“presse.vtu”,dataname=“u v w tot”,MeshTP, u,t vt, wt, [ut,vt,wt], order=Order, bin=false);

I also tried with a varf approach :

varf a ( [ut, vt, wt],[uut, vvt, wwt] ) = int3d ( MTPt ) ( Strain(uut,vvt,wwt)'CcStrain(ut,vt,wt) )
+ int2d(MTPt, 6)(Pressure*wwt)
+ on(5, wt = 0)
+ on(1, vt = 0)
+ on(2, ut = 0) ;

matrix A = a ( TP , TP, solver=sparsesolver ) ; // build the matrix
real[int] b = a(0, TP);
ut[] = A^-1*b;

which leads to the same result obviously.

I then tried to reduce the matrix A by removing the row and colum associated to my pinned nodes. Same for b. But I don"t get the result I expected …

Thank you for your help.

DENIS Yvan

2

Remark, in FreeFEM all boundary condition is taken on surface in 3d not on DoF to remove lot of error
.

If you want to remove only 3 rigid body modes in plan 0xy,

in the variational formution add

solve mechanics([ut,vt,wt],[uut,vvt,wwt])=
int3d(MeshTP)(Strain(uut,vvt,wwt)’*(CcStrain(ut,vt,wt))  
  + epsilon*( [ut,vt]'*[ uut,vvt] + (dx(vt)-dy(ut))*((dx(vvt)-dy(uut)) )

where epsilon is small but not too small (1e-6) and the solution do not depend of epsilon if epsilon is small.

The term [ut,vt]’*[ uut,vvt] imposes no translation in x and y and the

(dx(vt)-dy(ut))((dx(vvt)-dy(uut)) == curl([ut,vt,wt])_3curl[uut,vvt,wwt])_3 impose not rotation around axe z.

3

Thank you very much for your answer. I will give it a try and check the result.

May I ask why BC cannot be applied in DoF in FreeFEM ? If we are using the varf approach, we are using an implicit scheme which should allow such method ?

Regards,

Y

4

It is because , lot of person make crazy things and lots of bug.

5

You can always zero out the row except for diagonal element and set the rhs.
I guess by crazy the idea is that an isolated dof depends on the mesh which should
not be critical for the end result but at least maybe for sensitivity testing
you want to play some games.

6

@marchywka Let’s play a game then !

More seriously thanks for your answers. I will make some tries on my side.

7

@marchywka FYI, I’ve spent some times on the theory and managed to apply BC on nodes instead of FE faces. After making some tries and comparison it looks like it is working well. I think the method is not numerically optimized but has its merits! Thanks again for your advices.

image
image602×585 73.3 KB

Y.

1 Like
8

To make this you can follow the script