Mailman 3 May 2023 - SfePy

Non linear material
by J. Sanz 12 Jul '23

12 Jul '23

Hello, I have a problem with a simulation including non lineal material.The problem involves a cantilever beam. The free end is moved 20 mm in X direction. I have compared the results with Ansys. The stresses are different, in particular at the pastification zone. Yield point of the material Relastic = 130 MPa, Epsilon elastic =0.00173 The material is bilinear, it has two slopes, young 1 and young2. The slope changes in the yield point. The no linear material function is based on calculating the strains in every cell. It is assigned the second young modulus for those strains higher than epsilon elastic. If lower, first young modulus is set. You can find attached the files and comparison. The picture with the comparison has the same scale, so the colours represent the same stress range. My doubts are the next: 1-I am not sure if the material definition is right. Can you confirm this for me? 2-Order of the stress strain evaluation, I set it to 3 because it provides the closest results. 3-The output file .vtk contains only nodal solution. How to get the elemental solution as well? Thanks, Jesus Sanz

2 5

Fwd: Re: Problem with Modal Analysis
by Vladimir Lukes 30 May '23

30 May '23

-------- Forwarded Message -------- Subject: Re: [sfepy] Problem with Modal Analysis Date: Wed, 24 May 2023 15:36:40 +0200 From: Vladimir Lukes <vlukes(a)kme.zcu.cz> To: J. Sanz <jsanz451(a)gmail.com> Hi Jesus, try the new PRIMME eigenvalue solver which is now integrated in SfePy. Have a look at the following pull requests: https://github.com/sfepy/sfepy/pull/995 https://github.com/sfepy/sfepy/pull/982 Regards Vladimir On 19. 05. 23 9:29, J. Sanz wrote: > Hello, > I have calculated the modal frequencies of a part with Sfepy. The > results have been compared to Ansys. Unfortunately, the results > are different. > > I have repeated the calculation by changing the next parameters trying > to find stable results. > -Solver: eig.scipy/eigsh > > -Solver:eig.scipy_lobpcg - > -i_max(iterations) > -Mesh size. I have tested different mesh from 12k cells to 22k cells. > -Order > > I am not able to get results in a reasonable time by using > eig.scipy/eigsh . (as reference for 13k cell, time is around 10 min) > The calculation runs with the solver eig.scipy_lobpcg. In all the > tests, I am not able to reach good tolerance. So, the results are quite > dependent on the number of elements, i_max parameter, and order. On the > other hand, if I try to make it more accurate, the time processing is > too long. > > Could you give any hint or recommendation? > > The simulation condition is free, without boundary conditions.The first > modes from Ansys are 4746 Hz, 5118 Hz, 5426 Hz, 6906 Hz. The six first > modes corresponding to rigid motion have been discarded. > The modes with Sfepy are 8710 Hz, 12954 Hz, 13138 Hz , 15213 Hz > Ansys use tetrahedrons with quadratic formulation (midside nodes) and > fine mesh. > > You can find attached code and commands to run it. > To run the model write in the terminal: > python modal_analysis.py Dummy_Modal_4_1.mesh > --solver="eig.scipy_lobpcg,i_max:1000,largest:False" --young 74854 > --density 2.7e-9 --poisson 0.3 --order 1 > Units: mm, N, MPa, Tn/mm3 > Thanks, > Jesus S. > > > _______________________________________________ > SfePy mailing list -- sfepy(a)python.org > To unsubscribe send an email to sfepy-leave(a)python.org > https://mail.python.org/mailman3/lists/sfepy.python.org/ > Member address: vlukes(a)kme.zcu.cz

2 2

Problem with Modal Analysis
by J. Sanz 23 May '23

23 May '23

