
From:  Konstantinos Poulios 
Subject:  Re: [Getfemusers] Antwort: Re: 2D nonlinear plane stress 
Date:  Thu, 4 Jul 2019 21:39:38 +0200 
Hello Kostas,
your approach works perfectly for me.
I checked me code and it seems like there was a little typing error. So what I did at first works too. But your approach seems to be faster. (Some warnings occure in my code.)
I have one last question at the moment. This question concerns the line
md.add_nonlinear_term(mim, "Th*0.5*kappa*sqr(log(Det(F)))+"
Th*0.5*mu*(pow(Det(F),2/3)*Norm_sqr(F)3)")
Where did you get that _expression_ from? Is there any literatur, where I can find such expressions? Or can you give me an approach to deriving that _expression_?
Best regards,
Moritz Jena
Von: Konstantinos Poulios <address@hidden>
An: Moritz Jena <address@hidden>
Kopie: getfemusers <address@hidden>
Datum: 03.07.2019 23:20
Betreff: Re: Re: [Getfemusers] Antwort: Re: 2D nonlinear plane stress
Dear Jena,
Without having run your code I would say that your definition of the 3D deformation gradient looks correct. Another way of converting from 2x2 to 3x3 is
"(Id(3)+[[1,0,0],[0,1,0]]*Grad_u*[[1,0,0],[0,1,0]]')"
but it should be equivalent to your approach.
Could you attach a screenshot of the "strange" solution that you get?
BR
Kostas
On Wed, Jul 3, 2019 at 9:49 AM Moritz Jena <address@hidden> wrote:
Hello Kostas and all GetFEMUsers,
I finally managed to try your approach for my plane stress problem.
But unfortunately I got a problem with it.
When defining the 3D deformation gradient, I have to add a [2x2]  matrix to a [3x3] matrix. To work around this problem I tried to rewrite the displacement gradient to a [3x3]  matrix.
So instead of
md.add_macro("F", "Id(3)+Grad_u+epsZ*[0,0,0;0,0,0;0,0,1]")
it tried
md.add_macro("F", "Id(3)+[Grad_u(1,1),Grad_u(1,2),0;Grad_u(2,1),Grad_u(2,2),0;0,0,0]+epsz*[0,0,0;0,0,0;0,0,1]")
To test the approach, I created a little testscript with a simple beam and a force. The calculation works with this work around, but the deformation seems something strange. I think I did a mistake by accessing the matrixelements of the displacement gradient, but I'm not sure about it.
Can you, or anyone else, help me with my problem?
import getfem as gf
import numpy as np
# Parameter
l = 100
h = 10
b = 1.5
size = 1
E = 203000
nu = 0.3
Lambda = (E*nu)/(1pow(nu,2))
mu = E/(1+nu)
lawname = 'SaintVenant Kirchhoff'
# Create mesh
m = gf.Mesh('cartesian', np.arange(0, l+size , size), np.arange(0, h+size, size))
#MeshFem
mfu = gf.MeshFem(m,2)
mfd = gf.MeshFem(m,1)
mf1 = gf.MeshFem(m,1)
#assign FEM
mfu.set_fem(gf.Fem('FEM_QK(2,2)'))
mfd.set_fem(gf.Fem('FEM_QK(2,2)'))
mf1.set_fem(gf.Fem('FEM_QK(2,1)'))
#IM_methode
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS_PARALLELEPIPED(2,6)'))
# Some Infos
print("nbcvs = %d  nbpts = %d  qdim = %d  nbdofs = %d" % (m.nbcvs(), m.nbpts(), mfu.qdim(), mfu.nbdof()))
#regions
P = m.pts();
pidleft = np.compress((abs(P[0,:])<1e6),range(0,m.nbpts()))
pidrigth = np.compress((abs(P[0,:])>l1e6),range(0, m.nbpts()))
left = m.outer_faces_with_direction([1, 0], 0.5)
rigth = m.outer_faces_with_direction([1,0], 0.5)
fleft = m.faces_from_pid(pidleft)
frigth = m.faces_from_pid(pidrigth)
m.set_region(1, fleft)
m.set_region(2, frigth)
#Model
md = gf.Model('real')
md.add_fem_variable('u', mfu)
md.add_initialized_data('lambda', Lambda)
md.add_initialized_data('mu', mu)
md.add_fem_variable("epsz", mf1)
md.add_initialized_data('kappa', 68000)
md.add_initialized_data('Th', 0.2)
#material
md.add_macro("F", "Id(3)+[Grad_u(1,1),Grad_u(1,2),0;Grad_u(2,1),Grad_u(2,2),0;0,0,0]+epsz*[0,0,0;0,0,0;0,0,1]")
md.add_nonlinear_term(mim, "Th*0.5*kappa*sqr(log(Det(F)))+Th*0.5*mu*(pow(Det(F),2/3)*Norm_sqr(F)3)")
# fixed support left side
md.add_initialized_data('DirichletData1', [0,0])
md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu, 1, 'DirichletData1')
# force on rigth side
md.add_initialized_data('VolumicData', [0,0])
md.add_source_term_brick(mim, 'u', 'VolumicData')
md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu, 2, 'VolumicData')
# force array
F = np.array([[0,4], [0,10], [0,15], [0,25], [0,35], [0,50]])
nbstep = F.shape[0]
#iterative calc
for step in range(0, nbstep):
print(step)
md.set_variable('VolumicData', [F[step, 0],F[step,1]])
md.solve('noisy', 'max_iter', 50)
U = md.variable('u')
s1 = gf.Slice(('boundary',), mfu, 4)
s1.export_to_vtk('test_%d.vtk' % step, 'ascii', mfu, U , 'Displacement')
Von: Konstantinos Poulios <address@hidden>
An: Yves Renard <address@hidden>
Kopie: Moritz Jena <address@hidden>, getfemusers <address@hidden>
Datum: 05.02.2019 10:05
Betreff: Re: [Getfemusers] Antwort: Re: 2D nonlinear plane stress
Dear Moritz Jena,
To define a hyperelastic material law for a plane stress problem, the easiest way is to add one extra variable representing the out of plane strain. If mf1 is a scalar fem and mf2 is a vector fem with 2 components, then you can simply do:
md = gf.Model("real")
md.add_fem_variable("u", mf2) # displacements variable
md.add_fem_variable("epsZ", mf1) # out of plane strain variable
md.add_initialized_data(’kappa’, kappa) # initial bulk modulus
md.add_initialized_data(’mu’, mu) # initial shear modulus
md.add_initialized_data(’Th’, Th) # plate thickness
md.add_macro("F", "Id(3)+Grad_u+epsZ*[0,0,0;0,0,0;0,0,1]") # 3D deformation gradient
md.add_nonlinear_term(mim, "Th*0.5*kappa*sqr(log(Det(F)))+"
Th*0.5*mu*(pow(Det(F),2/3)*Norm_sqr(F)3)")
Best regards
Kostas
On Wed, Jan 16, 2019 at 9:40 PM Yves Renard <address@hidden> wrote:
Dear Jena,
For an hyperelastic law, the weak form of the static elastic problem can be written in the weak form language
"Def F := Id(meshdim)+Grad_u; (F * S) : Grad_test_u"
for u the displacement, F the deformation gradient and S has to contains the _expression_ of the second PiolaKirchhoff stress tensor (you can of course express it in term of PK1 also). For instance for a St Venant Kirchhoff law, you can write
"Def F := Id(meshdim)+Grad_u; Def E := 0.5*(F'*FId(meshdim)); (F * (lambda*Trace(E)+2*mu*E)) : Grad_test_u"
where E will be the Green Lagrange deformation tensor and lambda, mu the Lamé coefficients.
This gives you some examples of construction of hyperelastic laws. The weak form language gives you access to some standard operators (Trace, Det ...) see http://getfem.org/userdoc/model_nonlinear_elasticity.html#highlevelgenericassemblyversions
So, if you have the _expression_ of your law in plane stress, it should not be very difficult to implement it. But of course you need the _expression_ of the law in plane stress. On the construction itself of plane stress hyperelastic law, now, I don't know a good reference, unfortunately.
Best regards,
Yves
 Original Message 
From: "Moritz Jena" <address@hidden>
To: "yves renard" <address@hidden>
Cc: "getfemusers" <address@hidden>
Sent: Wednesday, January 16, 2019 2:17:01 PM
Subject: Antwort: Re: [Getfemusers] 2D nonlinear plane stress
Hello Yves,
thank you for your answer.
I'm afraid I'm not into the topic weak form language and I'm not sure
where to start with this problem.
I looked at the chapter in the documentation, however I don't know how to
describe a plane stress material model with it.
I also studied the examples, that come with the MATLABInterface. There
are a few examples, how to declare a material model with this weak form
expressions. But I still don't know, how to build this expressions.
Can you give me a approach for this problem or where I can find
expressions for such a problem? Is there any literature that you can
recommend?
Best regards,
Moritz
Von: Yves Renard <address@hidden>
An: Moritz Jena <address@hidden>
Kopie: getfemusers <address@hidden>
Datum: 09.01.2019 17:17
Betreff: Re: [Getfemusers] 2D nonlinear plane stress
Dear Moritz Jena,
No, unfortunately, plane stress versions of Hyperelastic laws has not been
implemented yet.
I would not be so difficult, but as to be made. If you need one in
particular and have the _expression_, it is no so difficult to describe it
with the weak form language.
Best regards,
Yves
 Original Message 
From: "Moritz Jena" <address@hidden>
To: "getfemusers" <address@hidden>
Sent: Tuesday, January 8, 2019 4:11:48 PM
Subject: [Getfemusers] 2D nonlinear plane stress
Dear GetFEMUsers,
I use the MATLABInterface of GetFEM to create a program that
automatically solves different models of the same problem.
The problem is threedimensional, but can be reduced by plain stress
approximation. (to reduce computing time).
I want to define a nonlinear material with the brick
gf_model_set(model M, 'add nonlinear elasticity brick',
[...])
For this nonlinear command it is specified in the description, that in 2D
always plain strain is used.
So my question is: Is there a way to define a nonlinear material with the
plain stress approximation? Or is it planned to install such an option in
a future release?
I hope you can help me with my problem.
Best regards,
Moritz Jena
neohookean.pdf
Description: Adobe PDF document
[Prev in Thread]  Current Thread  [Next in Thread] 