A simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam

Abstract

The application discloses a simulation method for calculating fluid-structure coupling dynamic response of a dielectric elastomer laminated beam, adopts an absolute node coordinate method to describe the motion of the dielectric elastomer laminated beam, deduces a coordinate transformation matrix between coupled beam units through a multi-unit coupling modeling method of the laminated structure, deduces a generalized elastic force array and a generalized stiffness array of the beam unit based on a Mooney-Rivlin hyperelastic constitutive model and a dielectric elastomer force-electricity constitutive, and obtains the dynamic equation of the dielectric elastomer laminated beam through unit assembly. α The application introduces an immersed boundary-lattice Boltzmann method, simulates a flow field through the lattice Boltzmann method, processes the interaction force between the dielectric elastomer laminated beam and the fluid through the immersed boundary method, establishes a fluid-dielectric elastomer laminated beam system fluid-structure coupling dynamic model, and obtains the displacement and velocity of the dielectric elastomer laminated beam based on the generalized-absolute node coordinate method. The application provides a dynamic model of the dielectric elastomer laminated beam considering the fluid-structure coupling effect, and can predict the motion trajectory of the dielectric elastomer laminated beam under the action of voltage driving and fluid resistance.

Description

Technical Field

[0001] The invention relates to a multi-body system dynamics modeling technology, in particular to a simulation method for calculating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam. Background Art

[0002] The multi-physics coupling problem of flexible bodies is an important issue that needs to be solved and improved in engineering applications. Dielectric elastomers, as flexible intelligent materials, have advantages such as high energy density and fast response speed. They are widely used in the design and manufacture of bionic robotic fish and flapping-wing robots. Flexible structures driven by dielectric elastomers are typical flexible multi-body systems. Their dynamic modeling process needs to consider factors such as the nonlinearity of the dielectric elastomer material, geometric nonlinearity, and fluid-structure coupling effects. Therefore, starting from the theory of dynamic modeling of flexible multi-body systems, considering the real physical field environment, establishing a dynamic model of dielectric elastomer flexible structures that takes into account fluid-structure coupling effects, and studying their multi-field coupling dynamic characteristics is of great significance.

[0003] In "Locomotion of a flapping flexible plate," Hua et al. used the immersed boundary-lattice Boltzmann method-finite element method to study the motion of a flapping flexible plate in a viscous, incompressible stationary fluid, clarifying the mechanism of plate dynamics. However, the flexible plate was described using a finite element model, which is not accurate enough for describing the motion of a flexible plate undergoing large deformations. In "Nonlinear dynamic analysis of anisotropic bimorphdielectric elastomer actuator for soft fish robots," Moumita et al. established a nonlinear dynamic model of an anisotropic convex conical dielectric elastomer actuator based on the Gent hyperelastic model and viscoelastic relaxation model. They studied the dynamic response of the dielectric elastomer actuator at different taper ratios and temperatures, but did not consider the effects of underwater fluid-structure interaction. Summary of the Invention

[0004] The present invention proposes a simulation method for calculating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam.

[0005] The technical solution for achieving the purpose of the present invention is: a simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam, characterized by comprising the following steps:

[0006] Step 1: Establish a physical model of a dielectric elastomer laminated beam. The dielectric elastomer laminated beam has a three-layer structure, wherein the upper and lower layers are dielectric elastomer driving layers driven by voltage, and the middle layer is a dielectric elastomer driven layer. Set the geometric parameters, material parameters, grid parameters of the dielectric elastomer laminated beam, and the driving parameters applied to the dielectric elastomer driving layer.

[0007] Step 2: Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam. The coordinate transformation matrix between the upper beam element and the middle beam element, as well as the coordinate transformation matrix between the lower beam element and the middle beam element, are derived using the multi-element coupling modeling method of the laminated structure. The generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element are derived based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of the dielectric elastomer.

[0008] Step 3: Based on the Newton-Euler equation, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element, and the coordinate transformation matrix, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the dielectric elastomer laminated beam are obtained; combined with the constraints on the dielectric elastomer laminated beam, the dynamic equations of the dielectric elastomer laminated beam are obtained;

[0009] Step 4, using the generalized-α method to solve the dynamic equations of the dielectric elastomer laminated beam, and obtain the displacement, velocity, and acceleration data of the nodes of the upper beam unit, the middle beam unit, and the lower beam unit of the dielectric elastomer laminated beam;

[0010] Step 5: Establish a physical model of the flow field around the dielectric elastomer laminated beam. The flow field is a three-dimensional flow field including an upper boundary, a lower boundary, a left boundary, a right boundary, a front boundary, and a rear boundary. Set the geometric parameters, fluid parameters, and grid parameters of the flow field.

[0011] Step 6: Perform uniform grid division on the flow field, establish the flow field evolution equation based on the lattice Boltzmann method, and calculate the macroscopic density and intermediate velocity of the fluid; determine the continuous kernel distribution function based on the immersed boundary method, interpolate the fluid velocity at the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam in the flow field through the continuous kernel distribution function, and calculate the boundary forces acting on the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam, interpolate the boundary forces of the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam to obtain the boundary forces acting on the flow field grid points, and correct the fluid velocity through the boundary forces of the flow field grid points;

