Dear Yves,
Thank you for your quick answer.
I tried H = [-2 1;0 0] but when the condition I found is not respected, the
boundary is completely blocked.
It only slides when this condition is respected : A^2=k*A.
For example :
A =
10.9899 -3.3151
-3.3151 1.0000
If this property is not true, the boundary is blocked. I came to this
conclusion trying a lot of combinations.
The problem is that I can't find the relation between the matrix definition and
the sliding vector.
I want to precise that I use GetFem without having installed the full version
(basic installation).
I attach my code below.
Simon
rot_mesh = gfMesh('from string', StringMesh);
%% Initialisation
gf_workspace('clear all')
% Numerical parameters
global E_rotor Nu_rotor
E = E_rotor; Nu = Nu_rotor;
lambda = E*Nu/((1+Nu)*(1-2*Nu));
mu = E/(2*(1+Nu));
%% Declaration of the mesh to GetFem
rot_mesh = gfMesh('from string', StringMesh);
%% Setting the var description
mim=gfMeshIm(rot_mesh); set(mim, 'integ',gfInteg('IM_TRIANGLE(6)'));
mfu=gfMeshFem(rot_mesh,2); set(mfu, 'fem',gfFem('FEM_PK(2,1)'));
mfd=gfMeshFem(rot_mesh); set(mfd, 'fem',gfFem('FEM_PK(2,1)'));
mf0=gfMeshFem(rot_mesh); set(mf0, 'fem',gfFem('FEM_PK(2,0)'));
mfdu=gfMeshFem(rot_mesh); set(mfdu, 'fem',gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
%% Boundary calculation
P=get(rot_mesh, 'pts');
pidlow=find(abs(P(2,:))<1e-6); %the bottom boundary can also be found using
'LowEdgeNodes' var
flow =get(rot_mesh,'faces from pid',pidlow);
ftop=get(rot_mesh,'faces from pid',(TopEdgeNodes'+1));
%% Assign boundary numbers
LOW_BOUND = 1; TOP_BOUND = 2;
rot_mesh.set_region(LOW_BOUND, flow);
rot_mesh.set_region(TOP_BOUND, ftop);
%% MODEL
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
gf_model_set(md, 'add initialized data', 'lambda', lambda);
gf_model_set(md, 'add initialized data', 'mu', mu);
gf_model_set(md, 'add isotropic linearized elasticity brick', mim, 'u',
'lambda', 'mu');
%% Centrifugal force source term
% 2D Centrifugal Force Model : F = Rho * W^2 * Distance_with_center .*[x or
y].*Projection_On_[x or y]
FXX = get(mfd, 'eval',
{[num2str(Rho_rotor),'.*',num2str(W),'.^2.*(x.^2+y.^2).^0.5
.*x./((x.^2+y.^2).^0.5)']});
FYY = get(mfd, 'eval',
{[num2str(Rho_rotor),'.*',num2str(W),'.^2.*(x.^2+y.^2).^0.5
.*y./((x.^2+y.^2).^0.5)']});
gf_model_set(md, 'add initialized fem data', 'VolumicData', mfd, [FXX;FYY]);
gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
% Referencing cource terms to the right number
gf_model_set(md, 'add initialized data', ['VolumicData'
num2str(Init_For_Mag_Term+i)], [Magn(i).PressureX,Magn(i).PressureY]);
gf_model_set(md, 'add source term brick', mim, 'u', ['VolumicData'
num2str(Init_For_Mag_Term+i)],Init_For_Mag_Term+i);
%% Sliding dirichlet conditions
% H matrix building : H_top is a projector (Hat) mathix : H^2=k*h and must be
symetrical
% The H matrix rules the sliding angle of the boundary
x4 = 1;x1 =TopAngle;k = x4+x1;x2 = sqrt(x1)*(-sqrt(k-x1));x3 = x2;
gf_model_set(md, 'add initialized data', 'H_LOW', [0 0;0 1]);
gf_model_set(md, 'add initialized data', 'VECTOR_LOW', [0;0]);
gf_model_set(md, 'add initialized data', 'H_TOP', [x1 x2;x3 x4]);
gf_model_set(md, 'add initialized data', 'VECTOR_TOP', [0;0]);
gf_model_set(md, 'add generalized Dirichlet condition with multipliers', mim,
'u', mfu, LOW_BOUND,'VECTOR_LOW', 'H_LOW');
gf_model_set(md, 'add generalized Dirichlet condition with multipliers', mim,
'u', mfu, TOP_BOUND,'VECTOR_TOP', 'H_TOP');
%% Solver
gf_model_get(md, 'solve');
U = gf_model_get(md, 'variable', 'u');
VM = gf_model_get(md, 'compute isotropic linearized Von Mises or Tresca', 'u',
'lambda', 'mu', mfdu);
Max_VM = max(VM); %Pa
Max_Displacement = max(U); %m
-----Message d'origine-----
De : Yves Renard <[email protected]>
Envoyé : mercredi 31 octobre 2018 09:43
À : SIMON AMEYE - U510180 <[email protected]>
Cc : getfem-users <[email protected]>
Objet : Re: [Getfem-users] Sliding angle for Dirichlet condition
>>> Real sender address / Reelle adresse d expedition :
>>> [email protected] <<<
**********************************************************************
Dear Simon,
A priori, H is not necessarilly a hat matrix. In that case, you can juste take
an orthognal vector, for instance [-2 1] and set H = [-2 1;0 0]], this should
work, I think.
Best regards,
Yves
----- Original Message -----
From: "SIMON AMEYE" <[email protected]>
To: "getfem-users" <[email protected]>
Sent: Tuesday, October 30, 2018 1:52:54 PM
Subject: [Getfem-users] Sliding angle for Dirichlet condition
C1-Non sensitive
________________________________
Hi all,
It has been a long time I have tried to understand the H matrix for Dirichlet
conditions.
I work with a 2D mesh, and I would like my boundary to slide freely according
to a vector, let's say [1 2].
I think I understood that H needs to be a Hat matrix.
How to do so ?
I tried to construct a hat matrix this way, but the "Angle" is not respected :
x4 = 1;
x1 = x4/(pi/2/Angle);
k = x4+x1;
x2 = sqrt(x1)*(-sqrt(k-x1));
x3 = x2;
gf_model_set(md, 'add initialized data', 'H_TOP', [x1 x2;x3 x4]);
gf_model_set(md, 'add initialized data', 'VECTOR_TOP', [0;0]); gf_model_set(md,
'add generalized Dirichlet condition with multipliers', mim, 'u', mfu,
TOP_BOUND,'VECTOR_TOP', 'H_TOP');
Simon
