[Getfem-users] reaction force for prescribed displacements

[Julia Guenther]

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);