### How to implement an implicit constitutive relation in elasticity?

309

views

1

I have a governing constitutive relation where I can only write strain as a function of the stresses and not vice versa.

i.e.

I have

def epsilon(u):

return 0.5*(grad(u) + grad(u).T)

but defining sigma(u) like,

def sigma(u):

return 2.0*mu*epsilon(u) + lambda*tr(epsilon(u))*Identity(2)

is

Rather my sigma and epsilon are related by an implicit constitutive relation, where epsilon is a function of sigma and not vice versa.

i.e.

I have

def epsilon(u):

return 0.5*(grad(u) + grad(u).T)

but defining sigma(u) like,

def sigma(u):

return 2.0*mu*epsilon(u) + lambda*tr(epsilon(u))*Identity(2)

is

**not possible.**

Rather my sigma and epsilon are related by an implicit constitutive relation, where epsilon is a function of sigma and not vice versa.

Community: FEniCS Project

### 2 Answers

4

Most probably, you need to consider the stress as an independent variable. Then, you will have two equations (balance of linear momentum and the constitutive relation) instead of one as in the classical approach when you can substitute for the stress in the balance of linear momentum from the constitutive relation. See the possible weak formulations for implicitly constituted fluids: http://ncmm.karlin.mff.cuni.cz/db/attachments/single/417

1

Let me also comment on this. Basically, you can distinguish 3 cases:

1.) You know \( \varepsilon(\sigma) \), it's inverse is a well defined function but you just cannot invert it by hand. Then you can do what Adam suggested, introduce new unknown field and specify relation between \( \varepsilon \) and \( \sigma \) weakly in an additional equation. If you then use e.g. Newton solver, it is a newton iteration which inverts your relation numerically (but on a global matrix level!).

2.) You know \( \varepsilon(\sigma) \) but it's inverse is not well defined for some values. Then you cannot principally invert the relation. Anyway, you can do again what Adam suggested. But be careful of what is going on. If you use Newton solver, you are trying to find a root for a function that has more than 1 root, so the solution is not unique! and depends on the initial guess. So instead of simple Newton you might consider some step controlled newton as arc-length method.

3.) You cannot express \( \varepsilon(\sigma) \) nor \( \sigma(\varepsilon)\) as a function, you have fully implicit relation between them. You wont be able to solve the problem without additional knowledge, additional constitutive law and notion of some state variable. This is typical example for cyclic loading-unloading in plasticity.

1.) You know \( \varepsilon(\sigma) \), it's inverse is a well defined function but you just cannot invert it by hand. Then you can do what Adam suggested, introduce new unknown field and specify relation between \( \varepsilon \) and \( \sigma \) weakly in an additional equation. If you then use e.g. Newton solver, it is a newton iteration which inverts your relation numerically (but on a global matrix level!).

2.) You know \( \varepsilon(\sigma) \) but it's inverse is not well defined for some values. Then you cannot principally invert the relation. Anyway, you can do again what Adam suggested. But be careful of what is going on. If you use Newton solver, you are trying to find a root for a function that has more than 1 root, so the solution is not unique! and depends on the initial guess. So instead of simple Newton you might consider some step controlled newton as arc-length method.

3.) You cannot express \( \varepsilon(\sigma) \) nor \( \sigma(\varepsilon)\) as a function, you have fully implicit relation between them. You wont be able to solve the problem without additional knowledge, additional constitutive law and notion of some state variable. This is typical example for cyclic loading-unloading in plasticity.

Please login to add an answer/comment or follow this question.

The much more common approach is to linearize the balance of linear momentum as:

$\int\text{grad}\delta u\cdot C_k\cdot\text{grad}\Delta u_{k+1}\mathrm{d}x=\int\sigma_k\cdot\text{grad}\delta u\mathrm{d}x$∫grad

δu·C_{k}·gradΔu_{k+1}dx=∫σ_{k}·gradδudx(in the special case of linear kinematics and neither surface tractions nor volumetric forces are present)and calculate stiffness $C_k$

C_{k}and stresses $\sigma_k$σ_{k}in integration points in order to get your next displacement update $\Delta u_{k+1}$Δu_{k+1}for a global iteration.How you calculate the stresses is up to you. And if it's a local newton routine for inverting your strain-stress relationship then so be it.