Fwd: SfePy  Electrothermomechanics problem
Hi Manuel,
I am forwarding your message to the mailing list, because I could not reach you using your address (Connection timed out error). I will answer to the list as well.
r.
 Forwarded Message  Subject: SfePy  Electrothermomechanics problem Date: Mon, 28 Aug 2017 12:31:33 +0200 From: Manuel Feuchter <manuel.feuchter@twtgmbh.de> To: cimrman3@ntc.zcu.cz
Hey Robert,
first of all: SfePy is an amazing project! I am getting to like SfePy more and more :)
I tried to start a thread at https://mail.python.org/mm3/archives/list/sfepy@python.org But unfortunately am I not able to. I recieve a "Server Error (500)" . So I try it via Email :)
Concerning SfePy: I get stuck with a few problems/questions: Background: I would like to set up a subsequent coupled system (so not fully coupled): like to apply the electric field as an time (time step) dependent function.)
 electric field evolves an electric current (at the end of the day, I would
 electric current results in liberated heat by joule heating
 the liberated heat drives the thermal expansion und results in thermal stresses
1:
Concerning example: http://sfepy.org/docdevel/examples/multi_physics/thermal_electric.html (multi_physics/thermal_electric.py) Equation: dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0 Your comment states: """ First solve the stationary electric conduction problem. Then use its results to solve the evolutionary heat conduction problem. """ Alright. However, just to understand it comprehensively? When you say you solve the "stationary electric conduction problem". Does this imply actually two equations? $(1) \textbf{E} = \nabla \phi$\\ $(2) \textbf{j} = \frac{1}{\varrho}\cdot\textbf{E}$ Is this right?
# 2:
Assuming I got it right, I would like to solve the electric conduction problem in terms of a quasistatic circumstance (external electric potential and internal resistivity will change). This means I would like to solve the electric conduction problem for every time step that I run the simulation and hence use the result (of $\textbf{j}$) for every next time step. I tried so, however I must have gotten sth wrong. The apparent electric field is not used?! Am I right? And how am I able to achieve so? (please find my code below)
# 3:
For the heat conduction equation: ... = dw_electric_source.2.Omega( m.electric_conductivity, s, phi ) Here, for the source term $\textbf{p} = \varrho \cdot \textbf{j}^2$ is solved? And how does the specific heat gets involved in your equation? Is the source term divided by the thermal conductivity [W/mK]? Does specific_heat to be named exactly like this to be involved in the 'dw_volume_dot.2.Omega( s, dT/dt )' term?
# 4:
How am I able to apply a time (time step) dependent function for the electric potential?!
Looking forward for your reply! Cheers and kind regards from Germany Manuel
def stress_strain_ext(out, pb, state, extend=False): """ Calculate and output strain and stress for given displacements. """ from sfepy.base.base import Struct from sfepy.mechanics import tensors
ev = pb.evaluate
strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.2.Omega(m.mech_solid, u)', mode='el_avg')
vms = tensors.get_von_mises_stress(stress.squeeze())
vms.shape = (vms.shape[0], 1, 1, 1)
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain, dofs=None)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress, dofs=None)
out['von_mises_stress'] = Struct(name='output_data', mode='cell',
data=vms, dofs=None)
return out
#! Presettings #!  t0 = 0.0 t1 = 100.0 n_step = 11 # material parameters and constants lam = 10.0 mu = 5.0 thermal_expandability = 1.25e5 T0 = 20.0 # Background temperature. specific_heat = 1.2 thermal_kappa_xyz = 2.0 electric_sigma_xyz = 1.5
#! Options #!  options = { 'ts': 'ts', 'nls': 'newton', 'ls': 'ls', 'output_dir' : output_dir, 'post_process_hook' : 'stress_strain_ext', }
#! Regions #!  regions = { 'Omega': 'all', 'Omega_Cell_Surface': ('vertices of surface', 'facet'), 'Mech_Mount': ('vertices in (x < 0.05)', 'facet'), 'Mech_Front': ('vertices in (x > 0.95)', 'facet'), 'Therm_inlet': ('vertices in (x < 90.05) & (x > 85.95) & (y > 9.95)', 'facet'), 'Therm_outlet': ('vertices in (x < 0.05)', 'facet'), 'Elec_Phi_0' : ('vertices in (x < 37.05) & (x > 32.95) & (y > 9.95)', 'facet'), 'Elec_Phi_1' : ('vertices in (x < 90.05) & (x > 85.95) & (y < 0.05)', 'facet'), }
#! Materials #!  eye_sym = np.array([[1], [1], [1], [0], [0], [0]], dtype=np.float64) materials = { 'm': ({ 'thermal_conductivity': thermal_kappa_xyz, 'electric_conductivity': electric_sigma_xyz, 'mech_solid' : stiffness_from_lame(3, lam=lam, mu=mu), 'alpha_therm' : (3.0 * lam + 2.0 * mu) * thermal_expandability * eye_sym },), 'mech_load': ({'.val': [0.0, 0.0, 0.0]},), }
#! Fields #!  fields = { 'displacement': ('real', 'vector', 'Omega', 1), 'temperature': ('real', 1, 'Omega', 1), 'potential' : ('real', 1, 'Omega', 1), }
#! Variables #!  variables = { 'u': ('unknown field', 'displacement', 2), 'v': ('test field', 'displacement', 'u'), 'T': ('unknown field', 'temperature', 0, 1), # 1 means, set history "on" 's': ('test field', 'temperature', 'T'), 'phi': ('unknown field', 'potential', 1, 1), 'psi': ('test field', 'potential', 'phi'), }
#! Initial Conditions #!  ics = { 'ic' : ('Omega', {'u.all' : 0.0, 'T.0': 0.0, 'phi.all': 0.0, }), }
#! Boundary Conditions #!  ebcs = { 'u0': ('Mech_Mount', {'u.all': 0.0}), 't0': ('Mech_Mount', {'T.0': 0.0}), 't1': ('Mech_Front', {'T.0': 0.0}), 'p0' : ('Elec_Phi_0', {'phi.0' : 0.0}), 'p1' : ('Elec_Phi_1', {'phi.0' : 100.0}), }
#! Equations #!  equations = { 'potential' : """dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0""", 'heat_conduct' : """dw_volume_dot.2.Omega( s, dT/dt ) + dw_laplace.2.Omega( m.thermal_conductivity, s, T ) = dw_electric_source.2.Omega( m.electric_conductivity, s, phi )""", 'thermo_mech': """dw_lin_elastic.2.Omega(m.mech_solid, v, u)  dw_biot.2.Omega(m.alpha_therm, v, T) = 0""", }
#! Solvers #!  solvers = { 'ls' : ('ls.scipy_direct', {}), 'newton' : ('nls.newton', { 'i_max' : 1, 'eps_a' : 1e10, 'eps_r' : 1.0, 'problem': 'nonlinear', }), 'ts': ('ts.simple', { 't0': t0, 't1': t1, 'dt': None, 'n_step': n_step, }), }
Schöne Grüße kind regards
i.A. Manuel Feuchter Competence Unit Engineering
Dr.Ing. Manuel Feuchter Head of Multibody & Structure Analysis
TWT GmbH  Science & Innovation Ernsthaldenstraße 17 D70565 Stuttgart EMail: manuel.feuchter@twtgmbh.de Mobil: +49 162 / 212 48 16
www.twtgmbh.de
Geschäftsführung: Dr. Dimitrios Vartziotis, Joachim Laicher (Stv.), Frank Beutenmüller (Stv.) Registergericht: Amtsgericht Stuttgart, HRB Nr. 212778 Umsatzsteuer: IDNr.: DE147841145
On 08/29/2017 02:31 PM, Robert Cimrman wrote:
Hi Manuel,
I am forwarding your message to the mailing list, because I could not reach you using your address (Connection timed out error). I will answer to the list as well.
r.
 Forwarded Message  Subject: SfePy  Electrothermomechanics problem Date: Mon, 28 Aug 2017 12:31:33 +0200 From: Manuel Feuchter <manuel.feuchter@twtgmbh.de> To: cimrman3@ntc.zcu.cz