[0012] Step 7: Visualize the obtained node displacement, velocity, and fluid velocity of the dielectric elastomer laminated beam unit to obtain a motion trajectory diagram of the dielectric elastomer laminated beam, a terminal node displacement-time curve diagram, and a flow field velocity cloud diagram;

[0013] A system for simulating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam is provided. The system implements the method for simulating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam to achieve fluid-solid coupling dynamic response simulation of the dielectric elastomer laminated beam.

[0014] A computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the method for simulating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam is implemented to achieve fluid-structure coupling dynamic response simulation of the dielectric elastomer laminated beam.

[0015] A computer-readable storage medium stores a computer program. When the computer program is executed by a processor, the method for simulating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam is implemented to achieve fluid-solid coupling dynamic response simulation of the dielectric elastomer laminated beam.

[0016] Compared with the prior art, the present invention has the following significant advantages: (1) The absolute node coordinate method is used to discretize the dielectric elastomer laminated beam. The mass matrix in the established dynamic equation is a constant matrix, and there is no Coriolis force and centrifugal force, which can accurately describe the deformation and movement of the dielectric elastomer flexible structure. (2) The laminated structure multi-unit coupling modeling method is used to reduce the number of overall unit node coordinates while ensuring the modeling accuracy of the coupled beam unit, thereby improving the calculation efficiency. (3) The established dynamic model takes into account factors such as the fluid-solid coupling effect, the nonlinearity of the dielectric elastomer material, and the geometric nonlinearity caused by the dielectric elastomer flexible structure during the deformation process. By changing the parameter settings, the fluid-solid coupling dynamic response of the dielectric elastomer laminated beam under different voltage drives and fluid resistance can be accurately calculated. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Schematic diagram of the laminated beam model.

[0018] Figure 2 Schematic diagram of electrical actuation of laminated beams.

[0019] Figure 3 Schematic diagram of the boundary constraint conditions of the coupling unit.

[0020] Figure 4(a) shows the overall solution process; Figure 4(b) shows the iterative flow chart of the generalized-α method for solving the equation.

[0021] Figure 5 Schematic diagram of the upper and lower layer driving voltages.

[0022] Figure 6 Schematic diagram of the three-dimensional flow field model.

[0023] Figure 7(a) is a diagram of the end node displacement of the laminated beam under voltage drive and fluid resistance of the embodiment; Figure 7(b) is a diagram of the overall motion change of the laminated beam model of the embodiment; Figure 7(c) is a velocity cloud diagram of the flow field of the embodiment. DETAILED DESCRIPTION

[0024] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0025] The present invention performs fluid-solid coupling dynamic response simulation on a dielectric elastomer laminated beam based on MATLAB, comprising the following steps:

[0026] Step 1: Establish a physical model of a dielectric elastomer laminated beam. The dielectric elastomer laminated beam has a three-layer structure. The upper and lower layers are dielectric elastomer driving layers driven by voltage, and the middle layer is a dielectric elastomer driven layer. Set the geometric parameters, material parameters, grid parameters of the dielectric elastomer laminated beam, and the driving parameters applied to the dielectric elastomer driving layer. Specifically,

[0027] (1) Geometric parameters:

[0028] Each layer of dielectric elastomer is L in length and W in width. The heights of the upper and lower dielectric elastomers are both h1, and the height of the middle driven layer is h2.

[0029] (2) Material parameters:

[0030] density of dielectric elastomer ρ, relative dielectric constant of dielectric elastomer ε;

[0031] (3) Grid parameters

[0032] The number of grids in each layer of dielectric elastomer laminated beam is N=N x ×1, the grid length is L e , the grid width is W e ;

[0033] (4) Driving parameters of the dielectric elastomer driving layer

[0034] The driving voltage of the upper dielectric elastomer is U1, and the driving voltage of the lower dielectric elastomer is U2, and the voltage frequency is f;

[0035] Step 2: Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam. The coordinate transformation matrix between the upper beam element and the middle beam element, as well as between the lower beam element and the middle beam element, is derived using the multi-element coupling modeling method of the laminated structure. Based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relationship of the dielectric elastomer, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element are derived. The specific method is as follows:

[0036] (1) Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam:

[0037] Establish a global coordinate system, with the length direction of the dielectric elastomer laminated beam as the X axis, the width direction as the Y axis, and the thickness direction as the Z axis. Take the midpoint of the left end surface of the middle layer of the dielectric elastomer laminated beam as the origin O, and discretize each layer of the dielectric elastomer laminated beam along the X axis according to the grid parameters as N=N x ×1 unit, then there are N units in total on the three layers. x ×3 units, each layer unit length L e =L / N x ;

[0038] Establish a local coordinate system for the element, with the length direction of the beam element as the x-axis, the width direction as the y-axis, and the thickness direction as the z-axis. The midpoint of the left end face of the beam element is the origin o. Each element has two nodes, located at both ends of the element. The position of any point on the beam element in the global coordinate system is expressed as:

[0039] r=[s1I s2I s3I s4I … s 11 I s 12 I s 13 I s 14 I]q (1)

[0040] Where S is the shape function, and the expression is:

[0041]

[0042] Where ξ=x / l,η=y / l,ζ=z / l, x, y and z are the coordinates of the point in the local coordinate system of the element, I is the third-order unit matrix, l is the length of the beam element, q=(q A ,q B ) is the unit coordinate matrix, q A ,q B are the unit node coordinate matrices respectively, and the expressions are:

