Dear Raphaël,
After Yves has now fixed this bug
<https://git.savannah.nongnu.org/cgit/getfem.git/commit/?id=57e5e2d431c76987900528a1b7f7badd0d0d3893>,
I came back to your case. Before trying to solve your original question, I
have prepared a version of the demo_wheel.py script (attached) where you
can simulate the contact with friction between the two bodies for
prescribed vertical displacement and rotation. Please recompile GetFEM with
Yves' fix and then try to run the attached script to make sure it works.
As an exercise you can then try to add the calculation of the moment Mz,
and the reaction forces Fx and Fy to the attached script.
Best regards
Kostas
On Tue, Apr 19, 2022 at 3:18 PM Raphaël Meunier <raphael.meuni...@gmail.com>
wrote:
> Dear Kostas,
>
> Thank you again for the kind answer.
>
> I tried the given script and it outputs an error message:
> [image: image.png]
> Which does not help me towards finding the right answer.
> So, I tried to describe the settings of my experiment in the attached
> document. Please, feel free to ask if it is not clear enough.
>
> Thank you in advance,
>
> Best regards,
>
> Raphaël Meunier
>
> Le sam. 16 avr. 2022 à 19:09, Konstantinos Poulios <
> logar...@googlemail.com> a écrit :
>
>> Dear Raphaël,
>>
>> It would help to have a description/sketch of the exact system you would
>> like to model. But if I understood it correctly, you wanted to prescribe a
>> constant force Fy in y-direction. If you introduce a degree of freedom yc
>> for the vertical displacement, and you couple your wheel to yc, then you
>> can easily apply a force on yc by multiplying the known force with the
>> variation of yc. Dividing with area is needed only if you add this term
>> inside an integral over an area "area". This is because the integral of a
>> constant is just as multiplying the constant with the area of the
>> integration domain. You could achieve the same without an integral using
>> the "model.add_explicit_rhs" function on the yc variable with rhs value
>> equal to Fy (or -Fy, I do not remember exactly).
>>
>> Regarding the term
>> "([Fx/area,0]*rot_der(phi))*(X-[xc;yc])*Test_phi"
>> I am not 100% sure about. This is the part based on intuition and some
>> partial understanding of what you want to achieve. You can give it a try,
>> but if it doesn't work please provide a proper sketch of your system so
>> that I can help you with deriving this equilibrium condition correctly.
>>
>> BR
>> Kostas
>>
>>
>>
>> On Fri, Apr 15, 2022 at 1:51 PM Raphaël Meunier <
>> raphael.meuni...@gmail.com> wrote:
>>
>>> Dear Kostas,
>>>
>>>
>>>
>>> Thank you for your response. As my intuition is not as sharp as yours, I
>>> might need a few explanations:
>>>
>>> - md.add_nonlinear_term(mim, "(u-rot(phi)*(X-[xc;yc])).mult",
>>> HoleZone) --> This order zero term seems pretty straightforward and
>>> much more readable than in my code
>>> - md.add_linear_term(mim, "Fy/area*Test_yc", HoleZone) --> No
>>> multiplier is involved in your formulation. Is this because we need to
>>> combine this term with the previous order 0 term
>>> -((rot(phi)*([0;1])).mult)*Test_yc (and check the sign) ?
>>> - md.add_linear_term(mim,
>>> "([Fx/area,0]*rot_der(phi))*(X-[xc;yc])*Test_phi", HoleZone) -->
>>> Could you explain how do you come up with this term ?
>>> - A more general question: Is it possible to formulate the whole
>>> condition as an order 0 term ?
>>>
>>>
>>> Thanks in advance,
>>>
>>> Best regards
>>>
>>> Raphaël
>>>
>>> Le mer. 6 avr. 2022 à 23:14, Konstantinos Poulios <
>>> logar...@googlemail.com> a écrit :
>>>
>>>> Dear Raphaël,
>>>>
>>>> Interesting question. Prescribing Fx is a bit strange condition. I
>>>> haven't checked it, but I think something like this should work:
>>>>
>>>> md.add_initialized_data("Fx", 10.)
>>>> md.add_initialized_data("Fy", 100.)
>>>>
>>>> md.add_initialized_data("xc", 0.)
>>>> md.add_variable("yc", 1)
>>>> md.set_variable("yc", yc)
>>>> md.add_variable("phi",1)
>>>> md.set_variable("phi", 0.)
>>>> md.add_fem_variable("u", mfu)
>>>> md.add_filtered_fem_variable("mult", mfu, HoleZone)
>>>>
>>>> md.add_macro("rot(x)", "[cos(x)-1,sin(x);-sin(x),cos(x)-1]")
>>>> md.add_macro("rot_der(x)", "[-sin(x),cos(x);-cos(x),-sin(x)]")
>>>>
>>>> md.add_initialized_data("area", gf.asm_generic(mim, 0, "1", HoleZone))
>>>> # Calculate the area of the region
>>>> md.add_nonlinear_term(mim, "(u-rot(phi)*(X-[xc;yc])).mult", HoleZone)
>>>> md.add_linear_term(mim, "Fy/area*Test_yc", HoleZone)
>>>> # check if a minus sign is needed
>>>> md.add_linear_term(mim,
>>>> "([Fx/area,0]*rot_der(phi))*(X-[xc;yc])*Test_phi", HoleZone) # check if a
>>>> minus sign is needed
>>>>
>>>> for time in [0,0.1,0.3,0.6,0.9]:
>>>> md.set_variable("xc", speed*time)
>>>> md.solve("noisy")
>>>>
>>>> A bit based on intuition but you can certainly learn something from
>>>> analyzing this suggestion. Let us discuss it further if you cannot get it
>>>> to work and have more questions. Try to avoid expressions like "speed*time"
>>>> and "alpha*time" in your expressions, they seem too specific to constant
>>>> velocities, they make the code more difficult to read, and they will not
>>>> help you to generalize the code later. Just solve for positions/angles
>>>> directly if that's what you want.
>>>>
>>>> FYI, the expression
>>>> (u-rot(phi)*(X-[xc;yc])).mult
>>>> will expand internally to something like
>>>> (u-rot(phi)*(X-[xc;yc])).Test_mult
>>>> +Test_u.mult
>>>> -((rot_der(phi)*(X-[xc;yc])).mult)*Test_phi
>>>> -((rot(phi)*([0;1])).mult)*Test_yc
>>>>
>>>> BR
>>>> Kostas
>>>>
>>>>
>>>> On Wed, Apr 6, 2022 at 4:54 PM Raphaël Meunier <
>>>> raphael.meuni...@gmail.com> wrote:
>>>>
>>>>> Dear Getfem-Users,
>>>>>
>>>>> I am trying to impose a rolling condition on a "tire". My work is
>>>>> directly inspired from the "Example of wheel in contact" tutorial
>>>>> <https://getfem.org/tutorial/wheel.html>. In this tutorial, one wants
>>>>> to prescribe the rim rigidity and the vertical force. From this, I just
>>>>> have to add a rolling condition. Specifically I want to find the rotation
>>>>> speed "alpha_rot" that enables me to impose some force Fx:
>>>>>
>>>>> gwfl_rot='-lambda_D.Test_u +
>>>>> ([cos(alpha_rot*time),sin(alpha_rot*time);-sin(alpha_rot*time),cos(alpha_rot*time)]*(X-[0;y_c])-(X-[0;y_c])\
>>>>> + [speed*time;alpha_y] - u).Test_lambda_D +(lambda_D.[1;0] +
>>>>> Fx)*Test_alpha_rot + 1E-6*alpha_rot*Test_alpha_rot'
>>>>> gwfl_y='(lambda_D.[0;1] + Fy)*Test_alpha_y + 1E-6*alpha_y*Test_alpha_y'
>>>>> roll_gwfl=gwfl_y+gwfl_rot
>>>>>
>>>>> idBrick = model.add_nonlinear_term(mim, roll_gwfl, HoleZone)
>>>>>
>>>>>
>>>>> Where alpha_rot and alpha_y are the "parameters" I am trying to find.
>>>>> This works for up to ten timesteps but then it diverges... I tried to use
>>>>> two multipliers instead of one which leads to worse results.
>>>>> Am I on the right track ? Is there a better way to achieve this ?
>>>>>
>>>>> Thanks in advance,
>>>>>
>>>>> Best regards,
>>>>>
>>>>> Raphaël
>>>>>
>>>>
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Python GetFEM interface
#
# Copyright (C) 2015-2020 Yves Renard.
#
# This file is a part of GetFEM
#
# GetFEM is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version along with the GCC Runtime Library
# Exception either version 3.1 or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
# License and GCC Runtime Library Exception for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
#
# Contact of a deformable 'wheel' onto a plane deformable obstacle.
#
############################################################################
import numpy as np
import getfem as gf
gf.util_trace_level(1) # No trace for mesh generation nor for assembly
#
# Physical parameters
#
E = 21E6 # Young Modulus (N/cm^2)
nu = 0.3 # Poisson ratio
clambda = E*nu/((1+nu)*(1-2*nu)) # First Lame coefficient (N/cm^2)
cmu = E/(2*(1+nu)) # Second Lame coefficient (N/cm^2)
applied_force = 1E7 # Force at the hole boundary (N)
#
# Numerical parameters
#
h = 1 # Approximate mesh size
elements_degree = 2 # Degree of the finite element methods
r_aug = E; # Augmentation parameter for the augmented Lagrangian
#
# Mesh generation. Meshes can also been imported from several formats.
R_o = 15.
R_i = 8.
mo1 = gf.MesherObject('ball', [0., R_o], R_o)
mo2 = gf.MesherObject('ball', [0., R_o], R_i)
mo3 = gf.MesherObject('set minus', mo1, mo2)
print('Meshes generation')
mesh1 = gf.Mesh('generate', mo3, h, 2)
#mesh1.translate([0.,-0.1])
mesh2 = gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[30,10];NOISED=0;NSUBDIV=[%d,%d];' % (int(30/h)+1, int(10/h)+1));
mesh2.translate([-15.,-10.])
# Boundary selection
#
fb1 = mesh1.outer_faces_in_box([-8.1, 6.9], [8.1, 23.1]) # Boundary of the hole
fb2 = mesh1.outer_faces_with_direction([0., -1.], np.pi/4.5) # Contact boundary of the wheel
fb3 = mesh2.outer_faces_with_direction([0., -1.], 0.01) # Bottom boundary of the foundation
fb4 = mesh2.outer_faces_with_direction([0., 1.], 0.01) # Top boundary of the foundation
HOLE_BOUND=1; CONTACT_BOUND=2; BOTTOM_BOUND=3; TOP_BOUND=4;
mesh1.set_region(HOLE_BOUND, fb1)
mesh1.set_region(CONTACT_BOUND, fb2)
mesh1.region_subtract(CONTACT_BOUND, HOLE_BOUND)
mesh2.set_region(BOTTOM_BOUND, fb3)
mesh2.set_region(TOP_BOUND, fb4)
print(mesh2.region(TOP_BOUND))
#
# Definition of finite elements methods and integration method
#
mfu1 = gf.MeshFem(mesh1, 2)
mfu1.set_classical_fem(elements_degree)
mflambda = gf.MeshFem(mesh1, 2)
mflambda.set_classical_fem(elements_degree-1)
mflambda_C = gf.MeshFem(mesh1, 2)
mflambda_C.set_classical_fem(elements_degree-1)
mfu2 = gf.MeshFem(mesh2, 2)
mfu2.set_classical_fem(elements_degree)
mfvm1 = gf.MeshFem(mesh1, 1)
mfvm1.set_classical_discontinuous_fem(elements_degree)
mfvm2 = gf.MeshFem(mesh2, 1)
mfvm2.set_classical_discontinuous_fem(elements_degree)
mim1 = gf.MeshIm(mesh1, 4)
mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)'))
mim2 = gf.MeshIm(mesh2, 4)
#
# Model definition
#
md=gf.Model('real');
md.add_fem_variable('u1', mfu1) # Displacement of the structure 1
md.add_fem_variable('u2', mfu2) # Displacement of the structure 2
#md.add_initialized_data('cmu', [cmu])
#md.add_initialized_data('clambda', [clambda])
#md.add_isotropic_linearized_elasticity_brick(mim1, 'u1', 'clambda', 'cmu')
#md.add_isotropic_linearized_elasticity_brick(mim2, 'u2', 'clambda', 'cmu')
md.add_initialized_data('matparams', [clambda,cmu])
md.add_finite_strain_elasticity_brick(mim1, 'Compressible_Neo_Hookean', 'u1', 'matparams')
md.add_finite_strain_elasticity_brick(mim2, 'Compressible_Neo_Hookean', 'u2', 'matparams')
md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', elements_degree-1, BOTTOM_BOUND)
# Contact condition (augmented Lagrangian)
md.add_initialized_data('r', [1e0*r_aug])
md.add_initialized_data('f_coeff', [0.5])
md.add_filtered_fem_variable('lambda1', mflambda_C, CONTACT_BOUND)
release_dist = 1
contact_implementation = 0
if contact_implementation == 0:
md.add_raytracing_transformation("contact_trans", release_dist)
md.define_variable_group("u", "u1", "u2")
md.add_slave_contact_boundary_to_raytracing_transformation("contact_trans", mesh1, "u", CONTACT_BOUND)
md.add_master_contact_boundary_to_raytracing_transformation("contact_trans", mesh2, "u", TOP_BOUND)
md.add_nonlinear_term(mim1c, "-lambda1.Test_u1", CONTACT_BOUND)
md.add_nonlinear_term(mim1c, "Interpolate_filter(contact_trans, lambda1.Interpolate(Test_u,contact_trans), 1)", CONTACT_BOUND)
md.add_nonlinear_term(mim1c, "-(1/r)*lambda1.Test_lambda1", CONTACT_BOUND)
md.add_nonlinear_term(mim1c,
"Interpolate_filter(contact_trans,"
"(1/r)*Coulomb_friction_coupled_projection(lambda1,"
"Transformed_unit_vector(Grad_u1, Normal),"
"u1-Interpolate(u,contact_trans),"
"(Interpolate(X,contact_trans)"
"+Interpolate(u,contact_trans)-X-u1)"
".Transformed_unit_vector(Grad_u1, Normal),"
"f_coeff,"
"r).Test_lambda1, 1)", CONTACT_BOUND)
elif contact_implementation == 1:
ind = md.add_integral_large_sliding_contact_brick_raytracing("r", release_dist, "f_coeff", "1", 0)
md.add_slave_contact_boundary_to_large_sliding_contact_brick(ind, mim1c, CONTACT_BOUND, "u1", "lambda1", "")
md.add_master_contact_boundary_to_large_sliding_contact_brick(ind, mim2, TOP_BOUND, "u2", "")
else:
md.add_integral_contact_between_nonmatching_meshes_brick\
(mim1c, 'u1', 'u2', 'lambda1', 'r', 'f_coeff', CONTACT_BOUND, TOP_BOUND)
md.add_initialized_data('xc', 0.)
md.add_initialized_data('yc', R_o)
md.add_variable('dy', 1)
md.set_variable('dy', -0.3)
md.add_variable('phi',1)
md.set_variable('phi', 0.01)
md.disable_variable('dy')
md.disable_variable('phi')
md.add_initialized_data('area', gf.asm_generic(mim1, 0, '1', HOLE_BOUND)) # Calculate the area of the region
md.add_initialized_data('Fy', -applied_force)
md.add_macro("rot(x)", "[cos(x)-1,sin(x);-sin(x),cos(x)-1]")
md.add_macro("rot_der(x)", "[-sin(x),cos(x);-cos(x),-sin(x)]")
md.add_linear_term(mim1, 'Fy/area*Test_dy', HOLE_BOUND)
md.add_initialized_data('Mz', 200)
md.add_filtered_fem_variable('mult', mflambda, HOLE_BOUND)
md.add_nonlinear_term(mim1, '(u1-[0;dy]-rot(phi)*(X-[xc;yc])).Test_mult', HOLE_BOUND)
md.add_nonlinear_term(mim1, 'Test_u1.mult', HOLE_BOUND)
#
# Model solve
#
print('Solve problem with ', md.nbdof(), ' dofs')
md.solve('max_res', 1E-11, 'max_iter', 40, 'noisy' , 'lsearch', 'simplest', 'alpha min', 0.8)
# print('Contact multiplier ', md.variable('lambda1'))
#
# Solution export
#
U1 = md.variable('u1')
U2 = md.variable('u2')
#VM1 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u1', 'clambdastar', 'cmu', mfvm1)
#VM2 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u2', 'clambdastar', 'cmu', mfvm2)
VM1 = md.compute_finite_strain_elasticity_Von_Mises('Compressible_Neo_Hookean', 'u1', 'matparams', mfvm1)
VM2 = md.compute_finite_strain_elasticity_Von_Mises('Compressible_Neo_Hookean', 'u2', 'matparams', mfvm2)
mfvm1.export_to_vtk('displacement_with_von_mises1.vtk',
mfvm1, VM1, 'Von Mises Stresses',
mfu1, U1, 'Displacements',
md.mesh_fem_of_variable("lambda1"), md.variable("lambda1"), 'lambda1',
md.mesh_fem_of_variable("mult"), md.variable("mult"), 'mult')
mfvm2.export_to_vtk('displacement_with_von_mises2.vtk', mfvm2, VM2, 'Von Mises Stresses', mfu2, U2, 'Displacements')
print('You can view solutions with for instance:\nmayavi2 -d displacement_with_von_mises1.vtk -f WarpVector -m Surface -d displacement_with_von_mises2.vtk -f WarpVector -m Surface')