Hey Robert,
first of all: SfePy is an amazing project! I am getting to like SfePy more and more :)
Thanks!
I tried to start a thread at https://mail.python.org/mm3/archives/list/sfepy@python.org But unfortunately am I not able to. I recieve a "Server Error (500)" . So I try it via Email :)
Concerning SfePy: I get stuck with a few problems/questions: Background: I would like to set up a subsequent coupled system (so not fully coupled): like to apply the electric field as an time (time step) dependent function.)
 electric field evolves an electric current (at the end of the day, I would
 electric current results in liberated heat by joule heating
 the liberated heat drives the thermal expansion und results in thermal stresses
1:
Concerning example: http://sfepy.org/docdevel/examples/multi_physics/thermal_electric.html (multi_physics/thermal_electric.py) Equation: dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0 Your comment states: """ First solve the stationary electric conduction problem. Then use its results to solve the evolutionary heat conduction problem. """ Alright. However, just to understand it comprehensively? When you say you solve the "stationary electric conduction problem". Does this imply actually two equations? $(1) \textbf{E} = \nabla \phi$\\ $(2) \textbf{j} = \frac{1}{\varrho}\cdot\textbf{E}$ Is this right?
It actually solves the electrostatic problem $\nabla \cdot \textbf{j} = 0$, cf. [1]  the Laplace's equation. Note that SfePy cannot be used to solve problems involving timevarying magnetic fields (full Maxwell's equations), because the vector finite elements needed for that are not implemented. Only irrotational fields (rot E = 0) can be solved, as in that example.
[1] https://en.wikipedia.org/wiki/Electrostatics#Poisson_and_Laplace_equations
# 2:
Assuming I got it right, I would like to solve the electric conduction problem in terms of a quasistatic circumstance (external electric potential and internal resistivity will change). This means I would like to solve the electric conduction problem for every time step that I run the simulation and hence use the result (of $\textbf{j}$) for every next time step. I tried so, however I must have gotten sth wrong. The apparent electric field is not used?! Am I right? And how am I able to achieve so? (please find my code below)
You need to do the timeloop and the data setting yourself, similarly to what is done in the example  check its main() function. Specifically, the results from the electrostatic computation (line 129) are used to initialize the source term of the heat conduction (line 135), and then the time solver iterates over time steps (lines 138139).
Your case is not that simple, as you need to have all the equations available in every time step, so you cannot use just a single Problem instance for that. You can try creating three Problem instances, for the electrostatics, the heat conduction, and the elasticity, and pass data among them as needed in a time loop: after problem.setup_default_output() line, use problem.copy() to create the other two instances and problem.set_equations() to set their equations.
But it might be better to abandon the problem description file approach, and build your subproblems using the sfepy classes as in the interactive examples
 I am not sure. Could you send me your complete equations (preferably already in the weak form) discretized in time, so that the data dependencies are clear?
