Hello guys 
!
So i wanna calculate the components of a 2D stiffness tensor and i already have all the values of the stress (sigma11,sigma22,sigma12) and strains (epsilon11,epsilon22,epsilon12). How can i do that ?
Any ideas about how can we solve that ? 
Hi @CHAIkah,
It would help if you provide a little more detail about your problem. For example:
- Are you trying to solve a homogenization problem? (i.e. to find the effective stiffness tensor)
- What format do you have the stresses and strains in? For example, do you know the strains pointwise throughout some domain, e.g. as a FEniCS
Function, or did you compute the volume average so that the strains are stored in an array? - Are you studying plane stress or plane strain?
1 Like
Hey @conpierce8 !
Yes itâs a homogenization problem so Iâm trying to extract the effective stiffness tensor.
Itâs a 2D square divided into 10 phases each phase contains a different material under compression on the top and bottom boundaries.
I didnât specify in the code if itâs plane strain/stress
Itâs basically like what we have here Periodic Homogenization except that my model doesnât have any periodic BC.
Hi @CHAIkah,
Iâll be able to help you more effectively if you provide a code example. In the meantime, note the following:
- You need to study three different load cases (uniaxial strain in x, uniaxial strain in y, and shear strain) to get the full 2D stiffness tensor. (This is the approach taken in the periodic homogenization demo.)
- If you prescribe unit strain on the boundary, the 2D stiffness tensor components are simply equal to the volume-averaged stress. You can compute the volume averaged stress using
assembleas in the demo that you mentioned.
As an aside, you implicitly choose either plane strain or plane stress when you create the function sigma(v) that defines your constitutive model, as mentioned in the 2D linear elasticity demo.
1 Like
Hi @CHAIkah,
Since this isnât a working code, I canât provide a full solution. But, you should be able to use
sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
sigma_avg = assemble(sum(inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases))) / area
to compute the average stress sigma_avg. You should repeat this analysis three times with three different affine boundary conditions u = \bar{\varepsilon}^{(i)} x\quad \forall x \in \partial \Omega, where:
\bar{\varepsilon} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{12} \end{bmatrix} \\
\bar{\varepsilon}^{(1)} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix},\quad \bar{\varepsilon}^{(2)} = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}, \quad \bar{\varepsilon}^{(3)} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}
Then the average stress obtained for each case will give you the corresponding column of the 2D stiffness matrix.
you mean i need to study three different load cases (uniaxial strain in x , uniaxial strain in y , and shear strain and each time iâll get a column of the tensor ?
Thatâs correct. The âeffective stiffnessâ tensor C^{eff} is defined by the relation
\left< \sigma \right> = C^{eff} \left< \varepsilon \right>
where \left<\cdot\right> denotes the volume average. When you apply an affine boundary condition, the volume average of the strain \left< \varepsilon \right> is equal to the âboundaryâ strain \bar{\varepsilon}. Thus,
\left< \sigma \right> = C^{eff} \bar{\varepsilon}
If we take our 2D constitutive relation to be
C^{eff} = \begin{bmatrix} C^{eff}_{11} & C^{eff}_{12} & C^{eff}_{16} \\ C^{eff}_{12} & C^{eff}_{22} & C^{eff}_{26} \\ C^{eff}_{16} & C^{eff}_{26} & C^{eff}_{66} \end{bmatrix}
and apply a unit strain in the x-direction, for example, the volume average of the stress will be:
\left< \sigma \right> = \begin{bmatrix} C^{eff}_{11} & C^{eff}_{12} & C^{eff}_{16} \\ C^{eff}_{12} & C^{eff}_{22} & C^{eff}_{26} \\ C^{eff}_{16} & C^{eff}_{26} & C^{eff}_{66} \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} C^{eff}_{11} \\ C^{eff}_{12} \\ C^{eff}_{16} \end{bmatrix}
Repeating for \bar{\varepsilon} = [0\,1\,0]^T and \bar{\varepsilon} = [0\,0\,1]^T gives the second and third columns of the effective stiffness tensor, respectively.
i understand now. Thank you it works !
but @conpierce8 this is what i got and i donât understand why isnât symmetrical.
-1.572410583195420555e-03 2.216522418232252676e+03 1.108260422906548456e+03
4.414543445680684272e+03 -6.237213042085443552e-04 2.207271410980084511e+03
-6.237146402310145103e-04 1.941621083991628257e+04 9.708105108099545760e+03
@CHAIkah can you share your complete code? I suspect there may be an error with your boundary conditions, since the third column is a linear combination of the first two columns (within discretization error).
Hi @CHAIkah,
- Your affine Dirichlet boundary conditions should be applied on all boundaries.
- In order to directly extract the stiffness tensor components from the average stress, you need to apply a unit strain. Otherwise you will need to divide the stress by the strain to find the stiffness tensor components. To achieve a unit strain, you can apply the following affine BCs (corresponding to uniaxial strain in x, uniaxial strain in y, and pure shear strain, respectively):
u = [x\quad 0]^T\\
u = [0\quad y]^T\\
u = [y\quad x]^T
- Be careful with the ordering of your stresses. FEniCS gives you the four components of \left< \sigma \right> as a vector:
\left< \sigma \right> = [\sigma_{11} \quad \sigma_{12} \quad \sigma_{21} \quad \sigma_{22}]^T
Therefore, after you
AA = np.delete(AA, 2)
the ordering of the stresses in AA is:
AA[0] = sigma_xx
AA[1] = sigma_xy
AA[2] = sigma_yy
To achieve the familiar symmetric stiffness matrix, you must re-order AA (and BB and CC):
AA = np.delete(AA, 2)[[0, 2, 1]]
See the code below.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
fname = "lamin"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
#define materials
Em, num = 50e3 , 0.2
Er, nur = 210e3 , 0.3
material_parameters = [(Em, num), (Er, nur),(Em, num), (Er, nur),(Em, num),\
(Er, nur), (Em, num),(Er, nur), (Em, num),(Er, nur)]
nphases = len(material_parameters)
#define eps and sigma
def eps(v):
return sym(grad(v))
def sigma(v,i):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v))*Identity(2) + 2*mu*(eps(v))
# Define function space
W = VectorFunctionSpace(mesh,"Lagrange", 2)
w = Function(W)
v = TestFunction(W)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
Wint = sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)])
#uniaxial strain in y solution
bc_yy = Expression(("0", "x[1]"), degree=2)
bcs1 = [DirichletBC(W, bc_yy, facets, j) for j in [1,2,3,4]]
solve(Wint == 0, w, bcs1)
#compute eff stiffness tensor
area = 100
sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
A = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
AA =np.array(A)
AA = np.delete(AA, 2)[[0,2,1]]
#Shear strain solution
bc_xy = Expression(("x[1]", "x[0]"), degree=2)
bcs2 = [DirichletBC(W, bc_xy, facets, j) for j in [1,2,3,4]]
solve(Wint == 0, w, bcs2)
B = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
BB =np.array(B)
BB = np.delete(BB, 2)[[0,2,1]]
#uniaxial strain in x solution
bc_xx = Expression(("x[0]","0"), degree=2)
bcs3 = [DirichletBC(W, bc_xx, facets, j) for j in [1,2,3,4]]
solve(Wint == 0, w, bcs3)
C = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
CC =np.array(C)
CC = np.delete(CC, 2)[[0,2,1]]
#assemble tensor columns
dd = [[m, n, l]for m,n,l in zip(CC,AA,BB)]
np.savetxt("stiffnessmatrix.txt", dd)
On my machine, this produces the following output in stiffnessmatrix.txt:
1.526894178315140598e+05 3.272079349476299831e+04 -2.101344925113153399e-03
3.272079349479604207e+04 9.543267626480710169e+04 -4.449651732026041126e-03
-1.050665393418484094e-03 -2.224827810501039949e-03 6.996483562031936890e+04
which is symmetric to within discretization error.
2 Likes
After thinking about this problem some more, I have realized I should make one minor correction to my answer(s). The 2D strain vector is generally formulated as
\varepsilon = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \end{bmatrix} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \end{bmatrix}
Then the 2D elasticity relation is written as
\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{16} \\ C_{12} & C_{22} & C_{26} \\ C_{16} & C_{26} & C_{66} \end{bmatrix} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \end{bmatrix}
Thus, the DirichletBC for a unit shear strain should be
bc_xy = Expression(("x[1]/2", "x[0]/2"), degree=2)
rather than
bc_xy = Expression(("x[1]", "x[0]"), degree=2)
as I originally posted.
Yeah i rectified the mistake !
Thank you very much !
@conpierce8
if i wanna extract the displacement vector and strain and stress tensors in each subdomain, how can i do that from this code ?
this is what i tried to do for now:
for i in range(nphases) :
sig11 = assemble(sigma(w,i)[0,0]*dx(i))
sig22 = assemble(sigma(w,i)[1,1]*dx(i))
sig12 = assemble(sigma(w,i)[1,0]*dx(i))
eps11 = assemble(eps(w)[0,0]*dx(i))
eps22 = assemble(eps(w)[1,1]*dx(i))
eps12 = assemble(eps(w)[1,0]*dx(i))
ux = assemble(w[0]*dx(i))
uy = assemble(w[1]*dx(i))
and then print them using np.array
The above code should work, but note that you need to divide by the area of the subdomain in order to obtain the correct average stresses, strains, and displacements.
1 Like
Hi @conpierce8,
Thanks so much for your answering, which really helps me out on characterizing composites by Fenics. However, I am just a beginner of Fenics and I need to specify various stiffness matrices for each element. So would you mind providing a little more details on your geometry?
Hi @ShawnZ,
I donât have access to the original geometry since it appears the original poster has deleted some messages in this thread. Here is a working code and the corresponding meshes for a 2\times 2 square grid. The meshes were created with Gmsh but unfortunately I donât have the .geo file.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
# fname = "lamin"
fname = "grid"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
subdomains_xdmf = XDMFFile("subdomains.xdmf")
subdomains_xdmf.write(subdomains)
subdomains_xdmf.close()
facets_xdmf = XDMFFile("facets.xdmf")
facets_xdmf.write(facets)
facets_xdmf.close()
#define materials
Em, num = 50e3 , 0.2
Er, nur = 210e3 , 0.3
# material_parameters = [(Em, num), (Er, nur),(Em, num), (Er, nur),(Em, num),\
# (Er, nur), (Em, num),(Er, nur), (Em, num),(Er, nur)]
n_phases = 4
nu_random = np.random.normal(0.3, 0.04, 4)
material_parameters = [(Em, nu) for nu in nu_random]
nphases = len(material_parameters)
#define eps and sigma
def eps(v):
return sym(grad(v))
def sigma(v,i):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v))*Identity(2) + 2*mu*(eps(v))
# Define function space
W = VectorFunctionSpace(mesh,"Lagrange", 2)
w = Function(W)
v = TestFunction(W)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
Wint = sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)])
#uniaxial strain in y solution
bc_yy = Expression(("0", "x[1]"), degree=2)
w = project(bc_yy, W)
bcs1 = [DirichletBC(W, bc_yy, facets, j) for j in [1,2,3,4]]
solve(Wint == 0, w, bcs1)
#compute eff stiffness tensor
area = 100
sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
A = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
AA =np.array(A)
AA = np.delete(AA, 2)[[0,2,1]]
# #Shear strain solution
# bc_xy = Expression(("x[1]/2", "x[0]/2"), degree=2)
# w = project(bc_xy, W)
# bcs2 = [DirichletBC(W, bc_xy, facets, j) for j in [1,2,3,4]]
# solve(Wint == 0, w, bcs2)
# B = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
# BB =np.array(B)
# BB = np.delete(BB, 2)[[0,2,1]]
# #uniaxial strain in x solution
# bc_xx = Expression(("x[0]","0"), degree=2)
# w = project(bc_xx, W)
# bcs3 = [DirichletBC(W, bc_xx, facets, j) for j in [1,2,3,4]]
# solve(Wint == 0, w, bcs3)
# C = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
# CC =np.array(C)
# CC = np.delete(CC, 2)[[0,2,1]]
# #assemble tensor columns
# dd = [[m, n, l]for m,n,l in zip(CC,AA,BB)]
# np.savetxt("stiffnessmatrix.txt", dd)
grid.xml:
<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://www.fenicsproject.org">
<mesh celltype="triangle" dim="2">
<vertices size="13">
<vertex index="0" x="0.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="1" x="1.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="2" x="2.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="3" x="0.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="4" x="1.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="5" x="2.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="6" x="0.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="7" x="1.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="8" x="2.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="9" x="5.0000000000000000e-01" y="5.0000000000000000e-01" z="0.0000000000000000e+00"/>
<vertex index="10" x="1.5000000000000000e+00" y="5.0000000000000000e-01" z="0.0000000000000000e+00"/>
<vertex index="11" x="1.5000000000000000e+00" y="1.5000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="12" x="5.0000000000000000e-01" y="1.5000000000000000e+00" z="0.0000000000000000e+00"/>
</vertices>
<cells size="16">
<triangle index="0" v0="0" v1="1" v2="9"/>
<triangle index="1" v0="3" v1="0" v2="9"/>
<triangle index="2" v0="1" v1="4" v2="9"/>
<triangle index="3" v0="4" v1="3" v2="9"/>
<triangle index="4" v0="1" v1="2" v2="10"/>
<triangle index="5" v0="4" v1="1" v2="10"/>
<triangle index="6" v0="2" v1="5" v2="10"/>
<triangle index="7" v0="5" v1="4" v2="10"/>
<triangle index="8" v0="4" v1="5" v2="11"/>
<triangle index="9" v0="7" v1="4" v2="11"/>
<triangle index="10" v0="5" v1="8" v2="11"/>
<triangle index="11" v0="8" v1="7" v2="11"/>
<triangle index="12" v0="3" v1="4" v2="12"/>
<triangle index="13" v0="6" v1="3" v2="12"/>
<triangle index="14" v0="4" v1="7" v2="12"/>
<triangle index="15" v0="7" v1="6" v2="12"/>
</cells>
</mesh>
</dolfin>
grid_facet_region.xml:
<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_function type="uint" dim="1" size="28">
<entity index="0" value="1"/>
<entity index="1" value="4"/>
<entity index="2" value="0"/>
<entity index="3" value="1"/>
<entity index="4" value="0"/>
<entity index="5" value="0"/>
<entity index="6" value="0"/>
<entity index="7" value="2"/>
<entity index="8" value="0"/>
<entity index="9" value="0"/>
<entity index="10" value="4"/>
<entity index="11" value="0"/>
<entity index="12" value="0"/>
<entity index="13" value="0"/>
<entity index="14" value="0"/>
<entity index="15" value="0"/>
<entity index="16" value="0"/>
<entity index="17" value="0"/>
<entity index="18" value="0"/>
<entity index="19" value="2"/>
<entity index="20" value="0"/>
<entity index="21" value="0"/>
<entity index="22" value="3"/>
<entity index="23" value="0"/>
<entity index="24" value="3"/>
<entity index="25" value="0"/>
<entity index="26" value="0"/>
<entity index="27" value="0"/>
</mesh_function>
</dolfin>
grid_physical_region.xml:
<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_function type="uint" dim="2" size="16">
<entity index="0" value="1"/>
<entity index="1" value="1"/>
<entity index="2" value="1"/>
<entity index="3" value="1"/>
<entity index="4" value="2"/>
<entity index="5" value="2"/>
<entity index="6" value="2"/>
<entity index="7" value="2"/>
<entity index="8" value="3"/>
<entity index="9" value="3"/>
<entity index="10" value="3"/>
<entity index="11" value="3"/>
<entity index="12" value="4"/>
<entity index="13" value="4"/>
<entity index="14" value="4"/>
<entity index="15" value="4"/>
</mesh_function>
</dolfin>
1 Like
Thanks so much for your answer! @conpierce8
These really help me out. I took some time to generate the corresponding .geo successfully. However, I just went over the codes again and got two new questions:
(1) area = 100 which I suppose is pre-defined by scenario, but will it be related to the mesh size? E.g., for this solution, if taking a 4 by 4 mesh, shall we also change the area?
(2) In the file grid_physical_region.xml, four physical regions (1, 2, 3, 4) have been defined. But I am just confused about the line sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)]). range(nphases) will iterate as 0, 1, 2, 3 for sure, which sort of conflicts to the sequences of physical region. My feeling is to fix the dx but cannot make sure as a beginner of FEniCS.
Please pardon me if there are any confusing descriptions or format mistakes.
- You can use
area = assemble(1 * dx, domain=mesh)to compute the area of the mesh. - You are correct, it should be
since the domains are marked [1, 2, 3, 4], not [0, 1, 2, 3].Wint = sum ([inner(sigma(v, i), grad(w))*dx(i + 1) for i in range(nphases)])