Hello, I have calculated the modal frequencies of a part with Sfepy. The results have been compared to Ansys. Unfortunately, the results are different. I have repeated the calculation by changing the next parameters trying to find stable results. -Solver: eig.scipy/eigsh -Solver:eig.scipy_lobpcg - -i_max(iterations) -Mesh size. I have tested different mesh from 12k cells to 22k cells. -Order I am not able to get results in a reasonable time by using eig.scipy/eigsh . (as reference for 13k cell, time is around 10 min) The calculation runs with the solver eig.scipy_lobpcg. In all the tests, I am not able to reach good tolerance. So, the results are quite dependent on the number of elements, i_max parameter, and order. On the other hand, if I try to make it more accurate, the time processing is too long. Could you give any hint or recommendation? The simulation condition is free, without boundary conditions.The first modes from Ansys are 4746 Hz, 5118 Hz, 5426 Hz, 6906 Hz. The six first modes corresponding to rigid motion have been discarded. The modes with Sfepy are 8710 Hz, 12954 Hz, 13138 Hz , 15213 Hz Ansys use tetrahedrons with quadratic formulation (midside nodes) and fine mesh. You can find attached code and commands to run it. To run the model write in the terminal: python modal_analysis.py Dummy_Modal_4_1.mesh --solver="eig.scipy_lobpcg,i_max:1000,largest:False" --young 74854 --density 2.7e-9 --poisson 0.3 --order 1 Units: mm, N, MPa, Tn/mm3 Thanks, Jesus S.

2 3

Two problems in one process
by Yosef Meller 22 May '23

22 May '23

I'm trying to solve two problems consecutively in the same program. Each problem defines its own field and variables. Specifically, I do this for testing (same problem creation and operation code run consecutively with different parameters), just for the context. It fails when the second problem defines its field: ``` ValueError: order 0 already used! ``` which means that somewhere, there's a global mapping of numbers to variables which I need to clean. I don't know how to do it, and if there are other such global states that need to be cleaned. So my question is whether there's some "nuke everything" function, or a checklist of things to clean, or some other way to do it? My current workaround is to stick each run under a different process using the multiprocessing module. This works but I hope there's a less cumbersome way. Thanks, Yosef Meller.

2 3

Chaining boundary conditions
by Yosef Meller 22 May '23

22 May '23

Hi All, I have a boundary condition on a region that's defined by a function, and applies the same (time-dependent) deflection to all coordinates. Now I'm looking to precede this part of the movement by another part, which is more specific and changes by coordinate. This means that the previously simple part needs to provide not the same uniform deflection to all coordinates, but some accumulated deflection that adds to the last deflection in the first part of the movement. That is, previously I had d(x) = f(t), and now I want d(x) = d_0(x) + f(t - t_0), where 0 is the end of the 1st part. Chaining the two parts technically is simple, either manually or (probably) with the times keyword to EssentialBC. But the second part still needs to have d_0(x). I could, I guess, repeat the calculation for the last frame of the 1st part in each call to the 2nd-part function, but that seems excessive. What would be the best way to get that information in a BC function, or else to achieve that chaining by different means? Thanks, Yosef Meller.

2 3

Probing failed due to undefined integral
by 2232508＠tongji.edu.cn 13 May '23

13 May '23

