Hello everyone,
I have tried to compute the reaction force by using alternative 1, i.e.
bending the beam with a force on the right side. With alternative 1 and 2
I refer to the Matlab program of my previous post (compute a reaction
force). My idea was to check the correctness of my program, because in
this case the reaction force has the same value as the applied force
(actio = reactio). I used the multiplier version for the homogenous
Dirichlet condition on region 1 (fixation on the left side):
gf_model_set(md, 'add Dirichlet condition with multipliers', ...
mim, 'u', mfu, 1);
I did not change the remaining Matlab code of my previous post. I only
added my routine to compute the reaction force at the end of the code:
% get the tangent matrix
tangent_matrix = gf_model_get(md, 'tangent_matrix');
% get the multipliers
mult = gf_model_get(md, 'variable', 'mult_on_u');
% get the number of multipliers and DOFs
nb_mult = size(mult,2);
nb_dof = gf_mesh_fem_get(mfu,'nbdof');
% part of the tangent matrix concerning the multipliers
mult_matrix = zeros(nb_dof,nb_mult);
for i = 1:nb_dof
for j = 1:nb_mult
mult_matrix(i,j) = tangent_matrix(nb_mult+i,j);
end
end
% computing the nodal forces by multiplying the multipliers
% with the right part of the tangent matrix
nodalforce = -(mult_matrix*transpose(mult));
% extract the x-, y- and z-components
for i = 0:(size(nodalforce,1)/3)-1
nodalforce_x(i+1) = nodalforce(3*i+1);
nodalforce_y(i+1) = nodalforce(3*i+2);
nodalforce_z(i+1) = nodalforce(3*i+3);
end
% sum the components for the total reaction force
reactionforce_x = sum(nodalforce_x);
reactionforce_y = sum(nodalforce_y);
reactionforce_z = sum(nodalforce_z);
Is this correct?
I get the right reaction force, when I have a beam e.g. with length L =
10, width 1 and height 1 and a force F = [0 100 0]. When testing other
versions of the beam, e.g. changing the width from 1 to 2, I get a
reaction force [0 -200 0]. The same occurred when changing the height, but
not when changing the lenght. Is my assumption right that the Source Term
brick adds some kind of normalized forces, i. e. a force applied to a unit
square? When summing up the values of the right hand side, I also get a
greater value then the actually applied force.
Now I want to compute the reaction force using alternative 2, i.e.
bending the beam with a prescribed displacement. For this I add a
Dirichlet condition on region 2 using the penalization version:
diri = [0 -1.5 0];
gf_model_set(md, 'add initialized data', 'dirichletdata', diri);
gf_model_set(md, 'add Dirichlet condition with penalization', ...
mim, 'u', 1e10, 2);
When I use the multiplier version for both Dirichlet conditions the beam
doesn´t stay in the fixation, hence I use penalization. With my routine
for computing the reaction force I don´t get any sensible results. I think
the reason is that the penalty coefficients change the values in the
tangent matrix. Does anyone have a clue how can I fix this problem? Maybe
I have to change my complete routine or use multipliers/penalization in a
different way.
Thank you very much in advance,
Julia Guenther
PS: My complete current Matlab program:
%% ---------- Mesh ----------
L = 10;
m = gfMesh('cartesian', 0:.5:L, 0:.5:1, 0:.5:1);
%% ---------- FEM ----------
n = gf_mesh_get(m, 'dim');
k = 1; % Degree
% FEM u
f = gf_fem(sprintf('FEM_QK(%d,%d)',n,k));
mfu = gf_mesh_fem(m,3); gf_mesh_fem_set(mfu, 'fem', f);
% FEM Von Mises
f2 = gf_fem(sprintf('FEM_QK_DISCONTINUOUS(%d,1)', n));
mfvm = gf_mesh_fem(m,1); gf_mesh_fem_set(mfvm, 'fem', f2);
%% ---------- IM ----------
k_i = 2*k;
INTEG_TYPE = sprintf('IM_GAUSS_PARALLELEPIPED(%d,%d)',n,k_i);
im = gf_integ(INTEG_TYPE); mim = gf_mesh_im(m,im);
%% ---------- region ----------
P = get(m,'pts');
boundary_left = gf_mesh_get(m, 'faces from pid', find(abs(P(1,:))<1e-6));
boundary_right = gf_mesh_get(m, 'faces from pid', find(abs(P(1,:) -
L)<1e-6));
gf_mesh_set(m, 'region', 1, boundary_left);
gf_mesh_set(m, 'region', 2, boundary_right);
%% ---------- data ----------
E = 200000;
nu = 0.3;
% Lamé
lambda = nu*E/((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
%% ---------- 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');
% homogenous Dirichlet on the left (region 1)
% gf_model_set(md, 'add Dirichlet condition with penalization', ...
% mim, 'u', 1e10, 1);
gf_model_set(md, 'add Dirichlet condition with multipliers', ...
mim, 'u', mfu, 1);
alt = 2; %alternative 1 Force, alternative 2 Displacement
% alternative 1: Force on the right side (region 2)
if alt ==1
F = [0 -100 0];
gf_model_set(md, 'add initialized data', 'Kraft', F);
gf_model_set(md, 'add source term brick', mim, 'u', 'Kraft', 2);
else
% alternative 2: Displacement on the right side (region 2)
diri = [0 -1.5 0];
gf_model_set(md, 'add initialized data', 'dirichletdata', diri);
gf_model_set(md, 'add Dirichlet condition with penalization', ...
mim, 'u', 1e10, 2);
% gf_model_set(md, 'add Dirichlet condition with multipliers', ...
% mim, 'u', mfu, 2, 'dirichletdata');
end
%% ---------- solving ----------
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', mfvm);
figure
gf_plot(mfvm,VM, 'deformation',u, 'deformation_mf',mfu, 'deformed_mesh',
'on', 'cvlst', get(m, 'outer faces'), 'deformation_scale', 1,
'disp_options', 'off');
colorbar;
xlabel('x'); ylabel('y'); zlabel('z');
%% ---------- reaction force ----------
% get the tangent matrix
tangent_matrix = gf_model_get(md, 'tangent_matrix');
% get the multipliers
mult = gf_model_get(md, 'variable', 'mult_on_u');
% get the number of multipliers and DOFs
nb_mult = size(mult,2);
nb_dof = gf_mesh_fem_get(mfu,'nbdof');
% part of the tangent matrix concerning the multipliers
mult_matrix = zeros(nb_dof,nb_mult);
for i = 1:nb_dof
for j = 1:nb_mult
mult_matrix(i,j) = tangent_matrix(nb_mult+i,j);
end
end
% computing the nodal forces by multiplying the multipliers
% with the right part of the tangent matrix
nodalforce = -(mult_matrix*transpose(mult));
% extract the x-, y- and z-components
for i = 0:(size(nodalforce,1)/3)-1
nodalforce_x(i+1) = nodalforce(3*i+1);
nodalforce_y(i+1) = nodalforce(3*i+2);
nodalforce_z(i+1) = nodalforce(3*i+3);
end
% sum the components for the total reaction force
reactionforce_x = sum(nodalforce_x);
reactionforce_y = sum(nodalforce_y);
reactionforce_z = sum(nodalforce_z);
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
