Hi SfePy community,
I am attempting model stress and strain rate variations over a region (Omega) that contains subregions (Omega1 and Omega2) with different viscosities.
I have a general understanding of this mathematically, but i am still unsure the exact equations used by SfePy to calculate this. So far my model is largely based on the NaviesStokes example provided in SfePy with the circle_in_square mesh (example here: https://pastebin.com/VZrEKVFZ). Ultimately I would like to be able to model linear and power law viscous flow, is this possible with the existing equations in SfePy?
More specific questions:
What are default units for pressure (p) the values range from ~ 50 to 50? How can I change or specify units if i need to?
How do you calculate strain rate? Can I use the GradTerm class to do this?: sr = problem.evaluate('ev_grad.i.Omega(u)', mode='el_avg')
When I try to calculate Cauchy stress i get error:
ValueError: wrong arguments shapes for "+1.00e+00 * ev_cauchy_stress.4.Omega(fluid.viscosity, u)" term!
I think this is because ev_cauchy_stress requires a 3x3 matrix material and I have input a single viscosity value? What is the best way to solve this?
Thank you!
Hi Ben,
On 03/02/2018 01:21 AM, Ben Melosh wrote:
Hi SfePy community,
I am attempting model stress and strain rate variations over a region (Omega) that contains subregions (Omega1 and Omega2) with different viscosities.
I have a general understanding of this mathematically, but i am still unsure the exact equations used by SfePy to calculate this. So far my model is largely based on the NaviesStokes example provided in SfePy with the circle_in_square mesh (example here: https://pastebin.com/VZrEKVFZ). Ultimately I would like to be able to model linear and power law viscous flow, is this possible with the existing equations in SfePy?
More specific questions:
What are default units for pressure (p) the values range from ~ 50 to 50? How can I change or specify units if i need to?
SfePy does not handle units at all, so the units of the results are determined by units of the input parameters (material constants, loads, etc.). For example, SI units on input > SI units on output. Note that all the input units need to be consistent. Some consistent unit sets: http://www.engtips.com/viewthread.cfm?qid=355753
How do you calculate strain rate? Can I use the GradTerm class to do this?: sr = problem.evaluate('ev_grad.i.Omega(u)', mode='el_avg')
I do not understand what you mean by the strain rate, because the solution of your problem does not depend on time. Maybe you just need the ev_cauchy_strain applied to the velocity u?
When I try to calculate Cauchy stress i get error:
ValueError: wrong arguments shapes for "+1.00e+00 * ev_cauchy_stress.4.Omega(fluid.viscosity, u)" term!
I think this is because ev_cauchy_stress requires a 3x3 matrix material and I have input a single viscosity value? What is the best way to solve this?
Try using stiffness_from_lame() (see examples/linear_elasticity/linear_elastic.py)  in linear isotropic elasticity, the stress is, given the Lame parameters lambda, mu:
sigma = lambda tr(e) I + 2 mu e
and in incompressible fluid:
sigma = p I + 2 mu e(velocity)
so I guess you could set lambda = 0, and mu to the viscosity to get the shear stress part using ev_cauchy_stress.
r.
SfePy does not handle units at all, so the units of the results are determined by units of the input parameters (material constants, loads, etc.). For >example, SI units on input > SI units on output. Note that all the input units need to be consistent. Some consistent unit sets: http://www.eng>tips.com/viewthread.cfm?qid=355753
I do not understand what you mean by the strain rate, because the solution of your problem does not depend on time. Maybe you just need the >ev_cauchy_strain applied to the velocity u?
Try using stiffness_from_lame() (see examples/linear_elasticity/linear_elastic.py)  in linear isotropic elasticity, the stress is, given the Lame parameters >lambda, mu:
sigma = lambda tr(e) I + 2 mu e
and in incompressible fluid:
sigma = p I + 2 mu e(velocity)
so I guess you could set lambda = 0, and mu to the viscosity to get the shear stress part using ev_cauchy_stress.
r.
Thanks Robert,
I guess my main question regarding strain rate is weather there are equations in place to solve for nonnewtonian flow?
Ben
On 03/05/2018 09:09 PM, Ben Melosh wrote:
SfePy does not handle units at all, so the units of the results are determined by units of the input parameters (material constants, loads, etc.). For >example, SI units on input > SI units on output. Note that all the input units need to be consistent. Some consistent unit sets: http://www.eng>tips.com/viewthread.cfm?qid=355753
I do not understand what you mean by the strain rate, because the solution of your problem does not depend on time. Maybe you just need the >ev_cauchy_strain applied to the velocity u?
Try using stiffness_from_lame() (see examples/linear_elasticity/linear_elastic.py)  in linear isotropic elasticity, the stress is, given the Lame parameters >lambda, mu:
sigma = lambda tr(e) I + 2 mu e
and in incompressible fluid:
sigma = p I + 2 mu e(velocity)
so I guess you could set lambda = 0, and mu to the viscosity to get the shear stress part using ev_cauchy_stress.
r.
Thanks Robert,
I guess my main question regarding strain rate is weather there are equations in place to solve for nonnewtonian flow?
I forgot about that one :)
It is not there directly, but you can update the nonlinear parameter depending on the solution as in the examples [1, 2]. The examples do not deal with flows, but by analogy, the viscosity could be made dependent on the velocity from a previous pseudotime step.
Does that help?
r.
[1] http://sfepy.org/docdevel/examples/diffusion/poisson_field_dependent_materi... [2] http://sfepy.org/docdevel/examples/linear_elasticity/material_nonlinearity....
Yes this helps and appears to work, but leads me to a new issue.
In both the case of newtonian and a nonnewtonian (ie. nonlinear parameter depending on velocity) models, the post_process() and get_pars() functions evaluate an equation over a the entire mesh ('Omega' : 'all'). How do i evaluate different equations over regions of the mesh (e.g. 'Omega1' : 'cells of group 1').
For example, here is something i have tried:
stress = problem.evaluate('ev_cauchy_stress.i.Omega1(viscosity1.D, u) + ev_cauchy_stress.i.Omega2(viscosity2.D, u)', mode='el_avg') out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None)
this leads to an error that operands could not be broadcast together with different shapes.
In the case of a nonnewtonian simulation this is especially relevant and important as I want to simulate flow in two different regions with different viscosities. The get_pars() shown in example [2] in your reply only deals with one, how do i apply this to two regions of a single mesh?
Thanks Robert
b
On 03/14/2018 11:55 PM, Ben Melosh wrote:
Yes this helps and appears to work, but leads me to a new issue.
In both the case of newtonian and a nonnewtonian (ie. nonlinear parameter depending on velocity) models, the post_process() and get_pars() functions evaluate an equation over a the entire mesh ('Omega' : 'all'). How do i evaluate different equations over regions of the mesh (e.g. 'Omega1' : 'cells of group 1').
For example, here is something i have tried:
stress = problem.evaluate('ev_cauchy_stress.i.Omega1(viscosity1.D, u) + ev_cauchy_stress.i.Omega2(viscosity2.D, u)', mode='el_avg') out['cauchy_stress'] = Struct(name='output_data', mode='cell', data=stress, dofs=None)
this leads to an error that operands could not be broadcast together with different shapes.
Check the post_process() function in example [1]  you need to extend the values to the whole domain, if you want to save the values to the results file.
[1] http://sfepy.org/docdevel/examples/multi_physics/piezo_elasticity.html
In the case of a nonnewtonian simulation this is especially relevant and important as I want to simulate flow in two different regions with different viscosities. The get_pars() shown in example [2] in your reply only deals with one, how do i apply this to two regions of a single mesh?
It is true, that you get the coefficient (mu in the example [2]) in the whole domain, but in your case the array will have different values in Omega1 and Omega2. So will not be the viscosity also different automatically?
r.
participants (2)

Ben Melosh

Robert Cimrman