Hi sfepy development team! I try to probe the Cauchy stress of my FEM model. But errors occurred when defining the evaluation of the problem. The full exception went like: Traceback (most recent call last): File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\base\base.py", line 556, in getitem ii = self.names.index(ii) ValueError: 'ii' is not in list During handling of the above exception, another exception occurred: Traceback (most recent call last): File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\integrals.py", line 53, in get obj = self[name] File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\base\base.py", line 563, in getitem raise IndexError(msg) IndexError: 'ii' is not in list During handling of the above exception, another exception occurred: Traceback (most recent call last): File "D:\PycharmWorkSpace\Box girder\FEM.py", line 129, in stress = pb.evaluate('ev_cauchy_stress.ii.girder(m.D, u)', mode='evaluate', File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\problem.py", line 1829, in evaluate aux = self.create_evaluable(expression, File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\problem.py", line 1771, in create_evaluable out = create_evaluable(expression, self.fields, materials, File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\evaluate.py", line 230, in create_evaluable equations = Equations.from_conf({'tmp' : expression}, File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\equations.py", line 65, in from_conf eq = Equation.from_desc(name, desc, variables, regions, File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\equations.py", line 743, in from_desc terms = Terms.from_desc(term_descs, regions, integrals) File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\terms\terms.py", line 218, in from_desc term = Term.from_desc(constructor, td, region, integrals=integrals) File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\terms\terms.py", line 354, in from_desc integral = integrals.get(desc.integral) File "D:\ProgramFiles\Anaconda3\envs\FEM_39\lib\site-packages\sfepy\discrete\integrals.py", line 56, in get raise ValueError('integral %s is not defined!' % name) ValueError: integral ii is not defined! However, I have defined the integral ii in the code. The full code is listed below: from __future__ import absolute_import import numpy as np import pandas as pd from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Function, Equation, Equations, Problem) from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.terms import Term from sfepy.base.base import Struct from sfepy.postprocess.probes_vtk import Probe from six.moves import range hf = 1 wf = 2 hr = 5 wr = 1 hs = 1 ws = 6 theta = 0.25 * np.pi span_len = 2 p = -0.5 * 28.14 / (wf + (hr + hs) / np.tan(theta) + ws * 0.5) def get_bearing_point(coors, domain=None): coord_df = pd.DataFrame(coors) bottom_coord = coord_df[(coord_df[2] == coord_df[2].min()) & (coord_df[1] == coord_df[1].min())] # 梁端底面 min_x = bottom_coord[0].min() max_x = bottom_coord[0].max() x_1_4 = (max_x - min_x) / 4 + min_x x_3_4 = 3 * (max_x - min_x) / 4 + min_x bearing_point_1 = bottom_coord[abs(bottom_coord[0] - x_1_4) < 0.5] bearing_point_2 = bottom_coord[abs(bottom_coord[0] - x_3_4) < 0.5] flag = np.array( list(set( list(bearing_point_1.index) + list(bearing_point_2.index) )) ) return flag def linear_compress(ts, coor, pressure, mode=None, **kwargs): if mode == 'qp': val = np.tile(-1 * pressure, (coor.shape[0], 1, 1)) return {'val': val} pressure = Function('linear_compress', linear_compress, extra_args={'pressure': p}) functions = { 'get_bearing_point': (get_bearing_point,) } mesh = Mesh.from_file('box_girder.vtk') domain = FEDomain('domain', mesh) bounding_box = domain.get_mesh_bounding_box() min_y, max_y = bounding_box[:, 1] min_z, max_z = bounding_box[:, 2] eps_y = 1e-8 * (max_y - min_y) eps_z = 1e-8 * (max_z - min_z) bearing = domain.create_region('bearing', 'vertices by get_bearing_point') # 梁端固定支座 bearing = domain.create_region('bearing', 'vertices in y < %.10f' % (min_y + eps_y), 'facet') deck = domain.create_region('deck', 'vertices in z > %.10f' % (max_z - eps_z), 'facet') # 梁顶面 girder = domain.create_region('girder', 'all') field = Field.from_args('fu', np.float64, 'vector', girder, approx_order=2) u = FieldVariable('u', 'unknown', field) v = FieldVariable('v', 'test', field, primary_var_name='u') m = Material('m', D=stiffness_from_youngpoisson(3, 3.5e10, 0.2)) f = Material('f', val=[[0], [0], [-23520]]) load = Material('load', function=pressure) integral = Integral('ii', order=3) t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, girder, m=m, v=v, u=u) # 弹性能积分项 t2 = Term.new('dw_volume_lvf(f.val, v)', integral, girder, f=f, v=v) # 体积力积分项 t3 = Term.new('dw_surface_ltr(load.val, v)', integral, deck, load=load, v=v) # 面力积分项 eq = Equation('balance', t1 - t2 - t3) eqs = Equations([eq]) fix_u = EssentialBC('fix_u', bearing, {'u.all': 0.0}) ls = ScipyDirect({}) nls_status = IndexedStruct() nls = Newton({'i_max': 1, 'eps_a': 1e-10}, lin_solver=ls, status=nls_status) pb = Problem('elasticity', equations=eqs) pb.save_regions_as_groups('regions') pb.set_bcs(ebcs=Conditions([fix_u])) pb.set_solver(nls) status = IndexedStruct() vec = pb.solve(status=status) print('Nonlinear solver status:\n', nls_status) print('Stationary solver status:\n', status) pb.save_state('linear_elasticity.vtk', vec) out = {} stress = pb.evaluate('ev_cauchy_stress.ii.girder(m.D, u)', mode='evaluate') out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None) Besides the above question, I wonder if there is any tutorial that shows how to write postprocess code interactively. I would greatly appreciate a full example from problem defining to postprocessing with coding style shown in this page Tutorial — SfePy version: 2023.1+git.fdee68c7 documentation. By the way, the above problem was also subjected as issue#985 on GitHub. Evaluation failed · Issue #985 · sfepy/sfepy (github.com) I would greatly appreciate your support of any kind.

2 1

Re: How to get the value with specific coordinate.
by Robert Cimrman 03 May '23

03 May '23

On 5/2/23 15:02, 潘梁柱 wrote: > Thanks a lot, I have some question about post-process and export data file. > I've noticed that I can get "Mises_stress" by "get_von_mises_stress", but I don't know how to get the von_mises_strain. > Since .vtk file stored a tensor, how can I get the specific component of the cauchy strain. Sfepy-view only shows the picture of the magnitude of cauchy strain. Have a look at the beginning of post_process() function in [1]. The von Mises stress (= a scalar in each point) is obtained from the Cauchy stress tensor. > > > -----原始邮件----- > > 发件人: "Robert Cimrman" <cimrman3(a)ntc.zcu.cz> > > 发送时间: 2023-04-25 16:01:08 (星期二) > > 收件人: sfepy(a)python.org > > 抄送: > > 主题: [sfepy] Re: How to get the value with specific coordinate. > > > > Hello 潘梁柱, > > > > On 4/22/23 05:04, 潘梁柱 wrote: > > > I have got the vtk file, which contains cauchy-strain, after runing 'sfepy-run sf.py'in command line. But I want to get the strain at the point (0, 0),(0.5, 0.5). And it will be better if i can export a data file that contains strain with its corresponding coordinates. > > > > Try something like > > > > vu = pb.get_variables()['u'] > > vu.evaluate_at(nm.array([[0, 0], [0.5, 0.5]]), mode='cauchy_strain') > > > > in your postproc(). > > > > I still don't know how to get the strain value of the point [0.5, 0.5], after trying the above code and runing Sfepy-view command, I didn't see the strain value that I need appear in the Pyvista. Is the Sfepy-veiw step wrong? I just don't know where to see the "vu" defined in the code, is there any method I can save the value of "vu" in a text file? sfepy-view is just a simple results viewer. You need to put the above code into a post-process hook function, see again [1], where the function is called post_process(). This function is called after the primary variables (= those in 'variables') are computed, and it is meant for user-defined post-processing. In your case that would be computing the von Mises stress, evaluating data in given points, saving those into a file, etc. This is not done automatically - a user must write such a function explicitly. r. [1] http://sfepy.org/doc-devel/examples/linear_elasticity-nodal_lcbcs.html

1 0

Material units in Modal Simulation
by J. Sanz 02 May '23

02 May '23

Hello, I have a problem when running a modal simulation. There is not error, but the number of eingenvalues found is zero. May be the problem is realed to the units. The mesh is in mm, and young modulus is MPa (N,mm), poisson is adimensional. My question is the next, which units shoud have rho, Tons/mm3 or Kg/mm^3? Thanks, Jesus Sanz

2 3