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 ?
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:
Function
, or did you compute the volume average so that the strains are stored in an array?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:
assemble
as 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.
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:
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
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,
If we take our 2D constitutive relation to be
and apply a unit strain in the x-direction, for example, the volume average of the stress will be:
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,
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.
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
Then the 2D elasticity relation is written as
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.
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>
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.
area = assemble(1 * dx, domain=mesh)
to compute the area of the mesh.Wint = sum ([inner(sigma(v, i), grad(w))*dx(i + 1) for i in range(nphases)])
since the domains are marked [1, 2, 3, 4], not [0, 1, 2, 3].