[0043] qj=(R,R x ,R y ,R z )j j=A,B (3)

[0044] Where A is the left node of the unit, B is the right node of the unit, and R is the node coordinate in the local coordinate system of the unit. For the left node A, R = (0, 0, 0), and for the right node B, R = (L e ,0,0),R x , R y , R z are the first-order derivatives of R with respect to x, y, and z respectively;

[0045] (2) The coordinate transformation matrix between the upper beam unit and the middle beam unit, as well as between the lower beam unit and the middle beam unit, is derived using the multi-element coupling modeling method of the laminated structure:

[0046] The multi-element coupling modeling method for laminated structures needs to meet the following two basic assumptions:

[0047] 1) Before and after unit deformation, the end faces of the upper beam unit and the lower beam unit are coplanar with the end faces of the middle beam unit; 2) There is no relative sliding between the upper beam unit and the middle beam unit, as well as between the lower beam unit and the middle beam unit, on the contact surface, and the deformation of the two remains consistent;

[0048] Assume that the midpoint of the left end face of the upper beam element is A 1 The midpoint of the left end face of the middle layer beam element is A 2 The midpoint of the coincidence line between the left end faces of the upper beam unit and the middle beam unit is C 1 , C 1 The point belongs to the left end surface of the upper beam unit, and the midpoint of the coincidence line between the middle beam unit and the left end surface of the upper beam unit is C 2 , C 2 The point belongs to the left end surface of the middle layer beam unit; the midpoint of the right end surface of the upper layer beam unit is B 1 The midpoint of the right end face of the middle layer beam element is B 2 The midpoint of the coincidence line between the upper beam unit and the middle beam unit is D 1 , D 1 The point belongs to the right end surface of the upper beam unit. The midpoint of the coincidence line between the middle beam unit and the right end surface of the upper beam unit is D 2 , D 2 The point belongs to the right end surface of the middle flexible follower layer beam element.

[0049] Based on the above basic assumptions, there should be A 1 -C 1 -C 2 -A 2 Four points and B 1 -D 1 -D 2 -B 2The four points are always collinear, and the end face C 1 Dot and C 2 Point, D 1 With D 2 The points remain coincident during the deformation process. 1 with C 2 The coordinate constraint equation between two points and D 1 With D 2 Coordinate constraint equation between two points:

[0050]

[0051] Where r C1 、r C2 、r D1 、r D2 C 1 、C 2 、D 1 、D 2 The position of the point in the global coordinate system, C 1 、C 2 、D 1 、D 2 The absolute position slope of the point in the global coordinate system;

[0052] Substituting formula (1) into (4) we can obtain:

[0053]

[0054]

[0055] Right now:

[0056]

[0057]

[0058] Where S is the shape function, The shape function S of the upper beam element is respectively for x, y, z, yz, y 2 ,z 2 The partial derivative of The shape function S of the middle layer beam element is respectively for x, y, z, yz, y 2 ,z 2 The partial derivative of C 1 、C 2 、D 1 、D 2 The matrix obtained by substituting the spatial coordinates of the point in the global coordinate system into the shape function S, q 1 ,q2 are the unit coordinate matrices of the upper beam unit and the middle beam unit respectively;

[0059] Writing formulas (7) and (8) in matrix form yields:

[0060]

[0061] After transformation, we get:

[0062]

[0063] Where, T 12 is the absolute node coordinate transformation matrix between the upper beam unit and the middle beam unit. Similarly, the coordinate transformation matrix T between the middle beam unit and the lower beam unit can be obtained. 32 .

[0064] (3) Based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of dielectric elastic bodies, the generalized elastic force matrix, generalized stiffness matrix and mass matrix of the high-order three-dimensional two-node beam element are derived:

[0065] The strain energy density function of the Mooney-Rivlin model is:

[0066]

[0067] Where I1 and I2 are the first and second invariants of the right Cauchy-Green deformation tensor C, respectively. 10 and C 01 are material constants, representing the shear modulus and bulk modulus respectively, k is the incompressible constant, and J is the determinant of the deformation gradient tensor J.

[0068] Without considering the influence of temperature, the free energy density function of the dielectric elastomer can be expressed as:

[0069]

[0070] Where U(λ1, λ2, λ3) is the elastic strain energy density function of the dielectric elastomer. The Mooney-Rivlin hyperelastic constitutive model can be used, and then:

[0071]

[0072] is the electric field energy density function, expressed as:

[0073]

[0074] Where, is the nominal electric displacement of the dielectric elastomer, λ1, λ2, λ3 are the stretching ratios of the dielectric elastomer in the length, width, and height directions, respectively, and λ1λ2λ3=1, ε is the dielectric constant of the dielectric elastomer;

[0075] When a dielectric elastomer is driven by an electric field, its electromechanical constitutive relation is:

[0076]

[0077] Where, is the nominal electric field of the dielectric elastomer;

[0078] Substituting formula (15) into formula (14) yields:

[0079]

[0080] Combining the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of dielectric elastic body, the total strain energy density function of dielectric elastic body is obtained as follows:

[0081]

[0082] By performing volume integration of Equation (17) on any beam element, the total free energy of the beam element is obtained as:

[0083]

