Two Material Modal Analysis

[lbeardslee1963]

Hi sfepy,

I am trying to conduct modal analysis on a two material build up. I have modified the modal analysis example to fit my needs and have it working well for a single material. At the moment I am importing a file from Gmsh and using that as the model, ultimately I would like to define two regions in Gmsh and assign a different material property to each. I have looked through several examples and have found where to define the different materials, but I have not figured out how to define the different regions of the mesh so that I can assign the materials to them. Once I can figure this out the next question is how do I include the two materials in a single mass and stiffness matrix. Any help would be much appreciated!

[rc]

Hi Luke,

there is a tutorial (section Multipart models) written by @vlukes that shows how to define subdomains using freecad/openscad and gmsh. This manual section shows how to easily define material parameters that are constant per subdomain.

[lbeardslee1963]

Okay, I have read through the tutorials and understand how to define the two regions/subdomains, but how do I add those regions from mesh to omega so that I can define a material property for each? I have pasted what I think is the key piece of code below, I am assuming I need to create two regions, one for group 1, and group 2, then define two materials one for each region/group. Then how do I add these to the two terms defining the mass and stiffness matrix?

omega = domain.create_region('Omega', 'all')
field = Field.from_args('fu', np.float64, 'vector', omega,
                        approx_order=order)

u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')

mtx_d = stiffness_from_youngpoisson(dim, young, poisson)
m = Material('m', D=mtx_d, rho=density)
integral = Integral('i', order=2*order)

t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
eq1 = Equation('stiffness', t1)
eq2 = Equation('mass', t2)
lhs_eqs = Equations([eq1, eq2])

[rc]

Let's say you define two regions Omega1, Omega2. Then try something like:

mat = Material('m', values={'Omega1' : {'rho' : 1, 'D' : ...}, 'Omega2' : {'rho' : 2, 'D', ...}})

[lbeardslee1963]

Alright, so I thought I had it figured out but ran into one more road block. When defining two regions I am able to easily define them by assigning each to their own variable like this: omega = domain.create_region('Omega', 'cells of group 1') omega1 = domain.create_region('Omega1', 'cells of group 2')

That works great and I am able to verify that the correct cells are selected by using this method. The problem I have is that I want to pass the region to the T1 and T2 terms I showed above. So in order to create two regions assigned to one variable I followed the documentation for doing so and tried to define two regions instead like this: omega = domain.create_regions({'Omega': 'cells of group 1', 'Omega1': 'cells of group 2'})

When doing this I get an error saying that the regions do not have an attribute name. I am trying to follow the documentation to figure this problem out but cannot find anything. Could you point me in the right direction? It seems that there might be a difference in how to define the regions in interactive mode.

The documentation says to define regions like this:

`regions = {

: , }`

[rc]

It seems to me that you mix the declarative and imperative ways of region definitions. In the imperative style that you use, you can simply define

t11 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega1, m=m, v=v, u=u)
t12 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega2, m=m, v=v, u=u)
eq1 = Equation('stiffness', t11 + t12)

etc. Or IMO(*) you can use the whole domain region omega in the terms (as you do), while defining the material parameters using omega1, omega2.

(*) If it does not work, let us know.

[lbeardslee1963]

Okay, so I started by trying to use your suggestions which seems like a much cleaner approach to the problem. I set it up like this:

domain = FEDomain('domain', mesh)
omega = domain.create_region('omega', 'all')
omega1 = domain.create_region('omega1', 'cells of group 1')
omega2 = domain.create_region('omega2', 'cells of group 2')
dim = 3
mtx_d = stiffness_from_youngpoisson(dim, 68e9, .36)
mtx_d1 = stiffness_from_youngpoisson(dim, 75e9, .3)
m = Material('m', values={'omega1': {'rho': 2700, 'D': mtx_d}, 'omega2': {'rho': 6000, 'D': mtx_d1}})
order = 2
field = Field.from_args('fu', np.float64, 'vector', omega, approx_order=order)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')
integral = Integral('i', order=order * 2)
t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
eq1 = Equation('stiffness', t1)
eq2 = Equation('mass', t2)
lhs_eqs = Equations([eq1, eq2])
pb = Problem('modal', equations=lhs_eqs)
pb.time_update()
pb.update_materials()

When doing this I get an Index error rho at update_materials(). When I use the other method of defining the region as suggested I don't get any errors and I can solve an Eigen problem using the the mass and stiffness matrix, but the solution is incorrect. I end up getting 12 zeros as the initial values in the Eigen value solution. I set this problem up by substituting in the following piece of code:

t11 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega1, m=m, v=v, u=u)
t12 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega2, m=m1, v=v, u=u)
t21 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega1, m=m, v=v, u=u)
t22 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega2, m=m1, v=v, u=u)
eq1 = Equation('stiffness', t11+t12)
eq2 = Equation('mass', t21+t22)

I hope this provides enough information to find a solution. Thanks!

[rc]

Sorry, I gave you wrong advice: the material definition should be:

m = Material('m', values={'rho': {'omega1': 2700, 'omega2': 6000}, 'D': {'omega1': mtx_d, 'omega2': mtx_d1}})

-> a "transposition" of what I wrote.

As for the zeros, what are your boundary conditions?

[lbeardslee1963]

Ah ha, that works great. I still get 12 zeros, I am using free boundary conditions, so I would expect there to be 6. I think that maybe its the way I have defined the mesh. I defined two regions in Gmsh, but I think its treating them as two free bodies and not creating a common boundary between them. Would this be possible?

[rc]

Yes, that might be the case. You can use script/show_mesh_info.py -d <your mesh file> to see the number of mesh components and/or connect the components (provided their vertices coincide) using script/convert_mesh.py -m <mesh-in> <mesh-out>.

[lbeardslee1963]

Hi, sorry for taking so long to get back. This solved my issues! I just had to use the convert mesh. Thank you for all your help.

[rc]

Thanks for letting us know! Closing...