|
From: | Julia Guenther |
Subject: | [Getfem-users] compute a reaction force |
Date: | Tue, 17 Dec 2013 15:42:58 +0100 |
I’m using the Matlab interface of GetFEM and I tried to build a very simple three-dimensional beam. It is fixed on the left side. On the right side I have two alternatives to bend it:
1. A force
2. A prescribed displacement
My question is: How can I compute the reaction force on the left side, when I use alternative 2?
I´ve found something similar in an earlier post ( http://www.mail-archive.com/address@hidden ), but I´m not able to adapt it for my problem using the Matlab interface.
I added my matlab-code below. I’m pretty new to GetFEM and FEM at all, hence I´m not sure if I set all the conditions correctly, especially concerning the use of multipliers and penalization.
Thank you in advance,
Julia Guenther
%% -------------------------------
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);
%% ------------------------------
Regions --------------------------------
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 region 1
gf_model_set(md, 'add
Dirichlet condition with penalization',
...
mim, 'u',
1e10, 1);
alt = 2; %
1 Force, 2 Displacement
% alternative 1: Force
on region 2
if
alt ==1
F = [0 -100 0];
gf_model_set(md, 'add
initialized data', 'Force',
F);
gf_model_set(md, 'add
source term brick', mim, 'u',
'Force',
2);
else
% alternative 2: Displacement
on region 2
diri = [0 -1.5 0];
gf_model_set(md, 'add
initialized data', 'dirichletdata',
diri);
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);
f = figure;
gf_plot(mfvm,VM, 'deformation',u,
'deformation_mf',mfu,
'deformed_mesh',
'on',
'cvlst',
get(m, 'outer faces'),
'deformation_scale',
1);
colorbar;
xlabel('x')
ylabel('y')
zlabel('z')
set(f, 'Units',
'normalized',
'Position',
[0.2, 0.1, 0.7, 0.7]);
[Prev in Thread] | Current Thread | [Next in Thread] |