[0084] By varying Equation (18) according to the principle of virtual work, the generalized elastic force matrix of any beam element is obtained as follows:

[0085]

[0086] Differentiating the nodal coordinate matrix q of any beam element using formula (19) yields the generalized stiffness matrix of the beam element:

[0087]

[0088] Where, S z is the partial derivative of the shape function S with respect to z;

[0089] By differentiating Equation (1) with respect to time, the absolute velocity of any point on any beam element is obtained as:

[0090]

[0091] The kinetic energy of any beam element is expressed as:

[0092]

[0093] Where ρ is the density of the dielectric elastomer, M is the mass matrix, and the expression is:

[0094]

[0095] Step 3: Based on the Newton-Euler equation, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element, and the coordinate transformation matrix, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the dielectric elastomer laminated beam are obtained; combined with the constraints on the dielectric elastomer laminated beam, the dynamic equation of the dielectric elastomer laminated beam is obtained. The specific method is:

[0096] The coordinate transformation matrix is ​​used to couple the upper and lower beam elements and the middle beam element together, and the coupled generalized stiffness matrix is ​​obtained as follows:

[0097]

[0098] Where, is the element generalized stiffness matrix of the middle layer, are the generalized stiffness matrices of the upper and lower beam elements, T 12 is the coordinate transformation matrix between the upper beam element and the middle beam element, T 32 is the coordinate transformation matrix between the lower layer beam element and the middle layer beam element;

[0099] The coupled generalized elastic force matrix is:

[0100]

[0101] Where, is the generalized elastic force matrix of the middle layer, are the element generalized elastic force matrices of the upper beam element and the lower beam element respectively;

[0102] The coupled mass matrix is:

[0103]

[0104] Where, is the unit mass matrix of the middle layer, are the element mass matrices of the upper beam element and the lower beam element respectively;

[0105] All coupled mass matrices and generalized elastic force matrices are assembled into the overall mass matrix and generalized elastic force matrix in node order. Considering the constraints of the laminated beam, the Lagrange multiplier method is introduced into the constraint equations to obtain the dynamic equations of the dielectric elastomer laminated beam:

[0106]

[0107] Where M z is the overall mass matrix of the dielectric elastomer laminated beam, Qz is the overall generalized elastic force array of the dielectric elastomer laminated beam, is the generalized acceleration array of the dielectric elastomer laminated beam, λ is the Lagrange multiplier, Φ(q,t)=0 is the constraint, Φ q is the Jacobian matrix of the constraint;

[0108] Step 4: Use the generalized-α method to solve the dynamic equations of the dielectric elastomer laminated beam to obtain the displacement, velocity, and acceleration data of the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam. The specific method is as follows:

[0109] Parameter initialization in the generalized-α method:

[0110]

[0111] Where P is the spectral radius and η is the convergence error;

[0112] Formula (27) can be written in iterative form as:

[0113]

[0114] In iterative calculation, the iterations of step n+1 and step n satisfy:

[0115]

[0116] The error between iterations is:

[0117]

[0118] Displacement, velocity, and acceleration array updates:

[0119]

[0120] Judgment||R n+1 Is || less than the limit error? If so, execute the following formula and return to formula (30) to calculate the next incremental step

[0121]

[0122] If it is not less than the limited error, return to formula (31) and iterate again until the limited error is met, and then return to formula (30) to calculate the next incremental step.

[0123] Step 5: Set the geometric parameters, fluid parameters, and grid parameters of the flow field, specifically:

[0124] (1) Geometric parameters: The flow field length is L f , width W f , height H f ;

[0125] (2) Fluid parameters: Reynolds number is Re, characteristic velocity is υ;

[0126] (3) Grid parameters: number of grids N in the length direction of the flow field l , the number of grids in the width direction of the flow field N w , the number of grids in the height direction of the flow field N h ;

[0127] Step 6: Perform uniform grid division on the flow field, establish the flow field evolution equation based on the lattice Boltzmann method, and calculate the macroscopic density and intermediate velocity of the fluid; determine the continuous kernel distribution function based on the immersed boundary method, and interpolate the fluid velocity at the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam in the flow field through the continuous kernel distribution function, and calculate the boundary forces acting on the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam. Interpolate the boundary forces of the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam to obtain the boundary forces acting on the flow field grid points, and correct the fluid velocity through the boundary forces of the flow field grid points. The specific method is as follows:

[0128] (1) Uniform grid division of the three-dimensional flow field: Divide the flow field into N grids along the length direction l grids, divided into N along the width direction w grids, divided into N along the height direction h grids, and the grid length l f =L f / N l =W f / N w =H f / N h , the total number of grids obtained by flow field grid division is N=N l *N w *N h , the total number of grid points is N d =(N l +1)*(N w +1)*(N h +1);

[0129] (2) Based on the lattice Boltzmann equation, the macroscopic density and intermediate velocity of the fluid are calculated:

[0130] In the lattice Boltzmann method, the lattice Boltzmann equation is:

[0131]

[0132] Where, f α is the distribution function, x is the fluid spatial position, t is the time, e αis the discrete velocity, α is the number of parameters, for three-dimensional flow field, α ranges from 0 to 18, τ is the relaxation time, δ t is the time step, is the equilibrium distribution function, F α is the external force term;

[0133] The equilibrium distribution function expression is:

[0134]

