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 <[email protected]> wrote: > Hello Kostas and all GetFEM-Users, > > 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 test-script 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 matrix-elements 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)/(1-pow(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,:])<1e-6),range(0,m.nbpts())) > pidrigth = np.compress((abs(P[0,:])>l-1e-6),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 <[email protected]> > An: Yves Renard <[email protected]> > Kopie: Moritz Jena <[email protected]>, getfem-users < > [email protected]> > Datum: 05.02.2019 10:05 > Betreff: Re: [Getfem-users] 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 <*[email protected]* > <[email protected]>> 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 Piola-Kirchhoff 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'*F-Id(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#high-level-generic-assembly-versions* > <http://getfem.org/userdoc/model_nonlinear_elasticity.html#high-level-generic-assembly-versions> > > > 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" <*[email protected]* <[email protected]> > > > To: "yves renard" <*[email protected]* <[email protected]>> > Cc: "getfem-users" <*[email protected]* <[email protected]>> > Sent: Wednesday, January 16, 2019 2:17:01 PM > Subject: Antwort: Re: [Getfem-users] 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 MATLAB-Interface. 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 <*[email protected]* <[email protected]> > > > An: Moritz Jena <*[email protected]* <[email protected]> > > > Kopie: getfem-users <*[email protected]* <[email protected]>> > Datum: 09.01.2019 17:17 > Betreff: Re: [Getfem-users] 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" <*[email protected]* <[email protected]> > > > To: "getfem-users" <*[email protected]* <[email protected]>> > Sent: Tuesday, January 8, 2019 4:11:48 PM > Subject: [Getfem-users] 2D nonlinear plane stress > > Dear GetFEM-Users, > > I use the MATLAB-Interface of GetFEM to create a program that > automatically solves different models of the same problem. > > The problem is three-dimensional, 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 > > >