# 3:
For the heat conduction equation: ... = dw_electric_source.2.Omega( m.electric_conductivity, s, phi ) Here, for the source term $\textbf{p} = \varrho \cdot \textbf{j}^2$ is solved?
Yes, and it is assumed that the electrostatic potential $\phi$ is known. See [2] for the term definition.
[2] http://sfepy.org/docdevel/src/sfepy/terms/terms_electric.html#sfepy.terms.t...
And how does the specific heat gets involved in your equation? Is the source term divided by the thermal conductivity [W/mK]? Does specific_heat to be named exactly like this to be involved in the 'dw_volume_dot.2.Omega( s, dT/dt )' term?
No, it is just substituted into the string defining the equation. It is equivalent to writing '1.2 * dw_volume_dot.2.Omega( s, dT/dt )'. You can also use 'dw_volume_dot.2.Omega(m.specific_heat, s, dT/dt)', if you define 'specific_heat' in the material 'm'.
# 4:
How am I able to apply a time (time step) dependent function for the electric potential?!
For the electrostatic problem material parameters? Use a function for setting those. Unfortunately, the docs are focused on the problem description file approach, while you need to do this in a script (interactively), to be able to pass data among several problems. The only example that sets a material parameter by a function is in create_local_problem() function of [3]. Let me know if you need help with that.
[3] http://sfepy.org/docdevel/examples/diffusion/poisson_parallel_interactive.h...
Cheers! r.
<snip>
Hey Robert,
alright. Thank you very much in advance!
Schöne Grüße & kind regards from Germany
i.A. Manuel Feuchter Competence Unit Engineering
Dr.Ing. Manuel Feuchter Head of Multibody & Structure Analysis
TWT GmbH  Science & Innovation Ernsthaldenstraße 17 D70565 Stuttgart EMail: manuel.feuchter@twtgmbh.de Mobil: +49 162 / 212 48 16
www.twtgmbh.de
Geschäftsführung: Dr. Dimitrios Vartziotis, Joachim Laicher (Stv.), Frank Beutenmüller (Stv.) Registergericht: Amtsgericht Stuttgart, HRB Nr. 212778 Umsatzsteuer: IDNr.: DE147841145
Von: Robert Cimrman <cimrman3@ntc.zcu.cz> An: sfepy@python.org Datum: 29.08.2017 14:31 Betreff: [sfepy] Fwd: SfePy  Electrothermomechanics problem
Hi Manuel,
I am forwarding your message to the mailing list, because I could not reach you using your address (Connection timed out error). I will answer to the list as well.
r.
 Forwarded Message  Subject: SfePy  Electrothermomechanics problem Date: Mon, 28 Aug 2017 12:31:33 +0200 From: Manuel Feuchter <manuel.feuchter@twtgmbh.de> To: cimrman3@ntc.zcu.cz