[0135] Where u is the fluid velocity, c s is the lattice sound speed, ω α is the weight coefficient, ω α and e α The expressions are:

[0136]

[0137]

[0138] External force F α The expression is:

[0139]

[0140] Where f is the external force on the fluid;

[0141] Macroscopic density of the fluid ρ f From the distribution function we get:

[0142] ρ f =∑f α (39)

[0143] The intermediate velocity u of the fluid * From the fluid macroscopic density, discrete velocity and distribution function we get:

[0144]

[0145] (3) Based on the immersed boundary method, the continuous kernel distribution function is determined. The fluid velocity at the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam in the flow field is obtained by interpolation of the continuous kernel distribution function, and the boundary forces acting on the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam are calculated. The boundary forces acting on the grid points of the flow field are interpolated from the boundary forces of the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam. The fluid velocity is obtained by correcting the boundary forces of the flow field grid points:

[0146] In the immersed boundary method, the fluid velocity at the nodes of the upper beam element, middle beam element and lower beam element of the dielectric elastomer laminated beam can be calculated by the fluid intermediate velocity u* Interpolation yields:

[0147] U * =∑u * δ h (41)

[0148] Where, δ h is the continuous kernel distribution function, and its expression is:

[0149]

[0150] Where r is equal to the flow field grid length l f ;

[0151] The boundary forces acting on the nodes of the upper, middle, and lower beam elements of the dielectric elastomer laminated beam can be expressed as:

[0152]

[0153] Where U d is the velocity of the nodes of the upper beam element, middle beam element and lower beam element of the dielectric elastomer laminated beam;

[0154] The boundary force acting on the fluid grid point is:

[0155] f=∑Fδ h Δs (44)

[0156] Where Δs is the distance between grid points;

[0157] The corrected fluid velocity is obtained by the fluid grid point boundary force and the fluid intermediate velocity:

[0158]

[0159] Step 7: Visualize the node displacements, velocities, and fluid velocities of the upper, middle, and lower layers of the dielectric elastomer laminated beam. The specific method is as follows:

[0160] The node displacements, velocities, and fluid velocities of the upper, middle, and lower beam units of the dielectric elastomer laminated beam are visualized in the MATLAB APP DESIGNER to obtain the motion trajectory diagram of the dielectric elastomer laminated beam, the end node displacement-time curve diagram, and the flow field velocity cloud diagram.

[0161] Example

[0162] In order to verify the effectiveness of the solution of the present invention, the following examples are carried out for verification.

[0163] In this embodiment, the parameters shown in Table 1 are used, and the driving voltage is set to a sinusoidal voltage with a voltage amplitude of 5 kV and a frequency of 5 Hz.

[0164] Table 1 System parameter settings used in this embodiment

[0165]

[0166] The present invention uses the absolute nodal coordinate method and the immersed boundary-lattice Boltzmann method to perform fluid-solid coupling dynamic modeling on the dielectric elastomer laminated beam, and writes a program in MATLAB based on the dynamic equation of the dielectric elastomer laminated beam and the fluid-solid coupling algorithm. Its interface is shown in Figures 7(a)-(c). After entering the corresponding parameters and clicking the run button, the displacement-time curve of the end node of the dielectric elastomer laminated beam can be obtained. The overall motion trajectory change of the dielectric elastomer laminated beam over a period of time can be output, and the velocity change cloud map of the flow field can be intuitively displayed.

[0167] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0168] The above-described embodiments merely represent several implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present application. It should be noted that a person of ordinary skill in the art may make various modifications and improvements without departing from the spirit of the present application, and these modifications and improvements fall within the scope of protection of the present application. Therefore, the scope of protection of the present application shall be determined by the appended claims.

Claims

1. A simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam, characterized in that: The following steps are involved: Step 1: Establish a physical model of a dielectric elastomer laminated beam. The dielectric elastomer laminated beam has a three-layer structure, wherein the upper and lower layers are dielectric elastomer driving layers driven by voltage, and the middle layer is a dielectric elastomer driven layer. Set the geometric parameters, material parameters, grid parameters of the dielectric elastomer laminated beam, and the driving parameters applied to the dielectric elastomer driving layer. Step 2: Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam. The coordinate transformation matrix between the upper beam element and the middle beam element, as well as the coordinate transformation matrix between the lower beam element and the middle beam element, are derived using the multi-element coupling modeling method of the laminated structure. The generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element are derived based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of the dielectric elastomer. Step 3: Based on the Newton-Euler equation, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element, and the coordinate transformation matrix, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the dielectric elastomer laminated beam are obtained; combined with the constraints on the dielectric elastomer laminated beam, the dynamic equations of the dielectric elastomer laminated beam are obtained; Step 4, using the generalized-α method to solve the dynamic equations of the dielectric elastomer laminated beam, and obtain the displacement, velocity, and acceleration data of the nodes of the upper beam unit, the middle beam unit, and the lower beam unit of the dielectric elastomer laminated beam; Step 5: Establish a physical model of the flow field around the dielectric elastomer laminated beam. The flow field is a three-dimensional flow field including an upper boundary, a lower boundary, a left boundary, a right boundary, a front boundary, and a rear boundary. Set the geometric parameters, fluid parameters, and grid parameters of the flow field. Step 6: Perform uniform grid division on the flow field, establish the flow field evolution equation based on the lattice Boltzmann method, and calculate the macroscopic density and intermediate velocity of the fluid; determine the continuous kernel distribution function based on the immersed boundary method, interpolate the fluid velocity at the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam in the flow field through the continuous kernel distribution function, and calculate the boundary forces acting on the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam, interpolate the boundary forces of the nodes of the upper beam unit, middle beam unit and lower beam unit of the dielectric elastomer laminated beam to obtain the boundary forces acting on the flow field grid points, and correct the fluid velocity through the boundary forces of the flow field grid points; Step 7: Visualize the obtained dielectric elastomer laminated beam unit node displacement, velocity, and fluid velocity to obtain a motion trajectory diagram of the dielectric elastomer laminated beam, a terminal node displacement-time curve diagram, and a flow field velocity cloud diagram.

