Fix/Pin nodes as boundary condition

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 + nu31
nu23) / (E2E3Aa);
real C13 = (nu31 + nu21nu32) / (E2E3Aa);
real C22 = (1 - nu13
nu31) / (E1E3Aa);
real C23 = (nu32 - nu12nu31) / (E1E3Aa);
real C33 = (1- nu12
nu21) / (E1E2Aa);
real C44 = 2G23;
real C55 = 2
G13;
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

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.