Hey Robert,
first of all: SfePy is an amazing project! I am getting to like SfePy more and more :)
I tried to start a thread at https://mail.python.org/mm3/archives/list/sfepy@python.org But unfortunately am I not able to. I recieve a "Server Error (500)" . So I try it via Email :)
Concerning SfePy: I get stuck with a few problems/questions: Background: I would like to set up a subsequent coupled system (so not fully coupled): like to apply the electric field as an time (time step) dependent
 electric field evolves an electric current (at the end of the day, I would
function.) 2. electric current results in liberated heat by joule heating 3. the liberated heat drives the thermal expansion und results in thermal stresses
# 1:
Concerning example: http://sfepy.org/docdevel/examples/multi_physics/thermal_electric.html (multi_physics/thermal_electric.py) Equation: dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0 Your comment states: """ First solve the stationary electric conduction problem. Then use its results to solve the evolutionary heat conduction problem. """ Alright. However, just to understand it comprehensively? When you say you solve the "stationary electric conduction problem". Does this imply actually two
equations? $(1) \textbf{E} = \nabla \phi$\\ $(2) \textbf{j} = \frac{1}{\varrho}\cdot\textbf{E}$ Is this right?
# 2:
Assuming I got it right, I would like to solve the electric conduction problem in terms of a quasistatic circumstance (external electric potential and internal resistivity will change). This means I would like to solve the electric conduction problem for every time step that I run the simulation and hence use the result (of $\textbf{j}$) for every next time step. I tried so, however I must have gotten sth wrong. The apparent electric field is not used?! Am I right? And how am I able to achieve so? (please find my code below)
# 3:
For the heat conduction equation: ... = dw_electric_source.2.Omega( m.electric_conductivity, s, phi ) Here, for the source term $\textbf{p} = \varrho \cdot \textbf{j}^2$ is solved? And how does the specific heat gets involved in your equation? Is the source term divided by the thermal conductivity [W/mK]? Does specific_heat to be named exactly like this to be involved in the 'dw_volume_dot.2.Omega( s, dT/dt )' term?
# 4:
How am I able to apply a time (time step) dependent function for the electric potential?!
Looking forward for your reply! Cheers and kind regards from Germany Manuel
def stress_strain_ext(out, pb, state, extend=False): """ Calculate and output strain and stress for given displacements. """ from sfepy.base.base import Struct from sfepy.mechanics import tensors
ev = pb.evaluate
strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = ev('ev_cauchy_stress.2.Omega(m.mech_solid, u)',
mode='el_avg')
vms = tensors.get_von_mises_stress(stress.squeeze())
vms.shape = (vms.shape[0], 1, 1, 1)
out['cauchy_strain'] = Struct(name='output_data', mode='cell',
data=strain, dofs=None)
out['cauchy_stress'] = Struct(name='output_data', mode='cell',
data=stress, dofs=None)
out['von_mises_stress'] = Struct(name='output_data', mode='cell',
data=vms, dofs=None)
return out
#! Presettings #!  t0 = 0.0 t1 = 100.0 n_step = 11 # material parameters and constants lam = 10.0 mu = 5.0 thermal_expandability = 1.25e5 T0 = 20.0 # Background temperature. specific_heat = 1.2 thermal_kappa_xyz = 2.0 electric_sigma_xyz = 1.5
#! Options #!  options = { 'ts': 'ts', 'nls': 'newton', 'ls': 'ls', 'output_dir' : output_dir, 'post_process_hook' : 'stress_strain_ext', }
#! Regions #!  regions = { 'Omega': 'all', 'Omega_Cell_Surface': ('vertices of surface', 'facet'), 'Mech_Mount': ('vertices in (x < 0.05)', 'facet'), 'Mech_Front': ('vertices in (x > 0.95)', 'facet'), 'Therm_inlet': ('vertices in (x < 90.05) & (x > 85.95) & (y > 9.95)',
'facet'), 'Therm_outlet': ('vertices in (x < 0.05)', 'facet'), 'Elec_Phi_0' : ('vertices in (x < 37.05) & (x > 32.95) & (y > 9.95)',
'facet'), 'Elec_Phi_1' : ('vertices in (x < 90.05) & (x > 85.95) & (y < 0.05)',
'facet'), }
#! Materials #!  eye_sym = np.array([[1], [1], [1], [0], [0], [0]], dtype=np.float64) materials = { 'm': ({ 'thermal_conductivity': thermal_kappa_xyz, 'electric_conductivity': electric_sigma_xyz, 'mech_solid' : stiffness_from_lame(3, lam=lam, mu=mu), 'alpha_therm' : (3.0 * lam + 2.0 * mu) * thermal_expandability * eye_sym },), 'mech_load': ({'.val': [0.0, 0.0, 0.0]},), }
#! Fields #!  fields = { 'displacement': ('real', 'vector', 'Omega', 1), 'temperature': ('real', 1, 'Omega', 1), 'potential' : ('real', 1, 'Omega', 1), }
#! Variables #!  variables = { 'u': ('unknown field', 'displacement', 2), 'v': ('test field', 'displacement', 'u'), 'T': ('unknown field', 'temperature', 0, 1), # 1 means, set history "on" 's': ('test field', 'temperature', 'T'), 'phi': ('unknown field', 'potential', 1, 1), 'psi': ('test field', 'potential', 'phi'), }
#! Initial Conditions #!  ics = { 'ic' : ('Omega', {'u.all' : 0.0, 'T.0': 0.0, 'phi.all': 0.0, }), }
#! Boundary Conditions #!  ebcs = { 'u0': ('Mech_Mount', {'u.all': 0.0}), 't0': ('Mech_Mount', {'T.0': 0.0}), 't1': ('Mech_Front', {'T.0': 0.0}), 'p0' : ('Elec_Phi_0', {'phi.0' : 0.0}), 'p1' : ('Elec_Phi_1', {'phi.0' : 100.0}), }
#! Equations #!  equations = { 'potential' : """dw_laplace.2.Omega( m.electric_conductivity, psi, phi ) = 0""", 'heat_conduct' : """dw_volume_dot.2.Omega( s, dT/dt ) + dw_laplace.2.Omega( m.thermal_conductivity, s, T ) = dw_electric_source.2.Omega( m.electric_conductivity, s, phi )""", 'thermo_mech': """dw_lin_elastic.2.Omega(m.mech_solid, v, u)  dw_biot.2.Omega(m.alpha_therm, v, T) = 0""", }
#! Solvers #!  solvers = { 'ls' : ('ls.scipy_direct', {}), 'newton' : ('nls.newton', { 'i_max' : 1, 'eps_a' : 1e10, 'eps_r' : 1.0, 'problem': 'nonlinear', }), 'ts': ('ts.simple', { 't0': t0, 't1': t1, 'dt': None, 'n_step': n_step, }), }
Schöne Grüße kind regards
i.A. Manuel Feuchter Competence Unit Engineering
Dr.Ing. Manuel Feuchter Head of Multibody & Structure Analysis
TWT GmbH  Science & Innovation Ernsthaldenstraße 17 D70565 Stuttgart EMail: manuel.feuchter@twtgmbh.de Mobil: +49 162 / 212 48 16
www.twtgmbh.de
Geschäftsführung: Dr. Dimitrios Vartziotis, Joachim Laicher (Stv.), Frank Beutenmüller (Stv.) Registergericht: Amtsgericht Stuttgart, HRB Nr. 212778 Umsatzsteuer: IDNr.: DE147841145
SfePy mailing list sfepy@python.org https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/
participants (2)

Manuel Feuchter

Robert Cimrman