2. The simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 1: Set the geometric parameters, material parameters, grid parameters of the dielectric elastomer laminated beam and the driving parameters applied to the dielectric elastomer driving layer, specifically: (1) Geometric parameters: Each layer of dielectric elastomer is L in length and W in width. The heights of the upper and lower dielectric elastomers are both h1, and the height of the middle driven layer is h2. (2) Material parameters: density of dielectric elastomer ρ, relative dielectric constant of dielectric elastomer ε; (3) Grid parameters The number of grids in each layer of dielectric elastomer laminated beam is N=N x ×1, the grid length is L e , the grid width is W e ; (4) Driving parameters of the dielectric elastomer driving layer The driving voltage of the upper dielectric elastomer is U1, and the driving voltage of the lower dielectric elastomer is U2, and the voltage frequency is f.

3. The simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 2: Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam. The coordinate transformation matrix between the upper beam element and the middle beam element, as well as between the lower beam element and the middle beam element, is derived using the multi-element coupling modeling method of the laminated structure. Based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relationship of the dielectric elastomer, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element are derived. The specific method is as follows: (1) Based on the absolute nodal coordinate method, a high-order three-dimensional two-node beam element is selected to discretize the dielectric elastomer laminated beam: Establish a global coordinate system, with the length direction of the dielectric elastomer laminated beam as the X axis, the width direction as the Y axis, and the thickness direction as the Z axis. Take the midpoint of the left end surface of the middle layer of the dielectric elastomer laminated beam as the origin O, and discretize each layer of the dielectric elastomer laminated beam along the X axis according to the grid parameters as N=N x ×1 unit, then there are N units in total on the three layers. x ×3 units, each layer unit length L e =L / N x ; Establish a local coordinate system for the element, with the length direction of the beam element as the x-axis, the width direction as the y-axis, and the thickness direction as the z-axis. The midpoint of the left end face of the beam element is the origin o. Each element has two nodes, located at both ends of the element. The position of any point on the beam element in the global coordinate system is expressed as: r=[s1Is2Is3Is4I…s 11 Is 12 Is 13 Is 14 I] q (1) Where S is the shape function, and the expression is: Where ξ=x / l,η=y / l,ζ=z / l, x, y and z are the coordinates of the node in the local coordinate system of the element, I is the third-order unit matrix, l is the length of the beam element, q=(q A ,q B ) is the unit coordinate matrix, q A ,q B are the unit node coordinate matrices respectively, and the expressions are: q j =(R,R x ,R y ,R z ) j j=A,B (3) Where A is the left node of the unit, B is the right node of the unit, and R is the node coordinate in the local coordinate system of the unit. For the left node A, R = (0, 0, 0), and for the right node B, R = (L e ,0,0),R x , R y , R z are the first-order derivatives of R with respect to x, y, and z respectively; (2) The coordinate transformation matrix between the upper beam unit and the middle beam unit, as well as the coordinate transformation matrix between the lower beam unit and the middle beam unit, are derived using the multi-element coupling modeling method of the laminated structure: The multi-element coupling modeling method for laminated structures needs to meet the following two basic assumptions: 1) Before and after unit deformation, the end faces of the upper beam unit and the lower beam unit are coplanar with the end faces of the middle beam unit; 2) There is no relative sliding between the upper beam unit and the middle beam unit, as well as between the lower beam unit and the middle beam unit, on the contact surface, and the deformation of the two remains consistent; Assume that the midpoint of the left end face of the upper beam element is A 1 The midpoint of the left end face of the middle layer beam element is A 2 The midpoint of the coincidence line between the left end faces of the upper beam unit and the middle beam unit is C 1 , C 1 The point belongs to the left end surface of the upper beam unit, and the midpoint of the coincidence line between the middle beam unit and the left end surface of the upper beam unit is C 2 , C 2 The point belongs to the left end surface of the middle layer beam unit; the midpoint of the right end surface of the upper layer beam unit is B 1 The midpoint of the right end face of the middle layer beam element is B 2 The midpoint of the coincidence line between the upper beam unit and the middle beam unit is D 1 , D 1 The point belongs to the right end surface of the upper beam unit. The midpoint of the coincidence line between the middle beam unit and the right end surface of the upper beam unit is D 2 , D 2 The point belongs to the right end surface of the middle flexible follower layer beam element; Based on the above basic assumptions, there should be A 1 -C 1 -C 2 -A 2 Four points and B 1 -D 1 -D 2 -B 2 The four points are always collinear, and the end face C 1 Dot and C 2 Point, D 1 With D 2 The points remain coincident during the deformation process, from which we can conclude that C 1 with C 2 The coordinate constraint equation between two points and D 1 With D 2 Coordinate constraint equation between two points: Where, C 1 、C 2 、D 1 、D 2 The position of the point in the global coordinate system, C 1 、C 2 、D 1 、D 2 The absolute position slope of the point in the global coordinate system; Substituting formula (1) into (4) we can obtain: Right now: Where S is the shape function, The shape function S of the upper beam element is respectively for x, y, z, yz, y 2 ,z 2 The partial derivative of The shape function S of the middle layer beam element is respectively for x, y, z, yz, y 2 ,z 2 The partial derivative of C 1 、C 2 、D 1 、D 2 The matrix obtained by substituting the spatial coordinates of the point in the global coordinate system into the shape function S, q 1 ,q 2 are the unit coordinate matrices of the upper beam unit and the middle beam unit respectively; Writing formulas (7) and (8) in matrix form yields: After transformation, we get: Where, T 12 is the absolute node coordinate transformation matrix between the upper beam unit and the middle beam unit. Similarly, the coordinate transformation matrix T between the middle beam unit and the lower beam unit is obtained. 32 ; (3) Based on the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of dielectric elastic bodies, the generalized elastic force matrix, generalized stiffness matrix and mass matrix of the high-order three-dimensional two-node beam element are derived: The strain energy density function of the Mooney-Rivlin hyperelastic constitutive model is: Where I1 and I2 are the first and second invariants of the right Cauchy-Green deformation tensor C, respectively. 10 and C 01 are material constants, representing the shear modulus and bulk modulus respectively, k is the incompressibility constant, and J is the determinant of the deformation gradient tensor J; Without considering the influence of temperature, the free energy density function of the dielectric elastomer is expressed as: Where U(λ1, λ2, λ3) is the elastic strain energy density function of the dielectric elastomer. The Mooney-Rivlin hyperelastic constitutive model can be used, and then: is the electric field energy density function, expressed as: Where, is the nominal electric displacement of the dielectric elastomer, λ1, λ2, λ3 are the stretching ratios of the dielectric elastomer in the length, width, and height directions, respectively, and λ1λ2λ3=1, ε is the dielectric constant of the dielectric elastomer; When a dielectric elastomer is driven by an electric field, its electromechanical constitutive relation is: Where, is the nominal electric field of the dielectric elastomer; Substituting formula (15) into formula (14) yields: Combining the Mooney-Rivlin hyperelastic constitutive model and the electromechanical constitutive relation of dielectric elastic body, the total free energy density function of dielectric elastic body is obtained as follows: By performing volume integration of Equation (17) on any beam element, the total free energy of the beam element is obtained as: By varying Equation (18) according to the principle of virtual work, the generalized elastic force matrix of any beam element is obtained as follows: Differentiating the coordinate matrix q of any beam element using formula (19) yields the generalized stiffness matrix of the beam element: Where, S z is the partial derivative of the shape function S with respect to z; By differentiating Equation (1) with respect to time, the absolute velocity of any point on any beam element is obtained as: The kinetic energy of any beam element is expressed as: Where ρ is the density of the dielectric elastomer, M is the mass matrix, and the expression is: M=∫ V ρS T SdV (23)。 4. The simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 3: Based on the Newton-Euler equation, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the high-order three-dimensional two-node beam element, and the coordinate transformation matrix, the generalized elastic force matrix, generalized stiffness matrix, and mass matrix of the dielectric elastomer laminated beam are obtained; combined with the constraints on the dielectric elastomer laminated beam, the dynamic equation of the dielectric elastomer laminated beam is obtained. The specific method is: The coordinate transformation matrix is ​​used to couple the upper and lower beam elements and the middle beam element together, and the coupled generalized stiffness matrix is ​​obtained as follows: Where, is the element generalized stiffness matrix of the middle layer, are the generalized stiffness matrices of the upper and lower beam elements, T 12 is the coordinate transformation matrix between the upper beam element and the middle beam element, T 32 is the coordinate transformation matrix between the lower layer beam element and the middle layer beam element; The coupled generalized elastic force matrix is: Where, is the generalized elastic force matrix of the middle layer, are the element generalized elastic force matrices of the upper beam element and the lower beam element respectively; The coupled mass matrix is: Where, is the unit mass matrix of the middle layer, are the element mass matrices of the upper beam element and the lower beam element respectively; All coupled mass matrices and generalized elastic force matrices are assembled into the overall mass matrix and generalized elastic force matrix in node order. Considering the constraints of the laminated beam, the Lagrange multiplier method is introduced into the constraint equations to obtain the dynamic equations of the dielectric elastomer laminated beam: Where M z is the overall mass matrix of the dielectric elastomer laminated beam, Q z is the overall generalized elastic force array of the dielectric elastomer laminated beam, is the generalized acceleration array of the dielectric elastomer laminated beam, λ is the Lagrange multiplier, Φ(q,t)=0 is the constraint, Φ q is the Jacobian matrix of the constraints.

