|
From: | SIMON AMEYE |
Subject: | [Getfem-users] TR: RE: Sliding angle for Dirichlet condition |
Date: | Wed, 31 Oct 2018 10:32:01 +0000 |
Dear all, I finally found thanks to Yves. In my case, the Dirichlet condition matrix
matrix needs to validate the condition :
H^2=k*H If not, the sliding is blocked. H = [x1 x2; x3 x4] Moreover, [x1 x2] needs to be a vector orthogonal with the sliding direction of the boundary with the Dirichlet condition. Here is a way I set my H matrix, with ‘Angle’ the angle between the sliding direction and [1 0]. x2 = -1; x1 =tan(Angle); x4 = (x2./sqrt(x1))^2; x3 = x2; Thank you again for your help, Simon De : SIMON AMEYE - U510180
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 <address@hidden>
Envoyé : mercredi 31 octobre 2018 09:43 À : SIMON AMEYE - U510180 <address@hidden> Cc : getfem-users <address@hidden> Objet : Re: [Getfem-users] Sliding angle for Dirichlet condition >>> Real sender address / Reelle adresse d expedition : >>> address@hidden <<< ********************************************************************** 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" <address@hidden> To: "getfem-users" <address@hidden> 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 |
[Prev in Thread] | Current Thread | [Next in Thread] |