Hi Kathrin,
On 08/25/2017 04:42 PM, Kathrin Sobe wrote:
Hi Robert,
thanks a lot for your help.
2017-08-25 14:37 GMT+02:00 Robert Cimrman <cimrman3@ntc.zcu.cz>:
Hi Kathrin,
On 08/25/2017 10:37 AM, Kathrin Sobe wrote:
Hi Robert,
I try to solve a linear elastic problem. I give a example in a python file that is a modified version of the example given in http://sfepy.org/doc-devel/examples/linear_elasticity/linear _elastic_tractions.html
The problem I want to solve is the following: I have a object, defined as a 3D mesh of tetrahedrons. In the test example I take the medium size 3D cube. The real 3D object is much bigger and has a more structured surface.
On the top of the object there are certain forces. These forces are defined as 3D vectors, one for each surface triangle. Additionally the object has certain material parameters given by density, poisson ratio and elastic modulus.
Working through the linear elastic examples I figured out two different types of solutions:
- modeling the forces at the surface as a material with a function (manual 4.3.9)
- modeling the forces at the surface as a boundary condition
So the first question is in general: Is this the right way to model my problem? At the moment I think solution type number one is the better approach?
Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary conditions = applied forces/tractions/pressure in the weak form of the elasticity PDE. Because you need a position-dependent value of the load, you need to use a material parameter given by a function.
ok, than I will concentrate on that approach. Additionally, I want to model the gravity. For that I added the rho parameter to the material defintion and defined a mass equation like this: # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr verformt: ist klar solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846), 'rho': 0.501}) load = Material('load', function=linear_tension) # material defined by values or function
integral = Integral('i', order=2) t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega,solid=solid, v=v, u=u) t2 = Term.new('dw_surface_ltr(load.val, v)', integral, surface, load=load, v=v) t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, v=v, u=u) eq1 = Equation('elasticity', t1 - t2) eq2 = Equation('mass', t3) eqs = Equations([eq1, eq2])
I am not sure if this is the right way to add the gravity to the problem.
Use the dw_volume_lvf term for volume forces just as you use dw_surface_ltr for surface forces. The material parameter here is a force per unit volume. The mass matrix would be useful in a dynamic simulation, not here.
Going through the first approach I have the following more detailed
problems:
- Region definition: I need a way to define the surface triangles of the 3D object. Because of the structure of the object, this cannot be done with coordinate ranges.
Is it going to be applied on the whole surface? Anyway, you can select the whole surface of a mesh with 'vertices of surface', for example:
surf = domain.create_region('surf', 'vertices of surface', 'facet')
Check the manual section 4.3.3 (on regions) for additional possibilities of defining regions and combining existing regions using various operations.
Or, if your mesh comes from some mesh generator/cad software, you can mark various vertex / cell groups there and use those for region definitions. This depends on the mesher/format you use - let me know if you need help with that.
The surface extraction is exactly what I need, cool. I found that out as well. Thanks for clarifying the indexing. I will try that.
Just as an information: Not all surface triangles have a force vector, only the triangles in the upper region of the geometry.
Then you may need to remove the extra faces from the surface region (or simply define the force zero there).
- Material definition with force vectors: the 3D force vectors are given
for each surface triangle (defined at the centroid of each triangle). In the material function examples I found only the definition of forces at quadrature points and only as scalar values.
You can use vectors there too, for example:
val = nm.tile([[-0.8], [0.1], [0.2]], (coor.shape[0], 1, 1))
To have a constant value per element, you need to know in which element or face each quadrature point is. This information is available in the material function (for example, kwargs['term'].region.faces contains the list of faces of the region, etc.), but the actual code depends on how your material function is represented - is it a table of (cell-value) pairs, or something else?
At the moment my force defintion is exacly as you expected: cell-value pairs. The cell is the surface triangle and the value is the force vector.
Is it that easy: Find the cell for each quadrature point and return the cell force vector as its load? Do I need to split the force vector of the cell between its quadrature points? Is there always exacly one cell for each quadrature point?
You need to provide a force per unit area vector (so yes, we may call it splitting) in each quadrature point. Each quadrature point is in exactly one face, yes. The number of quadrature point per face can be found as follows:
qps = t2.get_physical_qps() print qps.shape (88, 3, 3)
= 88 faces, 3 QPs per face, 3 coordinates (in 3D)
So this tells you, that the coor argument of linear_tension has 3 rows for each face (in the order of top.faces (= t2.region.faces).
r.