Re: [Getfem-users] Displacement fields in uniaxial tensile test

[Konstantinos Poulios]

Dear Xin Liu

The boundary condition that you have applied on the left side of your block is not the one you are showing in your figure.

You need to use the function

add_normal_Dirichlet_condition_with_multipliers

and somehow also restrict your structure in the vertical direction. You can for example apply the same boundary condition to the bottom side.

Best regards

Kostas

On Fri, May 24, 2019 at 3:29 AM Xin Liu <address@hidden> wrote:

Dear GetFEM-Users,

I wrote a simple 2D uniaxial tensile problem (see the mesh and coordinates in the attachment) mainly based on the example "thermo_elasticity_electrical_coupling".

The code is:

"

#include "getfem/getfem_mesh.h"

#include "getfem/getfem_mesh_fem.h"

#include "getfem/getfem_import.h"

#include "getfem/getfem_generic_assembly.h"

#include "gmm/gmm.h"

#include "getfem/getfem_model_solvers.h"

#include "getfem/getfem_mesher.h"

#include "getfem/getfem_fem.h"

using std::endl; using std::cout;

using bgeot::base_small_vector;

typedef getfem::model_real_plain_vector plain_vector;

int main(){

getfem::mesh m;

getfem::import_mesh("gmsh:./Simple.msh",m);

getfem::mesh_region border_faces;

getfem::outer_faces_of_mesh(m, border_faces);

getfem::mesh_region fb1 = getfem::select_faces_of_normal(m, border_faces,

base_small_vector( 1., 0.), 0.01);

getfem::mesh_region fb2 = getfem::select_faces_of_normal(m, border_faces,

base_small_vector( -1., 0.), 0.01);

m.region(1) = fb1; m.region(2) = fb2;

// Set finite element methods

getfem::pfem pf = getfem::fem_descriptor("FEM_QK(2,1)");

getfem::mesh_fem mfu(m); // Finite element for the elastic displacement

mfu.set_qdim(2);

mfu.set_finite_element(pf);

getfem::mesh_im mim(m); // Integration method

mim.set_integration_method(bgeot::dim_type(2));

getfem::model md;

md.add_fem_variable("u", mfu); // Displacement of the structure

// Membrane elastic deformation

double E = 21E6;

double nu = 0.3;

double F = 100E2;

double clambda = E*nu/((1+nu)*(1-2*nu)); // First Lame coefficient (N/cm^2)

double cmu = E/(2*(1+nu)); // Second Lame coefficient (N/cm^2)

double clambdastar = 2*clambda*cmu/(clambda+2*cmu); // Lame coefficient

md.add_initialized_scalar_data("cmu", cmu);

md.add_initialized_scalar_data("clambdastar", clambdastar);

getfem::add_isotropic_linearized_elasticity_brick

(md, mim, "u", "clambdastar", "cmu");

getfem::add_Dirichlet_condition_with_multipliers

(md, mim, "u", bgeot::dim_type(2-1), 2);

md.add_initialized_fixed_size_data("Fdata", base_small_vector(F,0.));

getfem::add_source_term_brick(md, mim, "u", "Fdata", 1);

gmm::iteration iter(1E-9, 1, 100);

getfem::standard_solve(md, iter);

//

// Solution export

plain_vector U(mfu.nb_dof());

gmm::copy(md.real_variable("u"), U);

bgeot::base_vector coeff;

getfem::slice_vector_on_basic_dof_of_element(mfu,U,1,coeff);

std::vector<std::string> v = { "Point a ", "Point b ", "Point c ", "Point d "};

cout << "Displacements for element 1: \n";

for (int i=0; i < coeff.size()/2; i++)

cout << v[i] << "u1: " << coeff[2*i] << "; u2： " << coeff[2*i+1] << endl;

return 0;

}

"

I output the displacement at points a, b, c, and d (see attachment).
Point a  u1: 0.000145356;  u2： 1.79962e-05
Point b  u1: 0.000161904;  u2： 6.42766e-05
Point c  u1: 0.000306614;  u2： 2.5169e-05
Point d  u1: 0.000313698;  u2： 7.23111e-05

I expect the u1 value of a and b is the same and c and d is the same in this problem since this is just an uniaxial tensile test. So does the difference come from numerical errors or some mistakes in the code?

Any comments/suggestions will be greatly appreicated.

Best regards,

--

======================================

Xin Liu

PhD Candidate, Graduate Research Assistant

Multiscale Structural Mechanics Group

https://sites.google.com/view/xinliu8809/

School of Aeronautics & Astronautics

Purdue University

West Lafayette, Indiana 47907-2045

======================================