5. The simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 4: Use the generalized-α method to solve the dynamic equations of the dielectric elastomer laminated beam to obtain the displacement, velocity, and acceleration data of the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam. The specific method is as follows: Parameter initialization in the generalized-α method: Where P is the spectral radius and η is the convergence error; Formula (27) can be written in iterative form as: In iterative calculation, the iterations of step n+1 and step n satisfy: The error between iterations is: Displacement, velocity, and acceleration array updates: Judgment||R n+1 Is || less than the convergence error? If so, execute the following formula and return to formula (30) to calculate the next incremental step. If it is not less than the convergence error, return to formula (31) and iterate again until the convergence error is satisfied, and then return to formula (30) to calculate the next incremental step.

6. The simulation method for calculating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 5: Set the geometric parameters, fluid parameters, and grid parameters of the flow field, specifically: (1) Geometric parameters: The flow field length is L f , width W f , height H f ; (2) Fluid parameters: Reynolds number is Re, characteristic velocity is υ; (3) Grid parameters: number of grids N in the length direction of the flow field l , the number of grids in the width direction of the flow field N w , the number of grids in the height direction of the flow field N h .

7. The simulation method for calculating the fluid-structure interaction dynamic response of a dielectric elastomer laminated beam according to claim 1, characterized in that: Step 6: Perform uniform grid division on the flow field, establish the flow field evolution equation based on the lattice Boltzmann method, and calculate the macroscopic density and intermediate velocity of the fluid; determine the continuous kernel distribution function based on the immersed boundary method, and interpolate the fluid velocity at the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam in the flow field through the continuous kernel distribution function, and calculate the boundary forces acting on the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam. Interpolate the boundary forces of the nodes of the upper beam unit, middle beam unit, and lower beam unit of the dielectric elastomer laminated beam to obtain the boundary forces acting on the flow field grid points, and correct the fluid velocity through the boundary forces of the flow field grid points. The specific method is as follows: (1) Uniform grid division of the three-dimensional flow field: Divide the flow field into N grids along the length direction l grids, divided into N along the width direction w grids, divided into N along the height direction h grids, and the grid length l f =L f / N l =W f / N w =H f / N h , the total number of grids obtained by flow field grid division is N=N l *N w *N h , the total number of grid points is N d =(N l +1)*(N w +1)*(N h +1); (2) Based on the lattice Boltzmann equation, the macroscopic density and intermediate velocity of the fluid are calculated: In the lattice Boltzmann method, the lattice Boltzmann equation is: Where, f α is the distribution function, x is the fluid spatial position, t is the time, e α is the discrete velocity, α is the number of parameters, for three-dimensional flow field, α ranges from 0 to 18, τ is the relaxation time, δ t is the time step, is the equilibrium distribution function, F α is the external force term; The equilibrium distribution function expression is: Where u is the fluid velocity, c s is the lattice sound speed, ω α is the weight coefficient, ω α and e α The expressions are: External force F α The expression is: Where f is the external force on the fluid; Macroscopic density of the fluid ρ f From the distribution function we get: r f =∑f α (39) The intermediate velocity u of the fluid * From the fluid macroscopic density, discrete velocity and distribution function we get: (3) Based on the immersed boundary method, the continuous kernel distribution function is determined. The fluid velocity at the nodes of the upper beam unit, the middle beam unit and the lower beam unit of the dielectric elastomer laminated beam in the flow field is obtained by interpolation of the continuous kernel distribution function, and the boundary forces acting on the nodes of the upper beam unit, the middle beam unit and the lower beam unit of the dielectric elastomer laminated beam are calculated. The boundary forces acting on the flow field grid points are interpolated from the boundary forces of the nodes of the upper beam unit, the middle beam unit and the lower beam unit of the dielectric elastomer laminated beam. The fluid velocity is obtained by correcting the boundary forces of the flow field grid points. Specifically, In the immersed boundary method, the fluid velocity at the nodes of the upper beam element, middle beam element and lower beam element of the dielectric elastomer laminated beam is calculated by the fluid intermediate velocity u * Interpolation yields: U * =∑u * d h (41) Where, δ h is the continuous kernel distribution function, and its expression is: Where r is equal to the flow field grid length l f ; The boundary forces acting on the nodes of the upper, middle, and lower beam elements of the dielectric elastomer laminated beam are expressed as: Where U d is the velocity of the nodes of the upper beam element, middle beam element and lower beam element of the dielectric elastomer laminated beam; The boundary force acting on the flow field grid point is: f=∑Fδ h Δs (44) Where Δs is the distance between grid points; The corrected fluid velocity is obtained by the boundary force of the flow field grid point and the intermediate velocity of the fluid:

8. The computational dielectric elastomer laminated beam fluid-structure coupling dynamic response simulation system according to claim 1, characterized in that: The fluid-solid coupling dynamic response simulation method of the dielectric elastomer laminated beam according to any one of claims 1 to 7 is implemented to realize the fluid-solid coupling dynamic response simulation of the dielectric elastomer laminated beam.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the method for simulating the fluid-solid coupling dynamic response of a dielectric elastomer laminated beam according to any one of claims 1 to 7 is implemented to achieve simulation of the fluid-solid coupling dynamic response of the dielectric elastomer laminated beam.

10. A computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, the method for simulating the fluid-structure coupling dynamic response of a dielectric elastomer laminated beam according to any one of claims 1 to 7 is implemented to achieve the fluid-structure coupling dynamic response simulation of the dielectric elastomer laminated beam.

Images