[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Getfem-users] problem with structural calculations in getfem
From: |
Yves Renard |
Subject: |
Re: [Getfem-users] problem with structural calculations in getfem |
Date: |
Wed, 7 Nov 2007 13:39:09 +0100 |
User-agent: |
KMail/1.9.5 |
On Tuesday 06 November 2007 17:50, Andriy Andreykiv wrote:
> Dear Yves,
>
> First of all thank you very much for debugging that two NonLin term
> problem. If I also may consult you regarding the following:
>
> While looking for bugs in my electrostatic-structural implementation I
> confrontated a problem with quite simple structural calculation in getfem.
> Although this calculation gives a result that looks nice, comparing the
> numbers whith commercial FEM package MSC. Marc gave very different
> solution. The problem is very simple (please see the picture in the
> attachment):
>
> A rectange with a size 0.8 by 0.1 is clamped at the sides.
> A constant distributed load P is applied on the bottom surface of this
> rectangle. In getfem, this load (pressure) is prescribed by source_term,
> with vector "V(#1)+=comp(Normal().vBase(#1))(i,:,i);", assembled on the
> bottom boundary. Then this vector is scaled with the actual value of the
> pressure. In both Marc and Getfem the corresponding nodal force is
> calculated on the underformed geometry of the body (hence, no follower
> force). The force causes the rectangle to subside downwards and the
> deflection of the middle part of the rectangle is measured.
>
> Elastic properties (used by MSC. Marc) are
> Young moduls: 20
> Poisson ratio: 0.3
> Plain strain assumed for 2D problem
>
> Corresponding Lame constants (used in Getfem program):
> Lambda: 11.5384
> Mu: 7.6923
>
> In both programs I used linear quad elements, splitting the rectangle 32X4.
> Both programs use 4 Gauss points for element integration. A similar problem
> was also solved in 3D for 0.8X0.1X0.1 paralellepiped structure, descretized
> by 32X4X4 linear hex elements. Bellow you can see small and finite strain
> results from the two programs:
>
> 2D, finite strain, distributed load P=-2.0:
> MSC. Marc deflection = 0.229; Getfem deflection 0.192651
>
> 2D, small strain, distributed load P=-0.4
> MSC.Marc deflection = 0.2667; Getfem deflection 0.1699
>
> 3D, finite strain, distributed load P=-2.0
> MSC. Marc deflection = 0.2375; Getfem deflection 0.0191458 !!!
>
> 3D, small strain, distributed load P=-0.4
> MSC. Marc deflection = 0.287; Getfem delflection 0.00400608 !!!
>
> Even if Marc results are not correct, the agreement between 2D and 3D
> solution is much better then in Getfem. Additionally, when doing either p-
> or h- refinement in Marc, the result changes only slightly, by increasing
> the deflection (2-3 %). When I refine the mesh in Getfem, the deflection
> considerably reduces which is really strange.
>
> The source of the program with the corresponding four parameter files
> (2D/3D and finite/small strain cases) can be found in the attachment.
>
> Thank you in advance,
> Andriy
>
>
>
>
> __________________________________________
> Andriy Andreykiv (PhD, MSc)
> Delft University of Technology
> Faculty of Mechanical Engineering
> Material Science and Marine Technology
> Group: Fundamentals of Microsystems
> Mekelweg 2
> 2628 CD Delft
> The Netherlands
>
> E-mail : address@hidden
> Tel. : (31) 15 2786818
> Fax. : (31) 15 2789475
> www : http://www-tm.wbmt.tudelft.nl/~andrico
> private: (31) 6 47376804
Dear Andriy,
I see at least one error in your program. When you run the following lines
plain_vector Unit_Pressure(mf_mechanic.nb_dof());
getfem::generic_assembly assem_boundary;
assem_boundary.set("V(#1)+=comp(Normal().vBase(#1))(i,:,i);");
assem_boundary.push_mi(mim_mechanic);
assem_boundary.push_mf(mf_mechanic);
assem_boundary.push_vec(Unit_Pressure);
assem_boundary.assembly(m_mechanic.region(SOUTH));
plain_vector Current_Pressure(Unit_Pressure);
gmm::scale(Current_Pressure,p_current);
pstruct_brick->source_term().set(Current_Pressure);
The vector you give to pstruct_brick->source_term().set(.) is already
assembled. But, it will be assembled a second time ! You dont have to give an
assembled vector to this function. You have to give a vector of degree of
freedom corresponding to your mf_data. For instance a more correct code will
be
plain_vector Current_Pressure(data_size*m_mechanic.dim());
gmm::fill(gmm::sub_vector(Current_Pressure,
gmm::sub_slice(1, data_size, 2)),
p_current);
pstruct_brick->source_term().set(Current_Pressure);
(the fact that the vector is defined outside the considered zone is not
important).
Yves.
--
Yves Renard (address@hidden) tel : (33) 04.72.43.87.08
Pole de Mathematiques, fax : (33) 04.72.43.85.29
INSA de Lyon, Universite de Lyon
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------